1
$\begingroup$

The measure of a given frequency $\omega$ in a signal $x(t)$ is:

$\frac{1}{N}\sum\limits^N_{t=0}x\left(t\right)e^{^{-i \omega t}}$

This is basically an average of the correlation between the signal and a complex sinusoid rotating at $\omega$ frequency over some time range (N).

My goal is to realize this as a IIR filter (the former definition is virtually a FIR filter), so that it can be implemented as a sort of resonator which is fed an input stream of data, from which the frequency over the last N samples can be measured, continuously?

Currently, I've made this working model:

$y\left(n\right)=e^{-i \omega n}x\left(n\right)$

$X\left(n\right)=y\left(n\right)+e^{\frac{-1}{N}}\left(X\left(n-1\right)-y\left(n\right)\right)$

Where $X(n)$ will be the complex result of the correlation over the last N samples, utilizing a moving average IIR filter. While this works, it is quite inefficient in my current, naive implementation:

for each sample:
    t0 = x * sample;
    real = t0 + pole * (real - t0);

    t0 = -y * sample;
    imag = t0 + pole * (imag - t0);

    t0 = x * c1 - y * c2;
    y = x * c2 + y * c1;
    x = t0;

where x and y corresponds to the cosine and the sine of the complex sinusoid, computed using a IIR oscillator filter. My hunch is that it will be possible to implement some standard filter with an adequate damping such that it will be equivalent to the former equation / code? And basically function like a resonator, where the amount of resonance is equal to the average of the correlation over the last N samples.

I've also been looking at the Goertzel algorithm, and possible ways to use that theory, even though it seems to suffer from scalloping/inability to resolve frequencies not quantized to integers between 0 and N?

Sorry for the redundancy, I'm still very new to filters and DSP in general, I'm just hoping someone can point me in the right direction and tell me I'm not grasping for straws :)

Thanks

Edit - appendix:

goertzel algorithm

Scalar sine, cosine, coeff, q0(0), q1(0), q2(0);

sine = sin(omega);
cosine = cos(omega);
coeff = 2.0 * cosine;

for (int t = 0; t < size; t++)
{
    q0 = coeff * q1 - q2 + data[t];
    q2 = q1;
    q1 = q0;
}

Scalar real = (q1 - q2 * cosine) / (size * 0.5);
Scalar imag = (q2 * sine) / (size * 0.5);

goertzel vs. correlation of sinusoid

enter image description here

$\endgroup$
  • $\begingroup$ I believe that it is indeed the Goertzel algorithm that you're looking for. It's not correct that the Goertzel algorithms needs quantized frequencies. Unlike the FFT, it is not restricted to any frequency grid. $\endgroup$ – Matt L. Oct 22 '14 at 21:36
  • $\begingroup$ The Sliding DFT might also be interesting for you. $\endgroup$ – Matt L. Oct 23 '14 at 9:42
  • $\begingroup$ @MattL. Okay, you're right about the quantization (it was an error in my code). I am however still having numerical/quantization troubles using the goertzel versus my original code (which is bruteforce correlation between signal and sinusoid) using single-precision. This happens with very closely-spaced frequencies (~0.5Hz). The difference can be seen in this picture: i.imgur.com/JJJePO3.png Also, I'm unsure how to convert the goertzel algorithm into an iir-filter (the usual implementations emulates fir). I've included the algorithm in the OP. $\endgroup$ – Shaggi Oct 26 '14 at 13:07
  • $\begingroup$ @MattL. Also, if i use the sliding dft, i would have to do windowing in the frequency domain (ie. by convolution), am i right? $\endgroup$ – Shaggi Oct 28 '14 at 18:20
  • $\begingroup$ No, not necessarily. Windowing can be used to reduce spectral leakage, but it's not an integral part of the SDFT algorithm. Have a look at this document. It gives a good overview of the SDFT algorithm compared to Goertzel. $\endgroup$ – Matt L. Oct 29 '14 at 8:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.