fft power spectrum woes

Question

I'm having trouble getting a frequency spectrum out of a fourier transform... I have some data:

That I have mean-centered, and doesn't seem to have too much of a trend...

I plot the fourier transform of it:

And I get something that is not nice....

Here is my code:

def fourier_spectrum(X, sample_freq=1):
    ps = np.abs(np.fft.fft(X))**2
    freqs = np.fft.fftfreq(X.size, sample_freq)
    idx = np.argsort(freqs)

    plt.plot(freqs[idx], ps[idx])

As adapted from code taken from here.

It seems to work for some naive sin wave data:

fourier_spectrum(np.sin(2*np.pi*np.linspace(-10,10,400)), 20./400)

So my questions are: I'm expecting a non-zero-almost-everywhere-spectrum, what am I doing wrong? If I'm not doing anything wrong, what features of my data are causing this? Also, if I'm not doing anything wrong, and fft is just unsuited for my data for some reason, what should I do to extract important frequencies from my data?

Just don't plot the 0 bin of the fft if there's too much "DC" component. Also by eyeball the data looks to have strong frequency components at very low frequencies - may look more "reasonable" by just looking at the 1st 50 or so bins — f5r5e5d
Subtract the mean of the signal before FFT—that removes the big component at “DC” (i.e., bin 0 of the FFT). Also/alternatively, plot the log10(abs(...)) (i.e., dB, or decibels) of the FFT output to get a better sense of what’s going on. — Ahmed Fasih
Also look at statistical spectral estimators like periodogram and related methods in scipy.signal. — Ahmed Fasih

Him Him · Accepted Answer · 2017-03-09T15:20:45

It turns out that I just didn't understand the units of the x-axis in the frequency spectrum, which is Hz. Because my sample spacings were on the order of a second, and my period was on the order of a day, the only units really visible on my frequency spectrum were ~1/s (at the edges) to about ~1/m (near the middle), and anything with a longer period than that was indistinguishable from 0. My misunderstanding stemmed from the graph on this tutorial, where they do conversions so that the x-axis units are in time, as opposed to inverse time. I rewrote my frequency_spectrum plotting function to do the appropriate "zooming" on the resulting graph...

def fourier_spectrum(X, sample_spacing_in_s=1, min_period_in_s=5):
    '''
        X: is our data
        sample_spacing_in_s: is the time spacing between samples
        min_period_in_s: is the minimum period we want to show up in our
            graph... this is handy because if our sample spacing is
            small compared to the periods in our data, then our spikes
            will all cluster near 0 (the infinite period) and we can't
            see them.  E.g. if you want to see periods on the order of
            days, set min_period_in_s=5*60*60 #5 hours
    '''
    ps = np.abs(np.fft.fft(X))**2
    freqs = np.fft.fftfreq(X.size, sample_spacing_in_s)
    idx = np.argsort(freqs)
    plt.plot(freqs[idx], ps[idx])
    plt.xlim(-1./min_period_in_s,1./min_period_in_s) # the x-axis is in Hz

fft power spectrum woes

1 Answers