Scipy/Numpy FFT Frequency Analysis

Question

I'm looking for how to turn the frequency axis in a fft (taken via scipy.fftpack.fftfreq) into a frequency in Hertz, rather than bins or fractional bins.

I tried to code below to test out the FFT:

t = scipy.linspace(0,120,4000)
acc = lambda t: 10*scipy.sin(2*pi*2.0*t) + 5*scipy.sin(2*pi*8.0*t) + 2*scipy.random.random(len(t))

signal = acc(t)

FFT = abs(scipy.fft(signal))
FFT = scipy.fftpack.fftshift(FFT)
freqs = scipy.fftpack.fftfreq(signal.size)

pylab.plot(freqs,FFT,'x')
pylab.show()

The sampling rate should be 4000 samples / 120 seconds = 33.34 samples/sec.

The signal has a 2.0 Hz signal, a 8.0 Hz signal, and some random noise.

I take the FFT, grab the frequencies, and plot it. The numbers are pretty nonsensical. If I multiply the frequencies by 33.34 (the sampling frequency), then I get peaks at about 8 Hz and 15 Hz, which seems wrong (also, the frequencies should be a factor of 4 apart, not 2!).

Any thoughts on what I'm doing wrong here?

HYRY · Accepted Answer · 2012-02-27 00:50:59Z

48

I think you don't need to do fftshift(), and you can pass sampling period to fftfreq():

import scipy
import scipy.fftpack
import pylab
from scipy import pi
t = scipy.linspace(0,120,4000)
acc = lambda t: 10*scipy.sin(2*pi*2.0*t) + 5*scipy.sin(2*pi*8.0*t) + 2*scipy.random.random(len(t))

signal = acc(t)

FFT = abs(scipy.fft(signal))
freqs = scipy.fftpack.fftfreq(signal.size, t[1]-t[0])

pylab.subplot(211)
pylab.plot(t, signal)
pylab.subplot(212)
pylab.plot(freqs,20*scipy.log10(FFT),'x')
pylab.show()

from the graph you can see there are two peak at 2Hz and 8Hz.

enter image description here

answered Feb 27, 2012 at 0:50

HYRY

97.8k28 gold badges197 silver badges192 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

Archie1986 Over a year ago

Thank you for such a complete answer. hyry, why did u choose to plot 20*scipy.log10(FFT) instead of FFT?

OldTinfoil Over a year ago

HYRY provided you a plot with the Y axis in the dB scale, and 20log10 provides the correct conversion for a magnitude spectrum.

tom10 · Accepted Answer · 2012-02-27 04:41:07Z

11

scipy.fftpack.fftfreq(n, d) gives you the frequencies directly. If you set d=1/33.34, this will tell you the frequency in Hz for each point of the fft.

answered Feb 27, 2012 at 4:41

tom10

69.5k11 gold badges133 silver badges143 bronze badges

Comments

Oliver Charlesworth · Accepted Answer · 2012-02-26 20:03:06Z

6

The frequency width of each bin is (sampling_freq / num_bins).

~~A more fundamental problem is that your sample rate is not sufficient for your signals of interest. Your sample rate is 8.3 Hz; you need at least 16Hz in order to capture an 8Hz input tone.¹~~

~~_{1. To all the DSP experts; I'm aware that it's actually BW that's relevant, not max frequency. But I'm assuming the OP doesn't want to do undersampled data acquisition.}~~

edited Feb 26, 2012 at 20:03

answered Feb 26, 2012 at 18:53

Oliver Charlesworth

274k34 gold badges591 silver badges687 bronze badges

2 Comments

nathan lachenmyer Over a year ago

I'm using 4000 samples for 120 seconds -- isn't that 33.3 Hz? That should be more than enough for it, and the numbers are still off...

Oliver Charlesworth Over a year ago

@asymptoticdesign: Ah, ok. Your question initially said 1000. Yes, that should be sufficient. Which bin index is the energy appearing in?

Collectives™ on Stack Overflow

Scipy/Numpy FFT Frequency Analysis

3 Answers 3

2 Comments

Comments

2 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

3 Answers 3

2 Comments

Comments

2 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related