Question

I have some arbitrary function, f(x), that has no period. I wish to find the magnitude of 1 frequency component in f(x). Could someone please recommend a fast way to do this?

So far, my ideas involve:

  1. Taking the inner product of f(x) with the desired frequency e.g. <f(x),cos(kx)>. This would essentially be finding the Fourier Coefficient at a given frequency. To find the fourier coeffieicients, I am meant to do it over one period; however, since the period is large, I suppose I would just do it over less than one period and hope that it's accurate enough...

  2. Doing an FFT of the signal, and then simply finding the frequency component I want from the FFT.

  3. Somehow perhaps doing something with a narrow band-pass filter.

Edit: Thank you for the suggestion.

I found the following algorithm written in Python (which is what I need) and slightly adjusted it to suit.

def Goetrzel(x, target_frequency, sample_rate):
    s_prev = 0
    s_prev2 = 0
    normalized_frequency = target_frequency / sample_rate

    wr = np.cos(2.0 * np.pi * normalized_frequency)
    wi = np.sin(2.0 * np.pi * normalized_frequency)

    coeff = 2.0 * wr
    for sample in x:
        s = sample + coeff * s_prev - s_prev2
        s_prev2 = s_prev
        s_prev = s

    XKreal = s_prev * wr - s_prev2
    XKimag = s_prev * wi

    XK = (XKreal + 1j*XKimag) / (len(x)/2.)

    #power = s_prev2 * s_prev2 + s_prev * s_prev - coeff * s_prev * s_prev2 ;
    return abs(XK), np.angle(XK)*180./PI

I test it with:

h = 0.01
z = np.arange(0., 300., h)
f = np.cos(2.*PI*3.*z)
print Goetrzel(f, 3., 1./h)

If I make h = 0.0007, the answer calculated by the Goetrzel algorithm seems to significantly change (increases massively). Do you have any idea why, or how I can make sure I pick the correct step?

So changing h from 0.06 to 0.07 changes the phase from -1.388e-11° to 21.6°.

Was it helpful?

Solution

Explore the Görtzel algorithm. It was designed to detect the presence of one or few (16) frequency components in a signal.

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top