% modified from Aster's script to create IIR filter figures for his notes
%  display one figure at a time with a caption then pause for input
%
%filter examples and figures for notes
clear all
close all

imp=zeros(512,1);
imp(50)=1;

[b,a]=butter(1,0.1);
[H,F]=freqz(b,a,2048,1);

figure(1)
subplot(2,1,1)
plot(log10(F),20*log10(abs(H)));
ylabel('| \Phi | (dB)')
ylim([-60,3])
title('1-pole Butterworth Low Pass Filter')
subplot(2,1,2)
plot(log10(F),unwrap(180/pi*angle(H)));
ylabel('Phase (deg)')
xlabel('log_{10}(f)')
figure(11)
plot(filter(b,a,imp))
ylabel('Amplitude')
xlabel('Time')
title('1-pole Butterworth Low Pass Filter')

disp(' ');
disp('Frequency response of a 1-pole Butterworth Low Pass Filter');
disp(' using the matlab butter(1,0.1) and freqz(b,a,2048,1) commands.');
disp(' ');
disp('The Butterworth filter is designed to be optimally flat in the passband');
disp(' that is, it has no ripple.');
disp(' ');
disp('Note that frequency is normalized and log10(0.5)=-0.301,');
disp(' and the corner at normalized frequency 0.1 in log space is -1.');
disp(' ');
pause;

figure(2)
subplot(2,1,1)
plot(log10(F),20*log10(abs(H.*H)));
ylabel('| \Phi | (dB)')
ylim([-60,3])
title('1-pole Butterworth Low Pass Filter, Time-forward/reversed')
subplot(2,1,2)
plot(log10(F),unwrap(180/pi*angle(H.*conj(H))));
ylabel('Phase (deg)')
xlabel('log_{10}(f)')
figure(22)
plot(filtfilt(b,a,imp))
ylabel('Amplitude')
xlabel('Time')
title('1-pole Butterworth Low Pass Filter, Time-forward/reversed')

disp('Zero phase 1-pole Butterworth Low Pass filter using matlab filtfilt command.');
disp('The phase distortion can be removed by applying the filter in reverse.');
disp('This effectively doubles the order of the filter but also causes the');
disp('time domain response to be acausal');
disp(' ');
pause;

close all

[b,a]=butter(8,0.1);
[H,F]=freqz(b,a,2048,1);

figure(1)
subplot(2,1,1)
plot(log10(F),20*log10(abs(H)));
ylabel('| \Phi | (dB)')
ylim([-60,3])
title('8-pole Butteworth Low Pass Filter')
subplot(2,1,2)
plot(log10(F),unwrap(180/pi*angle(H)));
ylabel('Phase (deg)')
xlabel('log_{10}(f)')
figure(2)
plot(filter(b,a,imp))
ylabel('Amplitude')
xlabel('Time')
title('8-pole Butterworth Low Pass Filter')

disp('8-pole Butterworth Low Pass Filter using matlab butter(8,0.1) and freqz(b,a,2048,1) commands.');
disp('The greater number of poles provides a steeper rolloff (greater slope)');
disp(' in the transition band at the expense of increased distortion in the');
disp('  phase response.');
disp(' ');
pause;

[b,a]=cheby1(1,3,0.1);
[H,F]=freqz(b,a,2048);


figure(1)
subplot(2,1,1)
plot(log10(F/(2*pi)),20*log10(abs(H)));
ylabel('| \Phi | (dB)')
ylim([-60,3])
title('1-pole Chebyshev Type 1 Low Pass Filter (R=3 dB)')
subplot(2,1,2)
plot(log10(F),unwrap(180/pi*angle(H)));
ylabel('Phase (deg)')
xlabel('log_{10}(f)')
figure(2)
plot(filter(b,a,imp))
ylabel('Amplitude')
xlabel('Time')
title('1-pole Chebyshev Type 1 Low Pass Filter (R=3 dB)')

disp('1-pole Chebyshev Type 1 Low Pass Filter with no more than 3db of ripple');
disp(' using the matlab cheby1(1,3,0.1) and freqz(b,a,2048) commands.');
disp(' ');
disp('The Chebyshev filter is maximally flat in the stop band and minimizes');
disp('the difference between the actual and ideal response in the passband');
disp('while minimizing the ripple in the passband.');
disp(' ');
pause;

[b,a]=cheby1(8,3,0.1);
[H,F]=freqz(b,a,2048);

figure(1)
subplot(2,1,1)
plot(log10(F/(2*pi)),20*log10(abs(H)));
ylabel('| \Phi | (dB)')
ylim([-60,3])
title('8-pole Chebyshev Type 1 Low Pass Filter (R=3 dB)')
subplot(2,1,2)
plot(log10(F),unwrap(180/pi*angle(H)));
ylabel('Phase (deg)')
xlabel('log_{10}(f)')
figure(2)
plot(filter(b,a,imp))
ylabel('Amplitude')
xlabel('Time')
title('8-pole Chebyshev Type 1 Low Pass Filter (R=3 dB)')

disp('8-pole Chebyshev Type 1 Low Pass Filter with no more than 3db of ripple');
disp(' using the matlab cheby1(8,3,0.1) and freqz(b,a,2048) commands.');
disp('Again note the increased distortion in the phase response');
disp(' ');
pause;

[b,a]=cheby1(8,1,0.1);
[H,F]=freqz(b,a,2048);

figure(2)
subplot(2,1,1)
plot(log10(F/(2*pi)),20*log10(abs(H)));
ylabel('| \Phi | (dB)')
ylim([-60,3])
title('8-pole Chebyshev Type 1 Low Pass Filter (R=1 dB)')
subplot(2,1,2)
plot(log10(F),unwrap(180/pi*angle(H)));
ylabel('Phase (deg)')
xlabel('log_{10}(f)')
figure(3)
plot(filter(b,a,imp))
ylabel('Amplitude')
xlabel('Time')
title('8-pole Chebyshev Type 1 Low Pass Filter (R=1 dB)')

disp('8-pole Chebyshev Type 1 Low Pass Filter with no more than 1db of ripple');
disp(' using the matlab cheby1(8,1,0.1) and freqz(b,a,2048) commands.');
disp('Again note the increased distortion in the phase response due to the');
disp('  more strict stop band requirement.')
