IIR FILTER DESIGN BY TRANSFORMATION

Begin with rational continuous-time filters as basis of IIR design:




Passband: (1-d1) < |Hc(jW) | < 1 for |W| < Wc

Stopband: |Hc(jW) | < d2 for Wr < |W|



BUTTERWORTH FIlTER DESIGN:








NOTE:

N=1--> 6 dbs/octave rolloff

N=2-->12 dbs/octave rolloff

N=3-->18 dbs/octave rolloff

N=4-->24 dbs/octave rolloff

N=8-->48 dbs/octave rolloff

POLES WHEN:


(s)2N= -1 (jWc)2N= -1 (j)2N (Wc)2N

= -1 (-1)N (Wc)2N= (-1)N+1 (Wc)2N

Therefore: sk = Wc(-1)(N+1)/2N

or where (-1)(N+1)/2N = 2N roots of +1 if N is odd,
= 2N roots of -1 if n is even:

n odd: (1)1/2N = (ej2pk)1/2N = ej(pk/N) , k = 0, 1, ...2N-1
n even: (-1)1/2N = (ej(-p+(2pk)1/2N
= ej(-p/2N + pk/N) , k = 0, 1, ...2N-1



Now: |HB(s)|2 = HB(s) HB(-s)

Since we want a stable system, we only want those poles that are in the left half plane:

sk = Wc ej(p/2N)(2k+N-1) k = 1, 2, ...N

e.g., N=3, k=1: s1 = Wc ej(p/6)(4) = Wc ej2p/3


N = EVEN NUMBER OF POLES

M = N/2= NUMBER OF CASCADES





where: HB(s) = H1(s) H2(s) H3(s)....HN/2(s)

sk+ sk* = Wcejfk + Wce-jfk = 2 Wccos Fk

sksk*= Wc2


DESIGN EXAMPLES

METHOD 1: Beginning with any two frequency requirements. Particularly useful if you want your passband to be better than the usual definition of Wc requires, e.g. down 1 dB rather than 3 dBs.
Example:fp = .1 , fr = .15
20 log10(1 - d1) = -1 dB d1 = 0.10875
1 - d1 = .89125
20 log10d2 = - 15 dB d2 = 0.17783
Then:Wp = .2p Wr = .3p








Turns out: N = 5.8858, so we need to pick N of 6. We could decide to always pick an even N.

Two choices for Wc: Pick one to match exactly and the other will be better:

If we go to passband equation: Wc = .7032, stopband will be better than expected.

METHOD 2:

ASSUMING THAT Wc IS GIVEN: JUST NEED TO SOLVE FOR N





METHOD 3: For large enough W (relative to Wc)





So first determine how many octaves between Wc and Wr and the amount of attenuation wanted, e.g.,

Wc = 1 and Wr = 4, and you want a 24-db attenuation at 4.






CHEBYSHEV FIlTER DESIGN:

TYPE I:

RIPPLES IN THE PASSBAND

TYPE II:


RIPPLES IN THE STOPBAND

To(X) = 1

T1(X) = X

T2(X) = 2X T1(X) - To(X) = 2X2 - 1

T3(X) = 2X T2(X) - T1(X) = 4X3 - 2X - X

TN(X) = 2X TN-1(X) - TN-2(X) etc.

SPECIAL FEATURES

:

TN(X) = cos(N cos -1(X)) = cosh(N cosh -1 X)


TN(1) = 1

ALWAYS EVEN OR ODD FUNCTIONS

ALWAYS HAVE N REAL ZEROS BETWEEN -1 and +1

THE FIGURE BELOW SHOWS FROM WHENCE COME THE RIPPLES

The values are created very simply with MATLAB

x=[-2:4/128:2];N=3; fx=cos(N*acos(x));



EVENS ONLY


Also note that all polynomials are squared in the Chebyshev designs so 1+f2(x) > 0, so no "divide by zeros"

see claliir.m for comparison of Chebyshev and Butterworth

TRANSFORMATION # 1: Impulse Invariance:


NOTE HOW: jW ----> ejw = ejWT, s = s + jW, z = esT

This implies that z=e s + jW= e s e jW = r e jW

So if s <0, i.e., s is in the left half of the s-plane, e s =r<1.

Therefore, the entire left half of the s-plane maps into the unit circle on the z-plane.

So if:





However, this gives us the relationship described at first. If we want our relationship to look like:

H(ejw) = Hc(jw/T); |w|<= p, use h(n) = T hc(nT)

EXAMPLE:

hc(t) = 2e-t u(t) + e-2t u(t)
h(n) = T2e-nT u(n) + Te-2nTu(n)
= 2T(e-T)nu(n)+T(e-2T)nu(n)



parallel: y1(n) = e-Ty1(n-1)+ 2Tx(n)
y2(n) = e-2Ty2(n-1)+ Tx(n)
y(n) = y1(n)+ y2(n)

cascade:

y(n)=(e-2T+e-T)y(n-1)-e-3Ty(n-2)+3Tx(n)-(2Te-2T+Te-T)x(n-1)

GENERAL 2nd ORDERS:




NOTE THE RELATIONSHIP BETWEEN THE POLES:

p (s-plane) ---> epT(z-plane)

The MATLAB residue and residuez routines are nice routines that can help you with this computation.

S-z PLANE TRANSFORMATIONS, i.e., H(z) = Hc(s)|s=f(z)

Would be nice if we could just make some substitution directly in H(s) to get H(z). To get ideas, we look at some simple digital operations:

E.G. Differentiation:



Simplest difference equation which does differentiation is:




The most widely used substitution comes from another operation designed to avoid aliasing. Jackson approach is to first compress the s plane

)


What does this do to our old relationships between w and W ?


So we still have the unit circle in the z plane covered by the imaginary axis in the x plane:

However, there is one BIG difference -- it is only covered once!

Because as W goes to inf, w goes to p




z = es'T or s' = ln z /T or Ts' = ln z




and since e ln y = y



Solve for z in terms of s:






This process is known as frequency warping. The process will change the parameters of our filter design. i.e. wc and wr must be prewarped according to the relationship

and for Low pass filters, we can actually get by with a smaller N.
Redo:fp = .1 , fr = .15
20 log10(1 - d1) = -1 dB d1 = 0.10875
1 - d1 = .89125
20 log10d2 = - 15 dB d2 = 0.17783
Then:Wp = .2p

pass=2tan(0.1p)

Wr = .3p

stop=2tan(0.15 p)







Turns out: N = 5.30466, instead of 5.8 so we still need to pick N of 6, but sometimes we might be able to reduce N. We could decide to always pick an even N.

ZERO ADJUSTED IMPULSE INVARIANT

We begin with a design for which we know what we want the DC component to be: For the lowpass filter, this is usually 1. e.g. the single pole RC circuit:






Then: Ha(1) = 1 or you could design the system to match at any given point:



which is always a check: WHY? Hint which is pretty handy at times.



SAME PROBLEM: BILINEAR TRANSFORMATION:






AS THE ALGEBRA IS SOMETIMES MESSY, JACKSON RECOMMENDS DOING THIS WITH ZERO-POLE PLOTS




use:

;

;

SO FOR OUR LITTLE ONE POLE EXAMPLE:



WHICH IS THE SAME AS BY DIRECT SUBSTITUTION.

A COMPARISON OF THE TWO APPROACHES

T=1, RC=1



MATLAB HAS SOME WONDERFUL SUBROUTINES FOR DOING IIR FILTER DESIGN. SUGGEST YOU DO THE MATLAB EXAMPLES IN THE CHAPTER--PARTICULARLY EXERCIE NO. 15 ON CLASSICAL FILTER DESIGNS.

THE FIGURES BELOW SHOW A COMPARISON

1. AN IMPULSE INVARIANT TRANSFORMATION OF THE BUTTERWORTH DESIGN ([b,a]=impinvar(bc, ac) which takes b's and a's from analog design to digital design. bc's and ac's can be obtained using "buttap". The zeros are not plotted on the zp plot because one is a very large real number. The analog signal has no finite zeros.

2. A BILINEAR TRANSFORMATION OF THE BUTTERWORTH DESIGN (butter(N,wn) where N is the number of poles and wn is the normalized cutoff frequency)

3. A BILINEAR TRANSFORMATION OF THE CHEBYSHEV DESIGN 1 (cheby1(N,R,wn) where N is the number of poles, wn is the normalized cutoff frequency and R is the db level of the ripples in the passband)

4. A BILINEAR TRANSFORMATION OF THE CHEBYSHEV DESIGN 2 (cheby2(N,R,wn) where N is the number of poles, wn is the normalized cutoff frequency and R is the 1/db level of the ripples in the passband)





CODE: Note- cheby2 does not meet design criteria because it has an even number of poles so no pole at -1. Also off at wn.

cldigiir.m

w=[-pi:2*pi/256:pi];

[zc,pc,kc]=buttap(6),%assumes poles are on unit circle;

pc=pi*pc/4;%to get poles on circle with radius pi/4;

bc=poly(zc),ac=poly(pc);

[b1,a1]=impinvar(bc,ac);

b1=(pi/4)^6*b1;

zs=roots(b1),ps=roots(a1)

H1=freqz(b1,a1,w);

max(abs(H1))

subplot(222);zplane([],a1);subplot(221);plot(w/pi,abs(H1),'y');

xlabel('norm w');title('IMPULSE INVARIANT BUTTERWORTH');

axis([-1 1 0 1.2]);

[b2,a2]=butter(6,1/4);

H2=freqz(b2,a2,w);

subplot(224);zplane(b2,a2);subplot(223);plot(w/pi,abs(H2),'m');

xlabel('norm w');title('BILINEAR TRANSFORMATION BUTTERWORTH');

axis([-1 1 0 1.2]);

figure(2);

[b3,a3]=cheby1(6,.5,1/4);

H3=freqz(b3,a3,w);

subplot(222);zplane(b3,a3);subplot(221);plot(w/pi,abs(H3),'c');

xlabel('norm w');title('BILINEAR TRANSFORMATION CHEBY1');

axis([-1 1 0 1.2]);

[b4,a4]=cheby2(6,20,1/4);

H4=freqz(b4,a4,w);

subplot(224);zplane(b4,a4);subplot(223);plot(w/pi,abs(H4),'r');

xlabel('norm w');title('BILINEAR TRANSFORMATION CHEBY2');

axis([-1 1 0 1.2]);

figure(3);

subplot(121);

plot(w/pi,abs(H1),w/pi,abs(H2),...

w/pi,abs(H3),w/pi,abs(H4));

axis([0 1 0 1.2]);title('AMP COMPARISON');xlabel('norm w');

subplot(122);

plot(w/pi,20*log10(abs(H1)),w/pi,20*log10(abs(H2)),...

w/pi,20*log10(abs(H3)),w/pi,20*log10(abs(H4)));

axis([0 1 -200 50]);title('dB COMPARISON');xlabel('norm w');

FREQUENCY TRANSFORMATIONS OF LOWPASS FILTERS

Common practice to design one lowpass filter and then continue to transform to do other jobs:



We can go back to the s-plane, do an s transformation,

come back to z in usual fashion:



Highpass is easiest to see:

which has zero at 0 and goes to 1 for large W

Your Text has a table by which you can skip the s-plane work and make appropriate Z-Transformations:

i.e. Z-1 = G(z-1)

The transforms have to follow some rules

1. rational functions

2. inside unit circle of Z plane must go to inside circle of z plane

3. unit circles must map to unit circles,

which imples that |G(e-jw) | = 1, i.e. ALL_PASS

and that transformations must have form:



Lowpass to Lowpass:


Lowpass to Highpass:


See text for Lowpass to Bandpass and Lowpass to Bandstop.

AND IF YOU EXAMINE THE BUTTERWORTH ROUTINES IN

MATLAB, YOU WILL FIND THEY ARE ALL DONE THIS WAY.

BELOW IS THE CODE FOR A CLASS DEMO

BASED ON BUTTERWORTH ROUTINES

%[]=clbutt97(T,Qc) Qc=pi/4 and T=1 are suggestions for a start.

function []=clbutt97(T,Qc)

echo on;

kv=[1:6]';

N=6;

normQcT=Qc*T/pi

spoles=Qc*exp(j*(pi/(2*N))*(2*kv+N-1))

[szs,sps,K]=buttap(6);

testray=[spoles Qc*sps] %Compare results of BUTTAP

pause;clg;

v=[-1.2,1.2,-1.2,1.2];

axis('square');

axis(v);plot(spoles,'x');axis('square');axis(Qc*v);pause;

hold on;

a=poly(spoles);

b=[Qc^N];

% IMPULSE INVARIANT OF THE ENTIRE SYSTEM

[A,p,K] = residue(b,a);

zpoles=exp(p*T)

TA=T*A;

TK=T*K;

[ib,ia]=residuez(TA,zpoles,TK);

zzeros=roots(ib), pause;

clzplane(ib,ia);axis(v);pause;

ws=[0 Qc*T pi];

Hs=freqz(ib,ia,ws);testws=abs(Hs), pause; %Check at select frequencies

[ibc,iac]=impinvar(b,a,1/T);

bray=[ib' T*ibc'],tray=[ia' iac'], pause; % Compare to Matlab's Routine

%BILINEAR SUBSTITUTION OF THE ENTIRE SYSTEM

% FIRST FOR DIRECT FORMS

% ASSUME WE DON'T PREWARP

[Bb,Ba] = bilinear (b,a,1/T);

clg; clzplane(Bb,Ba);pause;

[Hs2]=freqz(Bb,Ba,ws);

testws=abs(Hs2),pause;

%USING PREWARPED BUTTERWORTH DESIGN

Qb=2/T*tan(T*Qc/2) %PREWARP

normQbT=Qb*T/pi

spoles=Qb*exp(j*(pi/(2*N))*(2*kv+N-1));

aw=poly(spoles);

bw=[Qb^N];

[Bbw,Baw] = bilinear (bw,aw,1/T);

clg; clzplane(Bbw,Baw);pause;

[Hs3]=freqz(Bbw,Baw,ws);

testws=abs(Hs3),pause;

%FREQUENCY RESPONSE FOR DIRECT FORMS

w=[-pi:2*pi/64:pi-2*pi/64];

%Impulse Invariant

[Hi] = freqz(ib,ia,w);pause;

%BILINEAR SUBSTITIO