
% IIR filter Design
% Butterworth   - no ripple
% Chebycheff I  - ripple in pass band
% Chebycheff II - ripple in stop band
% Elliptic      - ripple in pass and stop band; narrow transition bands


Continuous
time domain <--> freq domain Laplace transform
                        ^
                        |  impulse invariance
                        |  bilinear transform
Discrete                v
time domain <--> freq domain: Z-transform

1. Design analog lowpass filter
2. Frequency transformation (e.g. lowpass to highpass)
3. Filter transformation (e.g. bilinear)
4. Desired digital filter

Frequency transformation in continuous or discrete domain

1. Approximational Derivative (more complex than analog)
2. Impulse invariance  (must be bandlimited, i.e. lowpass & bandpass only)

3. Bilinear transformation:
     2   1 - z^-1
s = ---  -----
     T   1 + z^-1

Continuous domain -->
% Comparison of analog filters
[z,p,k] = buttap(5);
w = logspace(-1,1,1000);
h = freqs(k*poly(z),poly(p),w);
subplot(411), semilogx(w,abs(h)), grid;

[z,p,k] = cheb1ap(5,0.5);
w = logspace(-1,1,1000);
h = freqs(k*poly(z),poly(p),w);
subplot(412), semilogx(w,abs(h)), grid;

[z,p,k] = cheb2ap(5,20);
w = logspace(-1,1,1000);
h = freqs(k*poly(z),poly(p),w);
subplot(413), semilogx(w,abs(h)), grid;

[z,p,k] = ellipap(5,0.5,20);
w = logspace(-1,1,1000);
h = freqs(k*poly(z),poly(p),w);
subplot(414), semilogx(w,abs(h)), grid;

% Different type and design
[b,a] = butter(5,0.4);            % Lowpass Butterworth
figure, freqz(b, a);
[b,a] = cheby1(4,1,[0.4 0.7]);    % Bandpass Chebyshev type I
figure, freqz(b, a);
[b,a] = cheby2(6,60,0.8,'high');  % Highpass Chebyshev type II
figure, freqz(b, a);
[b,a] = ellip(3,1,60,[0.4 0.7],'stop');% Bandstop elliptic
figure, freqz(b, a);
figure, zplane(b,a);

% Butterworth with different taps
[b,a] = butter(2,0.4);            % Lowpass Butterworth
[H1, W] = freqz(b, a);

%figure, zplane(b,a);

[b,a] = butter(4,0.4);            % Lowpass Butterworth
[H2, W] = freqz(b, a);

[b,a] = butter(8,0.4);            % Lowpass Butterworth
[H3, W] = freqz(b, a);

[b,a] = butter(15,0.4);            % Lowpass Butterworth
[H4, W] = freqz(b, a);

[b,a] = butter(50,0.4);            % Lowpass Butterworth
[H5, W] = freqz(b, a);

H = [H1 H2 H3 H4 H5];
figure, freqzplot(H, W, s);

% Chebyshev II with different ripple level
[b,a] = cheby2(20, 20,0.4);    % Lowpass Chebyshev type I
figure, freqz(b, a);
[b,a] = cheby2(20, 60,0.4);    
figure, freqz(b, a);
[b,a] = cheby2(20, 100,0.4);   figure, freqz(b, a);
figure, freqz(b, a);

% Exmple of ellipord
[n,wn] = ellipord(0.4,0.5,1,100)
[b,a] = ellip(n,1, 100, wn);
figure, freqz(b,a);

[n,wn] = ellipord(0.4,0.45,1,100)
[b,a] = ellip(n,1, 100, wn);
figure, freqz(b,a);

[n,wn] = ellipord(0.4,0.41,1,100)
[b,a] = ellip(n,1, 100, wn);
figure, freqz(b,a);


% Example of manual transformation: analog,freq,bilinear

fdatool

x = wavread('max.wav');
y = filter(Num, Den, x);
wavwrite(y, 22050, 'maxF1.wav');
y = filter(Num, Den, x);
wavwrite(y, 22050, 'maxF2.wav');

