1
$\begingroup$

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?

Improve this question
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$ Commented May 17, 2016 at 9:27
  • $\begingroup$ en.wikipedia.org/wiki/Goertzel_algorithm $\endgroup$ Commented May 17, 2016 at 13:50
  • $\begingroup$ @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? $\endgroup$ Commented May 17, 2016 at 14:10
  • $\begingroup$ 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$. $\endgroup$ Commented May 17, 2016 at 16:01
  • $\begingroup$ 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. $\endgroup$ Commented May 17, 2016 at 16:01

2 Answers 2

Reset to default
0
$\begingroup$

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.

enter image description here

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")
Improve this answer
$\endgroup$
6
  • $\begingroup$ 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... $\endgroup$ Commented May 18, 2016 at 9:58
  • $\begingroup$ 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. $\endgroup$ Commented May 18, 2016 at 12:09
  • 1
    $\begingroup$ 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. $\endgroup$ Commented May 19, 2016 at 11:27
  • $\begingroup$ 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? $\endgroup$ Commented May 19, 2016 at 12:19
  • $\begingroup$ 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. $\endgroup$ Commented May 19, 2016 at 12:27
0
$\begingroup$

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$

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.