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

조회 수: 6(최근 30일)
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

채택된 답변

David Goodmanson
David Goodmanson 2022년 11월 1일
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
  댓글 수: 1
Jiang Muqing
Jiang Muqing 2022년 11월 13일
Dear Goodmanson,
Very sorry for my late reply, thank you very much for you help! It completley solved my concern!

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

추가 답변(0개)

Community Treasure Hunt

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

Start Hunting!

Translated by