MATLAB Answers

FFT - am i doing anything wrong

232 views (last 30 days)
Adam
Adam on 20 Feb 2012
Edited: Cedric Wannaz on 19 Oct 2013
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?
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 ?

  0 Comments

Sign in to comment.

Accepted Answer

Wayne King
Wayne King on 20 Feb 2012
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);

  3 Comments

Adam
Adam on 21 Feb 2012
thank you, both of you... really helped :)
Adam
Adam on 21 Feb 2012
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
Dr. Seis on 21 Feb 2012
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).

Sign in to comment.

More Answers (1)

Dr. Seis
Dr. Seis on 20 Feb 2012
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));

  0 Comments

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by