232 views (last 30 days)

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 ?

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);

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).

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));

Opportunities for recent engineering grads.

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

Start Hunting!
## 0 Comments

Sign in to comment.