FFT - am i doing anything wrong

hi guys

I have a large vector with frequency data. That is measurements of the frequency on the power grid at a given time. Each sample has been sampled with an interval of 5 seconds. I want to make a fft analysis of these data, so i can see the components of this new frequency data, but i really can't figure out how. It seems simple if you follow the examples in the documentaion fft doc, but when i try with my data vector it doesn't get correct. Or maybe it is correct, but it seems strange i think? Am i doing anything wrong?

This is a plot of my signal http://imageshack.us/photo/my-images/16/freqsignal.jpg/

and this is the frequency spectrum http://imageshack.us/photo/my-images/708/freqspec.jpg/

and this is my code: data = importdata('data.txt'); freq = data(:,1)

figure(1)
Fs = 5;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = length(freq);                     % Length of signal
t = (0:L-1)*T;                % Time vector 
y = freq    % Sinusoids plus noise
plot(Fs*t(1:end),y(1:end))
title('freq vs time')
ylabel('Freq [mHz]')
xlabel('time (sec)')

figure(2)
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1))) 
title('frequency spectrum')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

I just took these lines from the example in fft doc. what does these lines mean?

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
f = Fs/2*linspace(0,1,NFFT/2+1);

and why does he use plot(f,2*abs(Y(1:NFFT/2+1))) instead of just plot(f,abs(Y)) ??

I hope somebody can explain me this, cause i am lost ?

Answer 1

[Wayne King]

Hi Adam, the main problem is that your data does not have zero mean.

This means that the zero frequency component (which isn't very interesting) is going to obscur what is probably more interesting.

first do this:

y = detrend(y,0);

Then take the Fourier transform of that zero-mean y.

To answer your other questions:

When you are dealing with real-valued data, the magnitudes of the Fourier transform are an even function of frequency so you only have to plot 1/2 the frequencies.

plot(f,2*abs(Y(1:NFFT/2+1)))

And the following two lines are just setting the number of points for the fft() and creating a meaningful frequency vector so you can plot the magnitudes as a function of frequency in Hz.

    NFFT = 2^nextpow2(L); % Next power of 2 from length of y
    f = Fs/2*linspace(0,1,NFFT/2+1);

[Adam]

thank you, both of you... really helped :)

[Adam]

One more thing please... what should the sample frequency Fs really be? cause when i change this value, the frequency scale also changes!

If i sample the real signal with an interval of T = 2 seconds, does that mean that Fs = 1/2?

[Dr. Seis]

Yes... T = 1/Fs.

Incidentally, why are you multiplying your time samples by Fs? When you do that your x-axis is just sample number... you are essentially doing (0:L-1)*T*Fs which is just (0:L-1).

Answer 2

[Dr. Seis]

Although the scaling performed on the result from FFT (in the example above, which is based on Matlab's FFT documentation) may be useful for some reason, it is not really appropriate for the purposes of performing meaningful scientific analysis (which depend on correct amplitudes).

I have shown in other posts that the result of the FFT should be scaled by "T" (equal to 1/Fs), but it doesn't hurt to reiterate. In the example from the documentation, the result of the FFT is normalized by "L" (equal to 1000), which is conveniently equal to "Fs". So the part where they divide the FFT result by "L" has the same effect as dividing by "Fs". However, as soon as "L" is changed to something other than "Fs" the scaling will not be correct. If the amplitudes are scaled appropriately, then the energy in the frequency domain should be equal to the energy in the time domain. I.e., if:

Fs = 1000;           % Sample Rate
N = 100;             % Number of samples
N2 = 2^nextpow2(N);  % Number of samples to next power of 2
T = 1/Fs;            % Time increment
F = Fs/N2;           % Frequency increment
y_time = randn(1,N); 
Y_FREQ = fft(y_time,N2)*T;

Then,

energy_time = sum(y_time.*conj(y_time))*T;
ENERGY_FREQ = sum(Y_FREQ.*conj(Y_FREQ))*F;

And,

"energy_time" should be equal to "ENERGY_FREQ" within machine precision.

I prefer to use "fftshift" to organize the amplitudes in the frequency domain before displaying:

Fs = 1000;           % Sample Rate
N = 100;             % Number of samples
N2 = 2^nextpow2(N);  % Number of samples to next power of 2
T = 1/Fs;            % Time increment
F = Fs/N2;           % Frequency increment
Nyq = Fs/2;          % Nyquist frequency
y_time = randn(1,N); 
Y_FREQ = fftshift(fft(y_time,N2))*T;

plot(-Nyq:F:Nyq-F , abs(Y_FREQ));