# Extract parameters of a single harmonic

## Context

I need to extract the parameters of a single harmonic in a signal.

For instance if

$$x(t) = A \sin (\omega_0t + d_0) + B \sin (\omega_2t + d_2)$$ then $f(x, \omega_0) = (A, d_0)$. This can be done offline.

Please consider the fact that $t$ does not cover a full number of periods!

Of course I used a fft() and then extracted the coefficient corresponding to the interesting frequency. The problem is that I compute the whole spectrum when I'm interested in only one value.

## What I tried and read

I first tried to use the definition of the DFT (or to project my signal on a sine/cosine plane) by using:

$$A e^{jd_0} = \frac{2}{N} \sum_t e^{j\omega_0 t} x(t) ,$$ with $N$ the number of samples of $t$.

This is what is proposed here: FFT for a single frequency

This works however only for a full number of periods of the interesting frequency: I can cut my signal to get reasonable results but I'd rather not (and it's not always possible)

What can I do else?

• tried something that I didn't thought of before: force w<sub>0</sub> to be a multiple of pi*fs/N, with fs my sampling rate. However it doesn't always give the same results as the fft. – ochurlaud May 17 '16 at 9:27
• en.wikipedia.org/wiki/Goertzel_algorithm – endolith May 17 '16 at 13:50
• @endolith The link was given in the other answer I linked. It works well but only for frequencies with the form given in the answer I proposed. How to do with an arbitrary frequency like 10.2521 when I don't have enough time sample to get this precision? – ochurlaud May 17 '16 at 14:10
• If $\omega_1$ and $\omega_1$ aren't too close you can try to shift the frequency (multiply by $exp(j2\pi f_ct)$) to move the frequency on to one of the FFT bins - when the sidelobes disappear then it is lined up exactly. Not sure how well this works with 2 sinusoids. Note there is also an ambiguity between the amplitude and phase i.e. $A$ and $d_0=\pi$ is the same as $-A$ and $d_0=0$. – David May 17 '16 at 16:01
• You could try to calculate a Phase Locked Loop with your specific characteristics. This will give you the actual frequency / phase of the incoming signal within a well defined range and also give you a "control signal" of which direction the changes occur towards. – A_A May 17 '16 at 16:01

There should be zero need to use exact periods of the signal.

I have changed your signal model a little to:

$$x(t) = A \cos (\omega_0t + d_0) + B \cos (\omega_2t + d_2)$$

and updated the sum to:

$$\frac{2}{N} \sum_t e^{-j\omega_0 t} x(t)$$

(i.e. the minus sign in front of $j$).

And I get sensible results illustrated below. The $\color{green}{\tt green}$ lines in the bottom plots show the actual values. The black circles are the estimates of $A$ and $d_0$ at each point in time.

R Code Below

#30846

A <- 1
omega_0 <- 2*pi*0.0143298579849
d_0 <- 2*pi*0.39238798235

B <- 1
omega_2 <- 2*pi*0.363209429
d_2 <- 2*pi*0.1835934324

T <- 1000
t <- seq(0,T-1)

x <- A*cos(omega_0*t + d_0) + B*cos(omega_2*t + d_2)

analytic <- function(x)
{
hilb <- c(1, rep(2, T-2), rep(0,T+1))
ana <- fft(hilb*fft(c(x, rep(0,T))), inverse=TRUE)/T/2
return (ana[1:T])
}

z <- analytic(x)
exponential <- exp(-1i*omega_0*t)

coeff <- rep(0,T)
for (time in t)
{
coeff[time] <- 2*sum(exponential[1:time]*x[1:time])/time
}

par(mfrow=c(2,2))
plot(x,type="l", lwd=5)
#lines(Re(z), col="red")
title('Original signal ')

plot(abs(fft(x)), type="l")
title("FFT of signal")

plot(Arg(coeff))
lines(c(0,T), c(d_0,d_0), col="green")
title("d_0 estimate as a funciton of time")

plot(abs(coeff))
lines(c(0,T), c(A,A), col="green");
title("A estimate as a function of time")

• Oh, this seems promising. I'll read that later this week and certainly accept it as the right answer. Can you add some references about your algorithm? (it's name would be enough). Is it Goerzel? It seems so... – ochurlaud May 18 '16 at 9:58
• I have just implemented your algorithm, with the two minor changes I mention in the answer! I changed the signal model slightly to make the generated phase be correct (otherwise it'll be out by $\pi/2$) and I added a minus sign to the exponent in the sum. Yes, it's just the Goertzel algorithm. – Peter K. May 18 '16 at 12:09
• Because I really need a sin/cos, using $exp(+j \omega_0 t)$ is preferable. Else it works perfect. I still don't get why it so perfectly works: I'll read some doc on this. I mark this as the correct answer, thx for the nice example. – ochurlaud May 19 '16 at 11:27
• Well, at the end it is really the same algorithm as I proposed since I can just take the last value (if I'm sure that it converged). So it works only if my sample number is high enough, which is not always the case. I guess that doing better is physically impossible? – ochurlaud May 19 '16 at 12:19
• You'll need to have enough samples and the second frequency will need to be far enough away from the frequency of interest. And as little other noise as possible. – Peter K. May 19 '16 at 12:27

I tried something that I didn't thought of before:

force $\omega_0$ to be a multiple of $2 \pi \frac{fs}{N}$, with $fs$ my sampling rate.

And of course, this is what happens in the fft.

I only ask myself if I can get the same results for $\omega_0\neq2 \pi \frac{k}{N} fs$