Changed system G to a (stable) closed loop system, results seem to match a lot better now! Only some transient behaviour on the start of the lsim data
Fourier transform vs lsim, results dont match.
1 view (last 30 days)
Show older comments
If a signal u(t), with Fourier transform U(e^jω) is applied to system G, then the Fourier transform of the output y(t) is given by:
Y(e^jω)=G(e^jω)U(e^jω).
I try to verify this in MATLAB, but something goes wrong.
Implementation:
-U(e^jω) is obtained easily enough using fft(u).
-For G(e^jω) I use freqresp(sys,ω) (and build ω in such a way that it is in the same format as the fft results: [dc positive-freqs negative-freqs]). I verify this by comparing the result to the bode plot (see commented lines in script)
To verify I perform the inverse Fourier transform y = ifft(Y) and compare the obtained y(t) to the lsim result, which is obtained with yprime(t)= lsim(G,u). Unfortunately these dont match. What am I missing? The lsim results appear to have some integrator effect in them???
clear
close all
G = tf([4700 4393 3.245e08],[1 7.574 1.202e5 0 0]);
f_sh = 1000; %sampling frequency
Gd = c2d(G,1/f_sh); %plant discretization
t_end = 10; %[s] simulation time end
t = linspace(0,t_end,t_end*f_sh+1); %time vector
f_ref = 0.5; %reference freq
u = sin(2*pi*f_ref*t); %reference signal
L = length(u);
dF = f_sh/L; %frequency bin step
w = dF:dF:f_sh/2; %one-sided frequency vector in Hz
w_eval = 2*pi*[0 w -fliplr(w)]; %frequency vector for freqresp, formatted to match fft results
H = freqresp(Gd,w_eval); %frequency response of system wrt freq vec
H = transpose(squeeze(H)); %remove singelton dimensions
% bode(Gd) %should be equal to the semilog plot
% figure
% semilogx(w,20*log10(abs(H(2:(L+1)/2)))) %plot in Hz
U = fft(u);
Y = H.*U;
y = ifft(Y); %frequency domain response
yprime = lsim(Gd,u); %simulation/time domain resonse
figure
plot(t,yprime,t,y) %y should equal yprime
legend('lsim response','fft(sys)*fft(ref) response')
Answers (0)
See Also
Categories
Find more on Digital Filter Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!