Why am I getting wrong spectrum applying FFT comparing with correct spectrum?

Hello! I am trying to calculate the spectrum of an exponentially dampened oscillating field but I got the wrong result comparing with correct spectrum (different peak amplitude and FWHM). However, I tried with sine function and got a correct result. So I am not sure which part goes wrong (fft part or definition of correct spectrum with MATLAB). Thank you very much!

clear all
clc
Fs = 100;                        % Sampling frequency                    
T = 1/Fs;                        % Sampling period       
duration= 50;                    %second
L=duration*Fs;                   %sampling length 
t = 0:1/Fs:duration-1/Fs;        % Time vector
S = 0.7*exp(-t/5).*exp(-i*2*pi*5*t); % Original signal
figure (1)
plot(t,S)                       % plot original signal in time domain
Y = fft(real(S));               % fast fourier transform
P2 = abs(Y/L);                  % return it back to correct amplitude
P1 = P2(1:L/2+1);               % one side spectrum
P1(2:end-1) = 2*P1(2:end-1);    
f = Fs*(0:(L/2))/L;             % Frequency vector
figure (2)
plot(f,P1)                      % plot spectrum in frequency domain
w = 2*pi*Fs*(0:(L/2))/L;        % angular Frequency vector
w1=(0:0.001:100); 
FFT=(0.7/(2*pi))*(1./((1/5)-i*(w1-2*pi*5))); % define correct spectrum
figure (4)
plot(w1,real(FFT))             % plot correct spectrum in angular frequency domain
figure (5)
plot(w,P1)                     % plot spectrum in angular frequency domain

Answer 1

[David Goodmanson]

Hi Jiang,

There are a couple of things going on here. The expression for FFT is not correct. Since you are taking the real part of S, the transform is on the cosine, (1/2)(exp(+...)+exp(-...)) . The answer, Integral (......) dt, is

0.7*(1/2)* (1./((1/t0)-i*(w1-w0)) + 1./((1/t0)-i*(w1+w0)));  %     with  t0 = 5, w0 = 2*pi*5

There is no factor of 2pi in the denominator of this expression.

The fft approximates this integral. Taking the fft does the sum, and that result is multiplied by T, the width of each interval. So Y gets a factor of T instead of (1/L). Dividing by L is used in lots of situations, but not this one.

With the changes, the plots agree.

Fs = 100;                        % Sampling frequency                    
T = 1/Fs;                        % Sampling period       
duration= 50;                    %second
L=duration*Fs;                   %sampling length 
t = 0:1/Fs:duration-1/Fs;        % Time vector
S = 0.7*exp(-t/5).*exp(-i*2*pi*5*t); % Original signal
figure (1)
plot(t,S)                       % plot original signal in time domain
Y = fft(real(S));               % fast fourier transform
% P2 = abs(Y/L);                  % return it back to correct amplitude
P2 = abs(Y)*T;                  
P1 = P2(1:L/2+1);               % one side spectrum
P1(2:end-1) = 2*P1(2:end-1);    
f = Fs*(0:(L/2))/L;             % Frequency vector
figure (2)
plot(f,P1)                      % plot spectrum in frequency domain
w = 2*pi*Fs*(0:(L/2))/L;         % angular Frequency vector

% one sided spectrum
w1=(0:0.001:100); 
w0 = 2*pi*5;
t0 = 5;
FFT = 0.7*(1/2)*( 1./((1/t0)-i*(w1-w0)) + 1./((1/t0)-i*(w1+w0)) );
FFT = 2*FFT;                    % one sided


% plot abs(FFT) instead of real(FFT) since it is a fairer comparison to abs(Y).       
figure(4)
plot(w1,abs(FFT))             % plot correct spectrum in angular frequency domain
figure (5)
plot(w,P1)                     % plot spectrum in angular frequency domain

[Jiang Muqing]

Dear Goodmanson,

Very sorry for my late reply, thank you very much for you help! It completley solved my concern!