Home > LSCF, Statistics > Lines, Sines, and Curve Fitting 15 – Periodic Noise

## Lines, Sines, and Curve Fitting 15 – Periodic Noise

2011 February 16

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