Applying Low-Pass filter to a signal with the Bilinear Method - MATLAB

Question

I have the following code:

Fs = 8000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
N=512;                        % Length of signal
t = (0:N-1)*T;                % Time vector


% Sum of a 1kHz sinusoid and a 2.5kHz sinusoid
x1=4*sin(2*pi*1000*t);
x2=sin(2*pi*2500*t); 
x=x1+x2;

% Frequency spectrum
y=fft(x);
fx = Fs/2*linspace(0,1,N/2+1);
plot(fx,2*abs(y(1:N/2+1))) 

i want to use the following low-pass filter on this signal:

[B A]=lp2lp([1],[1 1], 1.069*10^4)
[b a]=bilinear(B,A,Fs)
z=freqz(b,a,512,Fs);
axis ([0 Fs/2 -20 1])

any thought on how do it ?

Answer 1

Short Answer:

1. Using freqz: you only get the single-sided frequency response of the filter and hence you can only apply the filter on the single-sided spectrum of the signal, x to get the single-sided frequency spectrum of the output (filtered) signal o1:

   z = freqz( b, a, N/2+1, Fs);
   o1 = z' .* y(1:N/2+1);
   

2. Using filter: you can bypass using freqz if your sole goal is to apply the filter found earlier:

   o2 = filter( b, a, x); 
   

Optimal use of freqz:

freqz returns the the frequency response of the filter given the coefficients of its transfer function as described by mathworks. It is best used without returning a value, in which case, it automatically plots the magnitude and the phase response for the right half of the frequency spectrum upto half of the sampling frequency.

We can verify that the magnitude of the frequency response, plotted by the freqz command when called without a return value, is equivalent to the response plotted using the following code:

h = freqz(b,a,N/2+1,Fs);
plot( fx, 20*log10( abs(h) ));

Ouptut:

Moreover, we can see that the effect of the bilinear filters obtained using either of the two methods described above on the frequency response y of the original signal is precisely equivalent:

subplot(2,2,[1 2]); plot( fx, 2*abs( y(1:N/2+1))); ylim([0,2000]);
subplot(2,2,3); plot( fx, 2*abs( o2(1:N/2+1)));
subplot(2,2,4); plot( fx, 2*abs( o1));

P.S: Might be useful to go over FIR filtering overview in Matlab.