0
$\begingroup$

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???

$\endgroup$
2
$\begingroup$

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.

The Mathworks has a nice example comparing circular with linear convolution.

IIR filters are best implemented in the time domain. For their specs, they are usually much lower order than FIR filters and so we can afford to implement them directly.

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Niwol May 16 '17 at 12:50
  • $\begingroup$ 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). $\endgroup$ – Peter K. May 16 '17 at 13:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.