fft is not matching with analytical solution
조회 수: 8 (최근 30일)
이전 댓글 표시
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?
댓글 수: 0
채택된 답변
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
참고 항목
카테고리
Help Center 및 File Exchange에서 Spectral Measurements에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!