8 views (last 30 days)

Show older comments

Hi, I have a time series of temperature.I calculated fft for that. But when I calculate frequency, I am not sure whether it is correct or not. Can anyone please help me.Here is my code

t=10:10:2500;

Fs=0.1;

Y=fft(detrend(x,'constant'),256);

f=(0.1/256*(0:127));

Pyy=abs(Y.*conj(Y)/256);

plot(f,Pyy)

Eric
on 13 Oct 2014

Edited: Eric
on 13 Oct 2014

Here's what I would do:

N = 256; %Length of fft after padding

delta_t = 10;%seconds

delta_f = 1/(N*delta_t);%Hertz

f = (-fix(N/2):1:fix((N-1)/2))*delta_f; %Hertz; This handles the case where N is even or odd

Y = fftshift(fft(ifftshift(detrend(x,'constant')),N))/sqrt(N);

Pyy = abs(Y.^2);

plot(f,Pyy);

The interior ifftshift is irrelevant since you're calculating a power-spectral density at the end. It's utility is also diminished since you're padding the FFT with zeros at the beginning. Still, as I said, it's what I would do, though. If you're ever working with the phase of the Fourier transform of something it's usually advisable to shift prior to the FFT in addition to after the FFT.

All right, I guess I have to explain that. The Fourier transform of a real, even function is another even, real function. Consider the vector:

A = [0 0 1 1 1 0 0];

A is a real, even sequence about A(4). Calculate

B = fftshift(fft(ifftshift(A)));

and you get a real, even sequence about B(4) as expected. If you calculate

C = fftshift(fft(A));

you do not get a real function. This situation is subtler if you let A have an even length. Let

A = [0 0 1 1 0];

In that case

B = fftshift(fft(ifftshift(A)));

C = fftshift(fft(A));

both return a real function. However, the phase of C is still screwed up - just in a way that maintains a real function. The signs of every other data point are wrong. You can use the Fourier shift theorem to understand why. The DFT assumes that you are not working in a centered coordinate system. By working in a centered coordinate system you have shifted the data by N/2 data points. The Fourier shift theorem gives us:

FT{g(t-N/2)} = G(f)*exp[-j*2*pi*f*N/2]

where G(f)=FT{g(t)}. At the point where f=0, the phase term on the right-hand-side is 0. At the next data point over (i.e., when f=1/N), the phase term on the right is -pi. Thus by working in a centered coordinate system we introduce an alternating 0,-pi,0,-pi,... variation in the phase of the transform. By initially shifting the data we convert to a non-centered coordinate system to avoid this phase error.

Again, since in the end you're calculating a power-spectral density, this phase error is not going to impact the results. It's still good to know it's there, though.

Good luck,

Eric

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!