CODEWITHSUNDEEP
matlabfilterfftifft
Graham
Graham

Reputation: 581

Inverse FT of a filter in Matlab

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

Upvotes: 0

Views: 1470

Answers (1)

A. Donda
A. Donda

Reputation: 8467

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.

Upvotes: 2

Related Questions