%  This script solves the example in lecture notes comparing filter orders
%  for the Chebyshev IIR filters. High-pass filters of orders 5 and 10 are
%  designed with a cut-off of 500 Hz and a stopband ripple of 30 dB down.
%  The freqz command is used to examine the frequency response
%  (phase and magnitude), the grpdelay command is used to examine the group
%  delay for the filter, and finally a frequency swept signal (from 20Hz
%  to 2000Hz is generated and filtered with both filters for comparisons.
% 
%    Written by Kevin D. Donohue (donohue@engr.uky.edu) March 2004 
%


%% DSP parameters
fs = 8e3;  %  Sampling Frequency
fc = 500;  %  Cutoff Frequency
rs = 30; %  stop band ripple in dB
dt = 2;  % Time duration of frequency sweep signal in seconds
flow = 20;   %  Starting frequency of sweep in Hertz
fend = 2000; %  Ending frequency of sweep in Hertz


%%  Create filter coefficients
[b5,a5] = cheby2(5,rs,fc/(fs/2),'high');  %  Order 5 chebyshev type II, low pass
[b10,a10] = cheby2(10,rs,fc/(fs/2),'high');  %  Order 10 chebyshev type II, low pass

%%  Compute Frequency Response
[h10,f10] = freqz(b5,a5,1024,fs);
[h50,f50] = freqz(b10,a10,1024,fs);
%%  Plot magnitude
figure(1)
plot(f10,20*log10(abs(h10)), 'b',f50,20*log10(abs(h50)), 'r')
title('Compare 5 and 10 order Chebyshev type II low-pass filter magnitudes')
xlabel('Hertz')
ylabel('dB')
legend('0rder 5', 'Order 10')
%%  Plot Phase
figure(2)
%  unwrap phase to eliminate mod pi jumps
plot(f10,unwrap(angle(h10)), 'b',f50, unwrap(angle(h50)), 'r')
title('Compare 5 and 10 order Chebyshev type II low-pass filter phases')
xlabel('Hertz')
ylabel('radians')
legend('0rder 5', 'Order 10')

%%  Compute Group Delay
[gd10,fgd10] = grpdelay(b5,a5,1024,fs);
[gd50,fgd50] = grpdelay(b10,a10,1024,fs);
figure(3)
plot(fgd10,gd10, 'b',fgd50,gd50 , 'r')
title('Compare 5 and 10 order Chebyshev type II low-pass filter group delays')
xlabel('Hertz')
ylabel('samples')
legend('0rder 5', 'Order 10')

%%   Simulate signal - Frequency Sweep
t = [0:round(dt*fs)-1]/fs;  %  Create time axis for sweep signal
fsw = flow + ((fend-flow)/2)*[0:length(t)-1]/length(t);  %  Create Frequency ramp - Why divid by 2?
swp = sin(2*pi*t.*fsw);  % Generate sweep signal

%%  Filter signal
sigsw10 = filter(b5,a5,swp);
sigsw50 = filter(b10,a10,swp);
figure(4)
plot(t,sigsw10)
title('Filtered Frequency Sweep - Order 5')
xlabel('Seconds')
ylabel('Amplitude')
figure(5)
plot(t,sigsw50)
title('Filtered Frequency Sweep - Order 10')
xlabel('Seconds')
ylabel('Amplitude')