fft is not matching with analytical solution

clear;
clc;

% Generate the sine wave sequence
fs=8000; % Sampling rate
N=1000; % Number of data points
x1=sin(500*pi*[0:1:N-1]/fs)...
    +2*sin(500*1.5*pi*[0:1:N-1]/fs).*heaviside([0:1:N-1]/fs-300/fs)...
    +3*sin(500*3*pi*[0:1:N-1]/fs).*heaviside([0:1:N-1]/fs-600/fs)...
    +4*sin(500*4*pi*[0:1:N-1]/fs).*heaviside([0:1:N-1]/fs-800/fs);
t=[0:1:N-1]/fs;
x=x1;
% Apply the DFT algorithm
figure(1)
xf=abs(fft(x))/N; % Compute the amplitude spectrum
P=xf.*xf; % Compute power spectrum
f=[0:1:N-1]*fs/N; % Map the frequency bin to frequency (Hz)
subplot(2,1,1); plot(f,xf);grid
xlabel('Frequency (Hz)'); ylabel('Amplitude spectrum (DFT)');
subplot(2,1,2);plot(f,P);grid
xlabel('Frequency (Hz)'); ylabel('Power spectrum (DFT)');
figure(2)
% Convert it to one side spectrum
xf(2:N)=2*xf(2:N); % Get the single-side spectrum
P=xf.*xf; % Calculate the power spectrum
f=[0:1:N/2]*fs/N ; % Frequencies up to the folding frequency
% subplot(2,1,1);
plot(f,xf(1:N/2+1));
grid
xlabel('Frequency (Hz)'); ylabel('Amplitude spectrum (DFT)');
figure (3)
plot(t,x1);grid
xlabel('Time (s)');
ylabel('Amplitude');
title('Signal');

here i have created a signal and used fft to represent the signal in frequency domain. the analytical solution for the above signal is

F(ω)=f ̂(x(t))=1/2i* (δ(ω-250)+δ(ω+250))+2/2i* (δ(ω-375)+δ(ω+375))+3/2i* (δ(ω-750)+δ(ω+750))+4/2i* (δ(ω-1000)+δ(ω+1000))

which means that i can expect peaks at 250Hz,375Hz,750Hz and 1000Hz with amplitudes of 1,2,3 and 4 respectively in a single side amplitude spectra. i can see that the peaks occuring at this points but their amplitudes showing different values from the analytical solution. anyone can help me out by explaining why it is happening?

Answer 1

[Mathieu NOE]

Hello

I had to remove the heaviside as I dn't have the toolbox of it;

Anyway this code gives the right plot of fft amplitude for stationnary sine waves. Please not I changed Fs to 5000 so that with Nfft = 1000 I have a frequency vector spacing of 5 Hz , so all signal frequencies would fall exactly on one frequency bin of the fft. This ensure maximal amplitude accuracy (otherwise you can be affaected by leakage , the result will depend on which type of window you use)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Fs = 5000;
samples = 1000;
t = (0:samples-1)'*1/Fs;

signal=sin(500*pi*t)+2*sin(500*1.5*pi*t)+3*sin(500*3*pi*t)+4*sin(500*4*pi*t);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = samples;    
Overlap = 0.75;
w = hanning(NFFT);     % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display : averaged FFT spectrum 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq,fft_spectrum] = myfft_peak(signal, Fs, NFFT, Overlap);

figure(2),plot(freq,fft_spectrum,'b');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(' Amplitude')

function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples)) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select);
freq_vector = (select - 1)*Fs/nfft;
end

[subbu]

thank you for your effort but its look like u have rewritten the entire code. i will be happy if you can tell me whats the mistake that i have done in my code also if i make N=10^5(sample points) i am getting the output that i want. i thinik my code have some conceptual mistake that i am not sure

Answer 2

[Matt J]

Different professions define the scaling of the continuous Fourier Transform differently. If the definition that you use is,

then the appropriate scaling of the FFT is,

xf=abs(fft(x))/fs;

[subbu]

even if i make the changes that you mentioned, i am not getting the desired answer that i told earlier