MATLAB Answers

fft is not matching with analytical solution

조회 수: 12(최근 30일)
subbu
subbu 2021년 3월 9일
댓글: subbu 2021년 3월 9일
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?

채택된 답변

Mathieu NOE
Mathieu NOE 2021년 3월 9일
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
  댓글 수: 1
subbu
subbu 2021년 3월 9일
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

댓글을 달려면 로그인하십시오.

추가 답변(1개)

Matt J
Matt J 2021년 3월 9일
편집: Matt J 2021년 3월 9일
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;
  댓글 수: 1
subbu
subbu 2021년 3월 9일
even if i make the changes that you mentioned, i am not getting the desired answer that i told earlier

댓글을 달려면 로그인하십시오.

태그

Community Treasure Hunt

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

Start Hunting!

Translated by