# MATLAB: Fourier transform vs $\tt lsim$, results don't match

My (simplified) problem: Suppose I have a signal $u$, which passes through system/filter $G$, which gives output signal $y$. Then I have a theorem that states:

If a signal $u(t)$, with Fourier transform $U(e^{j \omega})$ is applied to system $G$, then the Fourier transform of the output $y(t)$ is given by:

$$Y(e^{j \omega})=G(e^{j \omega})U(e^{j \omega})$$

I try to verify this in MATLAB, but something goes wrong.

Implementation:

• $U(e^{j \omega})$ is obtained easily enough using fft(u).

• For $G(e^{j \omega})$ I use freqresp(sys,$\omega$) (and build $\omega$ 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 $y_{prime}(t)$ = lsim(G,u). Unfortunately these dont match. What am I missing?

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')


The lsim results appear to have some integrator effect in them???

To perform the FFT approach correctly, you will need the FFT length to be at least $N_U + N_H - 1$ where $N_U$ is the length of $u$ and $N_H$ is the length of $h$.
It looks like $h$ is an IIR filter, so the FFT approach probably won't work well.
• Could you elaborate somewhat? Currently signal u has a certain length(which is the fft length), which also determines the length of H (cant multiply if not the same length) so I dont see how I can increase the fft length to $N_U+N_H−1$ And why wouldnt the FFT approach work for an IIR filter? – Niwol May 16 '17 at 12:50
• The FFT approach does circular convolution. If the FFT length is $N$, then the circularity is length $N$. Any convolutions performed with an FFT of length $N$ will be circular convolutions of this length. If $N \ge N_U + N_H -1$ then the circular convolution result is the same as the standard convolution result. If that inequality is not satisfied, then time-domain aliasing occurs (i.e. the beginning and ends of the result are corrupted if you were expecting a linear convolution result). – Peter K. May 16 '17 at 13:45