How to obtain accurate fast fourier transformation

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

Answer 1

[Ivan van der Kroon]

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.

[Emerson De Souza]

Thank you a lot Ivan.

I wrote the changes as you have suggested:

1)fs=Datapoints/Length;

f = fs/2*linspace(0,1,NFFT/2+1);

2)NFFT=Datapoints;

And the output is as expected.

However, I have another question:

I replaced the Length from 10 to 0.90 and 100 seconds

and the output is wrong again.

The range to work correctly is:

1<Length<99

I thought that the fourier transformation is sensitive even

to fractions of a second and that higher length should

not affect at all. I guess that I still need to improve the code

but I don't know how. Do you know what I'm doing wrong?

Thank you in advance for your support

Emerson

[Teja Muppirala]

If the length is 100 seconds and your Datapoints is 1000, then your sample time is 0.1 seconds.

sin(5*(2*pi)*0.0) = 0

sin(5*(2*pi)*0.1) = 0

sin(5*(2*pi)*0.2) = 0

sin(5*(2*pi)*0.3) = 0

The 5 Hz sine wave does not show up at all because it is not sampled fast enough. This is called aliasing.

[Ivan van der Kroon]

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.