How to obtain accurate fast fourier transformation

23 views (last 30 days)
Hi
I tried to extract the amplitude and frequency of a wave equation
but the values are not accurate.
For a given frequency w1=5 and amplitude A1=1,
I've obtained w1=49.805 and A1=0.933.
And for a given frequency w2=1 and amplitude A2=5
I've obtained w2=9.765 and A2=4.623
It looks like the frequency is about 10 times
higher (but still not exact) and the amplitudes is about
10% lower than the expected values.
To see this, run the script below.
So I hope someone has an idea of what I'm doing wrong.
Thanks
Emerson
clear all;
close all;
%Wave properties
A1=1; % Amplitude1 (dimensionless)
w1=5; % Frequency1 (Hz)
A2=5; % Amplitude2 (dimensionless)
w2=1; % Frequency2 (Hz)
% Time properties
Datapoints = 1000; % Number of recorded data;
Length=10; % Length (seconds);
Step= Length/Datapoints; % Fraction of seconds
t = (0:Datapoints-1)*Step; % Time vector
% Sum of a w1-Hz sinusoid and a w2-Hz sinusoid
y = A1*sin(2*pi*w1*t) + A2*sin(2*pi*w2*t);
% Fourier Transformation
NFFT = 2^nextpow2(Datapoints); % Next power of 2 from length of y
Y = fft(y,NFFT)/Datapoints;
f = Datapoints/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('Amplitude - |Y(f)|')

Accepted Answer

Ivan van der Kroon
Ivan van der Kroon on 17 Jun 2011
For your frequencies, you make a mistake in your frequency array
fs=Datapoints/Length;
f = fs/2*linspace(0,1,NFFT/2+1);
As you have to multiply it by half of the sampling frequency and not the number of data points. As you 1000 points per 10 seconds your sampling frequency is 100Hz.
The fact that your amplitudes are not exact is because the frequency components of your signal are not exactly in your discrete frequency array. An fft matches your signal by finding components close to your signal, i.e. it decomposes the signal in cosines and sines that it is allowed to sum over.
If you want to have them exactly in your transformed signal you have to change the sampling points of the fft. You now give it NFFT=1024, because fft's are faster for powers of 2, but this means that your frequency stepsize will be fs/1024= 0.0977. As you can see 1 nor 5 is a multiple of this number. You can however do it with fs/1000=0.1, so put NFFT=Datapoints; and you will have exact peaks at your predefined frequencies.
  3 Comments
Ivan van der Kroon
Ivan van der Kroon on 20 Jun 2011
This is true. For me it is easiest to remember that the frequency spepsize is df=fs/NFFT=1/L (if N=NFFT).
The range of your frequency array depends on the stepsize in your temporal (or lateral, etc) array. Call it dt, then the frequency (in Hz) is between -0.5/dt and 0.5/dt.
For example, to observe 5Hz you need at least dt=0.1. If you keep N=1000, you need L<100s. Next you want your df to be smaller than 5 and preferably a multiple of 5. So you can pick L=0.2 or multiples. This value is actually a lower bound: you can pick 0.3 but then your transform does not look so nice.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by