Inverse FT of a filter in Matlab

Question

I am trying to find the inverse Fourier transform of a simple filter in Matlab. In the first case (sinc filter / "brick wall"), I use the ifft function to find the time-domain function, which is a sinc, centered at t = 0.

I now want to now find the time-domain function for a simple Chebyshev filter. However, for some reason, the code below seems to give me the impulse response, where the time axis is incorrect. Should I not expect a similar looking sinc function centered at t = 0?

fo = 10; %frequency of the sine wave
Fs = 100; %sampling rate
Ts = 1/Fs; %sampling time interval
t = -1+Ts:Ts:1-Ts; %sampling period
freq = -Fs/2:(Fs/length(t)):Fs/2;

%% Sinc with bandwidth = fo. This works!
y = 0.5*sinc(2*fo*t);
YfreqDomain1 = fft(y);
figure('Name','Brick wall sinc filter (freq)');
plot(freq(1:length(y)),2/length(y)*fftshift(abs(YfreqDomain1)))
y_ret1=ifft(YfreqDomain1,'nonsymmetric');
figure('Name','Brick wall sinc filter (time)');
plot(t,y_ret1);

%% Chebyshev with bandwidth fo. This gives me a strange result. 
[b,a] = cheby1(6,0.1,2*fo/Fs);  % 6th order, 0.1dB ripple
[YfreqDomain2 w] = freqz(b,a,length(t),'whole');
figure('Name','Chebyshev Filter (freq)');
plot(freq(1:length(YfreqDomain2)), 2/length(y)*fftshift(abs(YfreqDomain2)));
figure('Name','Chebyshev Filter (time)');
y_ret2=ifft(YfreqDomain2,'nonsymmetric');
plot(t,y_ret2);

I think the impulse response is the desired result. The filtered signal `y(t)` is the convolution between the input signal `x(t)` and the impulse response `h(t)`, so `h(t) = IDFT[H(w)]` is the impulse response. Or am I wrong there? — hbaderts, Dec 09 '14 at 14:56
The problem seems to be related to differences in the phase of the two signals. `plot(t,fftshift(ifft(abs(YfreqDomain2)))` looks correct and is comparable to the same command with YfreqDomain1. — Jim Quirk, Dec 09 '14 at 15:41
Not the cause of your problem, but also worth noting that YfreqDomain1 is 199x1 whereas YfreqDomain2 is 1x199 — Jim Quirk, Dec 09 '14 at 15:43
Thanks Jim - that worked. However, why is the plot not a symmetrical sinc-like function, as in the first case? It looks like an impulse response, and strangely, the period is twice the first case. hbaderts, that would be the case if x(t) was the Dirac function. However, in this case, we are simply doing a transform? — Graham, Dec 09 '14 at 15:46
To make sure we are seeing the same thing.... the plot command I gave above (assuming you add in the missing `)` at the end) should produce a symmetrical sinc-like function. As for why your original code shows a time-shifted version, that's currently a mystery to me — Jim Quirk, Dec 09 '14 at 16:11
Aha - my mistake. Thanks for the help! Not 100% sure why an abs() is needed though - surely this would eliminate phase information? — Graham, Dec 09 '14 at 16:21
Agreed. While the abs() creates the graph, it doesn't explain why your code didn't work. Hopefully someone else will know. — Jim Quirk, Dec 09 '14 at 16:32

A. Donda · Accepted Answer · 2014-12-09T20:42:34.847

There are two problems with your computation:

First, you evaluate the time domain filter coefficients on a very narrow time window, which truncates the filter and changes its characteristics.

Second, you do not properly keep track of which indices in the time domain and frequency domain vectors correspond to time point 0 and frequency 0, respectively. This can be done differently, I here chose to always have t(1) = 0 and f(1) = 0, and apply fftshift for plotting only.

Here's my corrected version of your code:

fo = 10;        % frequency of the sine wave
Fs = 100;       % sampling rate
Ts = 1 / Fs;    % sampling time interval
n = 1000;       % number of samples

% prepare sampling time vector such that t(1) = 0
t = (0 : n - 1)';
t = t - n * (t >= 0.5 * n);
t = t / Fs;

% prepare frequency vector such that f(1) = 0;
f = (0 : n - 1)' / n;
f = f - (f >= 0.5);
f = f * Fs;


%% sinc filter with bandwidth fo

% define filter in time domain
s = 2*fo/Fs * sinc(2*fo*t);

% transform into frequency domain
Hs = fft(s);
% since the filter is symmetric in time, it is purely real in the frequency
% domain. remove numeric deviations from that:
Hs = real(Hs);

subplot(2, 4, 1)
plot(fftshift(t), fftshift(s))
ylim([-0.1 0.25])
title('sinc, time domain')

subplot(2, 4, 2)
plot(fftshift(f), real(fftshift(Hs)), ...
    fftshift(f), imag(fftshift(Hs)))
ylim([-1.1 1.1])
title({'sinc, frequency domain', 'real / imaginary'})

subplot(2, 4, 3)
plot(fftshift(f), abs(fftshift(Hs)))
ylim([-0.1 1.1])
title({'sinc, frequency domain', 'modulus'})


%% Chebyshev filter with bandwidth fo

% define filter in frequency domain
[b,a] = cheby1(6, 0.1, 2*fo/Fs);  % 6th order, 0.1 dB ripple
Hc = freqz(b, a, n, 'whole', Fs);

% transform into time domain
c = ifft(Hc);

% determine phase such that phase(1) is as close to 0 as possible
phase = fftshift(unwrap(angle(fftshift(Hc))));
phase = phase - 2*pi * round(phase(1) /2/pi);

subplot(2, 4, 5)
plot(fftshift(t), fftshift(c))
title('Chebyshev, time domain')
ylim([-0.1 0.25])

subplot(2, 4, 6)
plot(fftshift(f), real(fftshift(Hc)), ...
    fftshift(f), imag(fftshift(Hc)))
ylim([-1.1 1.1])
title({'Chebyshev, frequency domain', 'real / imaginary'})

subplot(2, 4, 7)
plot(fftshift(f), abs(fftshift(Hc)))
ylim([-0.1 1.1])
title({'Chebyshev, frequency domain', 'modulus'})

subplot(2, 4, 8)
plot(fftshift(f), fftshift(phase))
title({'Chebyshev, frequency domain', 'phase'})

And here's the result:

As you can see, the sinc and Chebyshev filters are similar with respect to the modulus of the frequency response, but very different regarding the phase. The reason is that the Chebyshev filter is a causal filter, meaning that the coefficients in the time domain are constrained to be 0 for t < 0, a natural property for a filter that is implemented in a real physical system.

Inverse FT of a filter in Matlab

1 Answers1