require(knitr) require(formatR) require(dplyr) options(width=200) knitr::opts_chunk$set(cache=FALSE,prompt=FALSE,comment=">",message=FALSE,echo=TRUE,warning=FALSE,tidy=TRUE,strip.white=TRUE,size="small", fig.align = "center",fig.show='hold')
It will become increasingly difficult to use software like Excel and SPSS. Perhaps now is a good time to switch to R or Matlab. We do have a spreadsheet example of Standardised Dispersion Anaysis.
nlRtsa
First, download an install package nlRtsa
(under construction)
require(devtools) install_github("FredHasselman/nlRtsa") library(nlRtsa)
The nlRtsa
package contains a nice function in.IT()
that will load a list of packages and install them if they are not present on your system.
in.IT(c("signal","pracma"))
Below is an example of a signal built from sine components (y
) whose relative amplitudes are recovered in the powerspectrum.
# Sawtooth x <- seq(-3.2,3.2, length.out = 256) y <- 2*sin(10*x) - 1*sin(20*x) + (2/3)*sin(30*x) - (1/2)*sin(40*x) + (2/5)*sin(50*x) - (1/4)*sin(60*x) # Plot the sawtooth wave as constructed by the Fourier series above plot(x,y, xlab ='Time (a.u.)', ylab = 'Variable (a.u.)', main ='Sawtooth wave', type = "l") # Perform a Fast Fourier Transform and calculate the Power and Frequency Y <- fft(y) Pyy <- Y*Conj(Y)/256 f <- 1000/256*(0:127) # Plot the power spectrum of the sawtooth wave plot(f[1:50],Pyy[1:50], type="b",xlab='Frequency (a.u.)', ylab ='Power (a.u.)', pch=21, bg='grey60', main = 'Power Spectrum')
Now we do the same for a very noisy signal into which we insert a frequency.
# A time vector t <- pracma::linspace(x1 = 0, x2 = 50, n = 256) # There are three sine components x <- sin(2*pi*t/.1) + sin(2*pi*t/.3) + sin(2*pi*t/.5) # Add random noise! y <- x + 1*randn(size(t)) # Plote the noise. plot(t, y, type = "l", xlab = 'Time (a.u.)', ylab = 'Variable (a.u.)', main = 'A very noisy signal') # Get the frequency domain Y <- fft(y) Pyy <- Y*Conj(Y)/256 f <- 1000/256*(0:127) # Plot the power spectrum of this noise plot(f[1:50],Pyy[1:50], type="b",xlab='Frequency (a.u.)', ylab='Power (a.u.)', pch=21, bg='grey60', main = 'Power Spectrum')
We can use the power spectrum to estimate a self-affinity parameter, or scaling exponent.
ts1.txt
, ts2.txt
, ts3.txt
from blackboard. The text files contain time series.R
using any of the read
functions. TIP Use the function import()
in package rio
.Now we can do some basic data preparations:
log2(length of var)
)scale()
can do this for you, but you could also use mean()
and sd()
.Before a spectral analysis you should remove any linear trends (it cannot deal with nonstationary signals!)
pracma::detrend()
. stats::lm()
or stats::poly()
The function nlRtsa::fd.psd()
will perform the spectral slope fitting procedure.
R
and pressing F2
plot = TRUE
fd.dfa()
and fd.sda()
to estimate the self-affinity parameter and Dimension of the series. pracma::sample_entropy()
)growth.ac( ..., type = "logistic")
in nlRtsa
)Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.