FFT Fast Fourier Result not stable

I use FFT to analyze frequency and amplitude of a sequence of signal. But found out that to use the result of FFT to reconstruct the original signal is not always successful, even for two very closed signals. I am not sure if it is the limit of fft(x), or maybe I am not fully fully understand the function. The issue is shown by the following example:

Use signal 1: Signal=0.7*sin(2*pi*50*t), the reconstruction is successful:

Use signal 2: Signal=0.7*sin(2*pi*51*t), the only difference between two signal is frequency change from 50 to 51. The reconstruction is unsuccessful.

I try a few examples and find out:

1. Even frequency (ex. 50, 52, 58) usually have better reconstruction result from fft. Odd frequency (ex. 51, 53, 59) will have bad reconstruction.
2. I need a much longer signal sequence for odd freqency than even frequency to reconstruct the signal from fft. So it gives me trouble when I need to analyze very high frequency signal.

Here is the code I use for the above result.

Fs=1000; %Sampling frequency
T = 1/Fs;             % Sampling period       
L = 0.7;             % Length of signal by second
N = Fs*L;           % Total Sample number
t = (0:L*Fs-1)*T;        % Time vector
Signal=0.7*sin(2*pi*50*t); %Signal 1
%Signal=0.7*sin(2*pi*51*t); %Signal 2

%Compute the Fourier transform of the signal.
x = Signal;
y = fft(x);   
%Compute the two-sided spectrum P2. Then compute the single-sided spectrum P1 based on P2 and the even-valued signal length L.
P2 = abs(y/N); 
P1 = P2(1:N/2+1);
P1(2:end-1) = 2*P1(2:end-1);
%Define the frequency domain f and plot the single-sided amplitude spectrum P1. 
f=(0:(N/2))*Fs/N;
%Recreate Original Signal from result of Fourier transform
for time=1:N
    S(time)=0;
    for k=1:N/2+1
    S(time)=S(time)+P1(k)*sin(2*pi*f(k)*time*T);
    end
end
%plot result
figure
subplot(3,1,1);
plot(t,x);
title('Original');
subplot(3,1,2);
plot(f,P1);
title('Magnitude2')
subplot(3,1,3);
plot(t,x);hold on
plot(t,S);
title('Compare Original and Recreate Signal');

Please give me some help. Is there any way to improve the fft performance? Thank you.

[Walter Roberson]

N = Fs*L;           % Total Sample number

That will not generally be exactly an integer.

P1 = P2(1:N/2+1);

With N not generally being exactly an integer, then even if N is approximately even, N/2 might be just slightly less than what would be expected, so N/2+1 might be just slightly less than the expected integer, so 1:N/2+1 might be 1 shorted than you expect.

With N not generally being exactly an integer, then if N is approximately odd, N/2+1 would be approximately an integer plus 1/2, and 1:N/2+1 would only go as far as the integer without the 1/2 . But is that what is required, or should do you need it to go to the next location?

[HZ]

Thank you. N/2+1 is from Matlab's example, and it is for plot only.

Answer 1

[David Goodmanson]

Hi HZ,

It's good that you are working this stuff out. Not enought people do that.

WIth a few exceptions the result of an fft is a complex quantity. Once you do something like

P2 = abs(y/N); 
P1 = P2(1:N/2+1);

which happens all the time on this site, you have tossed away half the information. That's great for plotting but not great for anything else.

Once you are off into plot land, there is no getting back to fourier land. You only have abs(y), and you can't do the inverse transform with that. So it's better to ship the plotting results out to the plot command and stay with y.

Matlab ifft will get back to the signal, but I take your point that you want to do this explicitly with for loops. In the for loop calculation you have the correct time array, proportional to 0:N-1 as it should be and not 1:N. But the fourier coefficients aren't there.

Complex fourier coefficients mulitply exp(i* ...) and not sin(...) and cos)...). In your code if you replace the equivalent section with

%  Recreate Original Signal from result of Fourier transform
fall = (0:N-1)*Fs/N;         % all frequencies
S = zeros(1,N);
for timeind=1:N
   for k=1:N
      S(timeind)= S(timeind)+(y(k)/N)*exp(i*2*pi*fall(k)*(timeind-1)*T);
    end
end
S = real(S);                 % eliminate a tiny imaginary part due to numerical calculation

then the agreement is good for f = 51 or 59 or anything else. From personaI preference I changed the name of time to timeind since it is an index, and had to change timeind to (timeind-1) since the t array starts at t = 0.

[HZ]

Thank you. I have a a better understanding of FFT now.