python-FFT and PSD calculating

Question

I have a set of data that consists of 36002 items and I want to do FFT and PSD of it to know which frequency it includes and corresponding power density of the frequency.

My code is:

from __future__ import division
import numpy
import matplotlib.pyplot as plt

#read in the pressure p_dot and time t, they are [36002,] vector
nSteps=36002
p_dot=numpy.genfromtxt((r'E:\p_dot.dat'), delimiter=' ')[:,2]
t=numpy.genfromtxt((r'E:\t.dat'), delimiter=' ')[:,0]

T=(t[-1]-t[0])/nSteps # the interval between two data points
N=len(p_dot)//2+1 # FFT is symmetrical, so plot one half
Y=numpy.fft.fft(p_dot) # to compare with Yhann
hann=numpy.hanning(len(p_dot))
Yhann=numpy.fft.fft(hann*p_dot)
fa=1.0/T # scan frequency
X=numpy.linspace(0, fa/2, N, endpoint=True) # Nyquest frequency=fa/2
plt.close()
plt.subplot(2,1,1)
plt.plot(x,y)
plt.subplot(2,1,2)
plt.plot(X, 2.0*abs(Yhann[:N])/N)
plt.tight_layout()
plt.show()

But the results turned out incorrect, at least I think so. Apparently my data is periodical, but the FFT only has a spike at 0 Hz. See the result of first 1000 items.

result 1000 And the result of 36002 result 36002 So what's the problem? Thanks a lot.

Btw, can anyone explain how to use overlap window in python? And when should we use it? Since when I search the solution of my problem, I always see these methods and no idea how to use it. Thanks!

0 Hz means zero frequency, i.e. a constant signal, think DC component in the electrical world. It's not unusual that this is a large component in a measured signal, it depends on what you do measure. You could try and remove the mean from the original signal, and also use detrend functions in for example scipy.signal.detrend. Simplest way to know if you your code is doing the right thing is to give it a fake signal where you know what the outcome should be. — ahed87

bnaecker bnaecker · Accepted Answer · 2017-12-16T20:06:29

The reason for the peak at 0Hz is that your input signal has a non-zero mean. This is the DC component of the signal. There is another, smaller, peak at around 1Hz or so, which looks to be roughly the main frequency of your signal.

Window functions are usually applied when one computes power in small segments of the data, which are then combined or averaged to reduce the noise in the estimate of the spectrum. It's fine to apply it to the full signal, but it's less commonly-used that way. You can do all this windowing and averaging manually, but I'd really recommend using the scipy.signal.welch function (or similar) for estimating the PSD. It takes care of all the book-keeping, windowing, averaging, etc for you.

python-FFT and PSD calculating

1 Answers