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


---

The code is

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)
 }