[R] Fitting a sine wave using solver

[Thomas Petzoldt]

Ben Zuckerberg schrieb:
> Greetings,
> 
> I have several sets of oscillation data and would like to estimate the 
> parameters of a sine function to each set (and hopefully automate 
> this).  A colleague provided an excel sheet that uses solver to minimize 
> the RSS after fitting the sine function to each data set, but this 
> cumbersome and difficult to automate.  Is there a method in R for 
> fitting a given sine function to a supplied data using maximum 
> likelihood estimation (or minimizing the RSS).  Thanks in advance.
> 

Hi Ben,

why not using FFT? See one of the examples below.

HTH Thomas P.


### Some periodic data ----------------------------
n   <- 360
ord <- 180
x   <- 0:(n-1)
y <- sin((50 + x) * 2 * pi / n) + rnorm(n) * 0.1
plot(x, y,xlim=c(-10,370))

### example 1 -------------------------------------
## Fast Fourier Transform
pf <- fft(y)

## Plot highest possible order
pf[(n/2):n] <- 0
yy <- 2*Re(fft(pf, inverse=TRUE)/n) - Re(pf[1])/n
lines(x, yy, col="gray", lwd=1)

### example 2 -------------------------------------
## Fast Fourier Transform
pf <- fft(y)

## another, lower order
pf[4:n] <- 0
yy <- 2*Re(fft(pf, inverse=TRUE)/n) - Re(pf[1])/n
lines(x, yy, col="blue", lwd=2)


### example 3 -------------------------------------
## Fast Fourier Transform
pf <- fft(y)

## synthesize it the traditional way
n  <- length(y)
a0 <- Re(pf[1])/n
pf <- pf[-1]
a  <-  2*Re(pf)/n
b  <- -2*Im(pf)/n

harmonic <- function(x, a0, a, b, n, ord) {
   k <- (2 * pi * x/n) %*% t(1:ord)
   y <-    a0 + cos(k) %*% a[1:ord] +
                sin(k) %*% b[1:ord]
   y
}

lines(x, harmonic(x, a0, a, b, n, ord=2), col="red", lty=2, lwd=2)