% example of digital filtering with lowpass filter

w1 = 400*2*pi;
fs = 10000;
w2 = 4000*2*pi;
t = 0:1/fs:2;
x = 0.5*cos(w1*t)+0.5*cos(w2*t);
sound(x,fs)

% design of 2-pole Butterworth filter with wc = 1000 hz
[z,p,k] = buttap(2);  % creates a 2-pole Butterworth filter
[num,den] = zp2tf(z,p,k);  % poles & zeros to transfer function 
% need to transform the cutoff frequency to 1000*2*pi
wc = 1000*2*pi;
[num,den] = lp2lp(num,den,wc);
T = 1/fs;
disp('Digital Filter')
[numd,dend] = bilinear(num,den,1/T)

% filter the signal and display the results
y = filter(numd,dend,x);  % filters x
figure(1)
subplot(211), plot(t(1:200),x(1:200))
ylabel('x(t)')
title('Lowpass filter')
xlabel('Time (sec)')
subplot(212),plot(t(1:200),y(1:200))
xlabel('Time (sec)')
ylabel('y(t)')
Omega = -pi:pi/100:pi;
H = freqz(numd,dend,Omega);
figure(2)
subplot(211),plot(Omega,abs(H))
title('Frequency Response')
ylabel('|H|')
xlabel('\Omega')
hold on
H12 = freqz(numd,dend,[w1/fs w2/fs]);
plot([w1/fs w2/fs],abs(H12),'ro')
hold off

% create wave files 
wavwrite(x,fs,'d:\web_pages\book\applets\chap11\x')
wavwrite(y,fs,'d:\web_pages\book\applets\chap11\ylp')
pause
sound(y,fs)





