Lines, Sines, and Curve Fitting 15 – Periodic Noise
In this post, I am examining some examples of noisy signals.
I start with a simple sine wave, amplitude = 1, period = 40.
Add noise to the amplitude. This noise has a standard deviation of 0.5.
Or add noise to the period – making this “quasi-periodic.” This noise has a standard deviation of 10.
Or just add noise to each point in the series. This noise has a standard deviation of 1.
Make everything noisy.
We can easily detect the period of 40 from the clean gray sine in the previous figure.
1/sp$freq[sp$spec>2]
[1] 40
The Fourier transform zeros right into the 40 period.
But what does the Fourier transform detect in the very noisy signal?
1/sp$freq[sp$spec>2]
[1] 120.000000 60.000000 40.000000 30.000000 24.000000 13.333333
[7] 6.315789 5.714286 3.750000 3.636364 3.428571
By setting a higher threshold, I could have easily extracted just the 40 period, but I wanted to show the additional frequency peaks. Note that 120 is greater than the original source and is a simple multiple of the original. 30 is a simple fraction (3/4) of the source. As is 24 which is 3/5 of the original. Of course, 6.315789 doesn’t fit this pattern.
We can also reduce the number of high frequency ‘hits’ by using a smoothing function, in this case a two-sided moving average of 10.
1/sp$freq[sp$spec>2]
[1] 40.00000 30.00000 17.14286
getQuasiPeriod <- function(To,Tsd,Ao,Asd,No,Nsd,tlen,savg) { w <- NULL T <- NULL A <- NULL wx <- NULL periods = (tlen/To) * 2 # for wiggle room for (i in 1:periods) { Ti = To + rnorm(1,sd=Tsd) idx = Ti-1 T = c(T,rep(Ti,idx)) # set the freq for each point in the period Ai = Ao + rnorm(1,sd=Asd) A = c(A,rep(Ai,idx)) # set the amplitude for point in the each period wx = c(wx,rep((length(wx)+1),idx)) # used to reset the phase each period } x <- 1:length(wx) # get a series to work the function # if set, add noise to each point in the series y = A*sin(((x-wx)/T)*(2 * pi)) + (No + rnorm(length(wx),sd=Nsd)) # smooth the signal by 'savg' y = filter(y, rep(1/savg,savg), sides=2) y = y[!is.na(y)] # remove the NAs from the filtering y = y[1:tlen] # chop the longer series to the requested length return(y) }
-
2011 February 20 at 8:55 amLines, Sines, and Curve Fitting 16: Wavelets « The Whiteboard