# Filtering of signals

Document Sample

```					Filtering of signals

1

What is meant by a filter!
The DTFT is remembered again:

X (e jw ) 

n  

x[n]e  jwn 





1 x[n]  2

 X (e

jw

)e

jwn

dw

X[n] is expressed as a summation of sinusoids with scaled amplitude. Using a system with a frequency selective to these inputs, then it is possible to pass some frequencies and attenuate the others. Such a system is called a Filter.
2

The function of a filter is to remove unwanted parts of the signal, such as random noise, or to extract useful parts of the signal, such as the components lying within a certain frequency range.

Unfiltered signal or raw signal

Filtered signal

3

Example Choose frequency response of a system such that

H (e ) 
jw

1, 0,

for |w|≤wc for wc≤|w|≤π 0<w1<wc<w2< π

If x[n] =Acos(w1n) + Bcos(w2n),

y(n)=A|H(ejw1)|cos(w1n+Ө(w1))+ A|H(ejw2)|cos(w2n+Ө(w2))

Almost =0
y(n)=A|H(ejw)|cos(w1n+Ө(w1))

Which indicating the LPF effect of the LTI system
4

Design a HP digital filter that passes the 0.4rad/sec, and stops the 0.1 rad/sec frequency. h[n]=[ a b a]

H (e jw )  h[0]  h[1]e  jw  h[2]e  j 2 w
 a(1  e  j 2 w )  be  jw

 (2a cos(w)  b)e
H (e jw )  2a cos(w)  b

 e jw  e  jw   jw e  be  jw  2a    2   jw 
Ө(w)=-w

H (e j 0.1 )  2a cos(0.1)  b  0 H (e j 0.4 )  2a cos(0.4)  b  1
Solving for the two equations gives a=-6.76195, b=13.456335 y[n]=h[n]*x[n] y[n]=h[0]x[n]+h[1]x[n-1]+h[2]x[n-2]=ax[n]+bx[n-1]+ax[n-2]
5

y[n]=-6.76195(x[n]+x[n-2])+13.456335x[n-1]

If X[n]={cos(0.1n)+cos(0.4n)}u(n)
X1+X2 Output is almost equal to x2, the high frequency x1

Transient

Output of the filter

6

Classification of filters as analog or digital
Digital filters
A digital filter uses a digital processor to perform numerical calculations on sampled values of the signal. The processor may be a general-purpose computer such as a PC, or a specialized DSP (Digital Signal Processor) chip.

Analog filters
An analog filter uses analog electronic circuits made up from components such as resistors, capacitors and op amps to produce the required filtering effect. Such filter circuits are widely used in such applications as noise reduction, video signal enhancement, graphic equalizers in hi-fi systems, and many other areas.

7

Design of analog filters
We need to discuss some of the famous techniques for analog filter design due to two main reasons:

We need them as a prefilter or antialiasing filter before the
A/D conversion,

Some techniques for digital filter design are based on the
transformation of some analog techniques.
Ha(jΩ) 1+δp 1-δp Transition band Pass band ripples Stop band δs
8

Characteristics of a LPF

0

Ωp

Ωs

1-δp ≤ |Ha(jΩ)| ≤ 1+δp

0 ≤ |Ω| ≤ Ωp
Pass band Stop band

|Ha(jΩ)| ≤ δs

Ωs ≤ |Ω| ≤ ∞

k 

p s

Transition ratio, or selectivity parameter, it is larger than unity
Peak pass band ripple

ap=-20log10(1-δp) dB as=-20log10(δs) dB

Minimum stop band attenuation

Such filter is completely characterized by 9 Ωc, the 3dB point, Ωs, Ωp, δp, and δs.

Butterworth Approximation

1 H a ( j )  1  ( / c) 2 N
2
N=2 N=4

Ωc is the cut off frequency or the -3db cut ff frequency

N=10

It has a maximally flat magnitude at zero frequency. This clear 10 from (2N-1) differentiation of its function gives zeros.

MATLAB Butterworth functions [num,den] = BUTTER(N,Wn) designs an Nth order lowpass
digital Butterworth filter and returns the filter coefficients in length N+1 vectors (numerator) and (denominator). The coefficients are listed in descending powers of z. See MATLAB help for more Butterworth functions

%design of butterworth filter [num,den]=butter(4,1,'s'); w=0:0.1:5; h=freqs(num,den,w); gain=20*log10(abs(h)); subplot(2,1,1); plot(w,gain); grid; xlabel('Normalized frequency'); ylabel('Gain in dB');
11

Type 1 Chebyshev Approximation

1 H a ( j )  2 2 1   TN ( / p)
2

TN(Ω)= cos(Ncos-1Ω)
= cosh(Ncosh-1Ω)

|Ω| ≤ 1
|Ω| > 1

Type 1 Chebyshev filter. It is equiripple in the pass band and monotonically decreasing in the stop band

η2 represent the ripples in the pas band
CHEBY1 Chebyshev type I digital and analog filter design. [num,den] = CHEBY1(N,R,Wn) designs an Nth order lowpass digital Chebyshev filter with R decibels of ripple in the pass band. CHEBY1 returns the filter coefficients in length N+1 vectors (numerator) and (denominator).
12

N=3, Rp=1 dB, fp=1000Hz, wp=2*pi*fp
%design of type1 chebyshev filter [num,den]=cheby1(3,1,2000*pi,'s'); w=0:200:12000*pi; h=freqs(num,den,w); gain=20*log10(abs(h)); subplot(2,1,1); plot(w/(2*pi),gain); grid; ylabel('Gain in dB'); w=0:200:2400*pi; h=freqs(num,den,w); gain=20*log10(abs(h)); subplot(2,1,2); plot(w/(2*pi),gain); grid; xlabel('Frequency in Hz'); ylabel('Gain in dB');

13

Type 2 Chebyshev Approximation

H a ( j ) 
2

1  TN ( s /  p )  1    (  s / ) 
2

Type 2 Chebyshev filter. It is equiripple in the stop band and monotonically decreasing in the 2 pass band

CHEBY2 Chebyshev type II digital and analog filter design.

[B,A] = CHEBY2(N,R,Wn) designs an Nth order lowpass digital
Chebyshev filter with the stop band ripple R decibels down and stop band edge frequency Wn. CHEBY2 returns the filter coefficients in length N+1 vectors B (numerator) and A (denominator).

14

N=11, R=30 dB, fs=1000Hz, ws=2*pi*fs
Type2 Chebyshev filter

%design of type2 chebyshev filter [num,den]=cheby2(11,30, 2000*pi,'s'); w=0:200:12000*pi; h=freqs(num,den,w); gain=20*log10(abs(h)); subplot(2,1,1); plot(w/(2*pi),gain); grid; ylabel('Gain in dB');

15

Elliptic Approximation

1 H a ( j )  2 2 1   R N ( /  P )
2

Elliptic filter. It is equiripple in the stop band and in the pass band

RN(Ω) is a rational function satisfies the property RN(1/Ω)=1/RN(Ω)
ELLIP Elliptic or Cauer digital and analog filter design.

[B,A] = ELLIP(N,Rp,Rs,Wn) designs an Nth order lowpass digital
elliptic filter with Rp decibels of ripple in the passband and a stopband Rs decibels down. ELLIP returns the filter coefficients in length N+1 vectors B (numerator) and A (denominator).

16

N=7, pass band ripples 1 dB, Minimum stop band attenuation 30db cut off frequency 1000Hz
%design of elliptic filter [num,den]=ellip(7,1,30,2000*pi,'s'); w=0:200:12000*pi; h=freqs(num,den,w); gain=20*log10(abs(h)); subplot(2,1,1); plot(w/(2*pi),gain); grid; ylabel('Gain in dB'); w=0:200:2000*pi; h=freqs(num,den,w); gain=20*log10(abs(h)); subplot(2,1,2); plot(w/(2*pi),gain); grid; xlabel('Frequency in Hz'); ylabel('Gain in dB');

17

Classification of filters According to frequency response
H(ejw)

LPF
π

0

wc

w

H(ejw)
0 wc

HPF
π w

H(ejw)

BPF
w

0

π wc1 wc2

H(ejw)
0
18

BSF
π w

wc1 wc2

Pass band 0 ≤ w ≤ wc
LPF Stop band wc ≤ w ≤ π

HPF

Stop band 0 ≤ w ≤ wc Pass band wc ≤ w ≤ π Pass band wc1 ≤ w ≤ wc2

BPF

Stop band 0 ≤ w ≤ wc1, wc2 ≤ w ≤ π
Stop band wc1 ≤ w ≤ wc2

BSF

Pass band 0 ≤ w ≤ wc1, wc2 ≤ w ≤ π

19

Wc, wc1, and wc2 are called the cut off frequencies.

Design Of Digital Filters

20

1. A digital filter is programmable, i.e. its operation is
determined by a program stored in the processor's memory. This means the digital filter can easily be changed without affecting the circuitry (hardware). An analog filter can only be changed by redesigning the filter circuit.

2. Digital filters are easily designed, tested and
implemented on a general-purpose computer or workstation.

3. The characteristics of analog filter circuits (particularly
those containing active components) are subject to drift and are dependent on temperature. Digital filters do not suffer from these problems, and so are extremely stable with respect to both time and temperature. 21

4. Unlike their analog counterparts, digital filters can handle low
frequency signals accurately. As the speed of DSP technology continues to increase, digital filters are being applied to high frequency signals in the RF (radio frequency) domain, which in the past was the exclusive preserve of analog technology.

5. Digital filters are very much more versatile in their ability to
process signals in a variety of ways; this includes the ability of some types of digital filter to adapt to changes in the characteristics of the signal.

6. Fast DSP processors can handle complex combinations of filters
in parallel or cascade (series), making the hardware requirements relatively simple and compact in comparison with the equivalent analog circuitry.

22

1- Filter Characteristics Specification
Pass band w≤wp 1-δp≤|H(ejw)|≤1+δp δp Pass band deviation wp Pass band edge frequency Stop band ws≤w ≤π |H(ejw)|≤δs δs Stop band deviation ws Stop band edge frequency
23

Normalized LPF specs
k=Ωp/Ωs
k1  g A2  1

1 1 g 2

k, and k1 will be used to estimate the degree of IIR filter

1/A

24

Classification of filters according to impulse response length

Finite Impulse Response, FIR filters

Infinite Impulse Response, IIR filters

h[ n] 

M 1 n 0

a

n

u[ n]

h[ n]   a n u[ n]
n 0



25

2- Selection of filter type
FIR or IIR

1. FIR can have an exactly linear phase response. 2. FIR realized nonrecursively is always stable.
3. Quantization effects are less severe in FIR than in IIR. 4. FIR requires more coefficients for sharp cutoff than IIR. 5. Analog filters can be transformed into IIR. 6. FIR is easier to synthesize if CAD support is available.
An FIR system is always stable, But an IIR system may be stable or not, and it must be designed properly.
26

An originally stable IIR filter with precession coefficients may become unstable after implementation due to unavoidable quantization error in its coefficients. !!!

1 H ( z)  1  1.845z 1  0.850586z 1
Stable IIR filter

After quantization unstable IIR filter

H ( z) 
27

1 1  1.85z 1  0.85z 1

%Demonistration for the effect of quantization of filter coefficients L=100; %L is the length of the impulse response h[n] num=[1]; den=[1 -1.845 0.850586]; den2=[1 -1.85 0.85] [h1 t]=impz(num,den,L); This program draws subplot(2,1,1); the previous impulse stem(h1); ylabel('Amplitude'); response that shows xlabel('Time index n'); the effect of [h2 t]=impz(num,den2,L); quantization on the subplot(2,1,2); stem(h2); system stability. ylabel('Amplitude'); xlabel('Time index n');

28

Filter degree
1- IIR filter
Butterworth filter Chebyshev filter

log 10 (1 / k1) N  log 10 (1 / k )
Elliptic filter

cosh1 (1 / k1) N  cosh1 (1 / k )

2 log 10 ( 4 / k1) ,    0  2(  0) 5  15 (  0) 9  150 (  0)13 N  log 10 (1 /  )

0 

1 k'

2 1  k'





k'  1 k 2

2- FIR filter
N 
29

 20 log 10



 p  s  13



14 .6( ws  w p ) / 2

There are another approximation for very narrow pass band and very wide pass band Adawy

Example LPF with 1dB at wp=1kHz, and 40dB at ws=5kHz

 1 10 log 10  1 g 2 

   1  

g2=0.25895

 1  10 log 10  2   40 A 

A2=10000

1- Butterworth filter N 

log 10 (1 / k1) log 10 (1 / k )

N=3.281=4 N=2.6059=3 N=7

cosh1 (1 / k1) 2- Chebyshev filter N  cosh1 (1 / k )  20 log 10  p  s  13 3- FIR filter N  14 .6( ws  w p ) / 2 30





IIR Filter Design
Design using Analog Filters
1) Impulse Invariant Method (IIM) 2) Bilinear z-Transform (BZT) Method

31

1) Impulse Invariant Method (IIM)

In this method we require that the unit sample response h[n] of the desired causal digital transfer function H(z) be equal to the sampled version ha(t) sampled at uniform intervals of T seconds. h[n]=ha(nT) We have seen before that the freq. response H(ejw) of a sampled sequence ha(nT) is given by:

32

1  H (e jw )   H a ( j  jk sT ) T k   w / T

In the last equation put z=ejw, s=jΩ, and Ω=w/T will be s=(1/T)lin(z), so we can write this equation as follows:

1  H ( z )   H a ( s  j 2k ) T k  s (1/ T ) lin ( z )
Z transform of the sequence Laplace transform of the continuous impulse response
Take a look at this transformation

S=σ+jΩ

Z=eST=eσT ejΩT

S=(1/T)lin(z) or z=eST |z|=r=eσT Points on the unit circle r=1 in the z plane r<1 inside the unit circle of the z plane

Points on the jΩ axis means σ=0 in the s plane
33

For σ negative, left half of the s plane

jΩ
S Plane
3π/T

Im z

Z Plane

π/T

σ
-π/T

-1

1 Re z

-3π/T

Infinity of strips of width 2π/T in the s plane maps to the unit circle in the z plane. Left half of each strip maps 34 inside the unit circle

1  H (e jw )   H a ( j  jk sT ) T k   w / T
The resultant frequency response is given by the sum of the frequency responses of the original analog function and its shifted versions at ±2πk/T with the over all sum scaled by 1/T. So if Ha is band limited, there will be no aliasing, and if it is not band limited there will be aliasing as it was seen before.

So in Impulse Invariant method IIM
Stable analog filter maps to stable digital filter Frequency response of the digital filter is a multiple versions of the analog filter response centered at ±2πk/T, (sampling frequency).
35

2) Bilinear Transformation Method
Disadvantage of the previous method is the aliasing, and the many to one mapping where many strips in the s plane each of them maps to the entire z plane.
Bilinear Transformation Method

2  1  z 1 S   T  1  z 1 

   

 1  ST / 2  z   1  ST / 2 

H ( z )  Ha(s) s  2  1 z 1   
36

T  1 z 1   

Put T=2 to make things simple in the expression of S and Z

Put s=σ+jΩ

z

2

(1   ) 2   2  (1   ) 2   2

1    j z 1    j

For σ=0 , |z|=r=1
unit circle
Im z

S Plane

jΩ

Z Plane σ
Re z

For σ<0 LHP, |Z|<1
inside unit circle.

For σ>0 RHP, |Z|>1
outside unit circle.

One to one mapping
37

2  1  z 1 S   T  1  z 1 
2  w   tan  T 2

   

2 1  e  jw 2  w j   j tan   jw T 1 e T 2

w
π

W=2tan-1(ΩT/2)
0

Ω
+ve Ω maps into upper half of the unit circle, 0 to π

-π

-ve Ω maps into lower half of the unit circle, 0 to -π
38

The relation is highly nonlinear

w
π

W=2tan-1(ΩT/2)

w

0

Ω

|H(ejw)|

-π

|Ha(jΩ)| Ω

1- Given the digital filter specs, wp and ws, find there corresponding in the analog filter, use the relation between w and Ω 2- Design an analog filter that meets these specs. Find Ha(s) 3- Use bilinear transformation to find H(z) from Ha(s).

39

Example Design a digital LPF using bilinear transformation to satisfy the following: 1- Monotonic stop band and pass band, 2- 3.01dB cutoff frequency at 0.5π, and 15dB at w=0.75π
1- Monotonically decreasing magnitude means Butterworth filter

2- Calculation of corresponding analog frequencies considering T=1 :
Ωp=(2/T)tan(wp/2) = 2tan(0.5π/2)=2 Ωs=(2/T)tan(ws/2)=2tan(0.75π/2)=4.828 3- Calculation of
40

log 10 (1 / k1) N  log 10 (1 / k )

N=2

From tables of analog Butterworth filter :

H a (s) 

1 s 2  s2 2  4

Apply bilinear transformation

1  2z  z H ( z)  3.412  0.585 z 2

1

2

41

IIR digital filter design using MATLAB 1- BUTTER Butterworth digital and analog filter design.
[B,A] = BUTTER(N,Wn) designs an Nth order low pass digital Butterworth filter and returns the filter coefficients in length N+1 vectors B (numerator) and A (denominator). The coefficients are listed in descending powers of z. The cutoff frequency Wn must be 0.0 < Wn < 1.0, with 1.0 corresponding to half the sample rate. If Wn is a two-element vector, Wn = [W1 W2], BUTTER returns an order 2N band pass filter with pass band W1 < W < W2. [B,A] = BUTTER(N,Wn,'high') designs a high pass filter. [B,A] = BUTTER(N,Wn,'stop') is a band stop filter if Wn = [W1 W2]. When used with three left-hand arguments, as in [Z,P,K] = BUTTER(...), the zeros and poles are returned in length N column vectors Z and P, and the gain in scalar K. BUTTER(N,Wn,'s'), BUTTER(N,Wn,'high','s') and BUTTER(N,Wn,'stop','s') design analog Butterworth filters. In this case, Wn can be bigger than 1.0.

B( z ) b1  b2 z 1  ...... bN 1 z  N H ( z)   A( z ) 1  a 2 z 1  ........ a N 1 z  N 
42

[B,A]=butter(N,wn) [B,A]=cheby1(N,Rp,wn)

[z,p,k]=butter(N,wn) [z,p,k]=cheby1(N,Rp,wn)

[B,A]=cheby2(N,Rs,wn)
[B,A]=ellip(N,Rp,Rs,wn)

[z,p,k]=cheby2(N,Rs,wn)
[z,p,k]=ellip(N,Rp,Rs,wn)

Rp, ripples in the pass band, Rs, ripples in the stop band

The function freqz(B,A,w) can be used to plot the frequency response.

43

BUTTORD Butterworth filter order selection.
[N, Wn] = BUTTORD(Wp, Ws, Rp, Rs) returns the order N of the lowest order digital Butterworth filter that loses no more than Rp dB in the pass band and has at least Rs dB of attenuation in the stop band. Wp and Ws are the pass band and stop band edge frequencies, normalized from 0 to 1 (where 1 corresponds to pi radians/sample). For example, Low pass: Wp = .1, Ws = .2 High pass: Wp = .2, Ws = .1 Band pass: Wp = [.2 .7], Ws = [.1 .8] Band stop: Wp = [.1 .8], Ws = [.2 .7] BUTTORD also returns Wn, the Butterworth natural frequency (or, the "3 dB frequency") to use with BUTTER to achieve the specifications. [N, Wn] = BUTTORD(Wp, Ws, Rp, Rs, 's') does the computation for an analog filter, in which case Wp and Ws are in radians/second. When Rp is chosen as 3 dB, the Wn in BUTTER is equal to Wp in BUTTORD.

[N, Wn] = buttord(Wp, Ws, Rp, Rs) [N, Wn] = cheby1ord(Wp, Ws, Rp, Rs) [N, Wn] = cheby2ord(Wp, Ws, Rp, Rs) [N, 44 Wn] = ellipord(Wp, Ws, Rp, Rs)

Example
Determine the transfer function and plot the gain response of a 5th order elliptic IIR filter with a pass band edge at wp=0.4, pass band ripples Rp=0.5dB, and minimum stop band attenuation Rs=40dB. %elliptic IIR LP digital filter design Numerator polynomial N=5; 0.0528 0.0797 0.1295 0.1295 Rp=0.5; 0.0797 0.0528 Rs=40; wn=0.4; denominator polynomial [b,a]=ellip(5,0.5,40,0.4); 1.0000 -1.8107 2.4947 -1.8801 disp('Numerator polynomial'); 0.9537 -0.2336 disp(b); disp('denominator polynomial'); disp(a); 0.0528  0.0797 z 1  0.129 z 2  0.129 z 3  0.0797 z 4  0.0528 z 5 w=0:0.01/pi:pi; H ( z )  1  1.81 z 1  2.49 z 2  1.88 z 3  0.953 z 4  0.233 z 5 h=freqz(b,a,w); gain=20*log10(abs(h)); plot(w/pi,gain); grid; xlabel('normalized frequency'); ylabel('gain in dB'); 45 Adawy

Gain response of the previous filter

46

Pole zero representation of the previous filter
%elliptic IIR LP digital filter design (poles and zeros) N=5; Rp=0.5; Rs=40; wn=0.4; [b,a]=ellip(5,0.5,40,0.4); [z,p,k]=tf2zp(b,a); disp('Zeros are:'); disp(z); disp('Poles are:'); disp(p); disp('Gain constant is'); disp(k); zplane(b,a);
Zeros are: -1.0000 -0.3020 + 0.9533i -0.3020 - 0.9533i 0.0474 + 0.9989i 47 0.0474 - 0.9989i Poles are: 0.2787 + 0.8973i 0.2787 - 0.8973i 0.3812 + 0.6274i 0.3812 - 0.6274i 0.4909

Gain constant is 0.0528

Example
Determine the transfer function and plot the gain response of an 8th order Butterworth band pass digital filter with pass band edges 0.4 and 0.7.

%butterworth IIR BP digital filter N=8; M=N/2; w1=0.4; w2=0.7; wn=[w1 w2]; [b,a]=butter(M,wn); disp('numerator polynomial');disp(b); disp('denominator polynomial');disp(a); w=0.2*pi:0.01/pi:0.85*pi; h=freqz(b,a,w); gain=20*log10(abs(h)); plot(w/pi,gain); grid; 0.0186 xlabel('normalized frequency'); 0.0743 0 ylabel('gain in dB'); Denom pol. 48
1.0000 0.9780 1.9399 1.3386

0 -0.0743 0.0186

0

0.1114 Num pol.

0 -

1.6271

0.7349

0.5826 0.1386

Pole zero representation of the previous filter
%butterworth IIR BP digital filter design (poles and zeros) N=8; w1=0.4; w2=0.7; wn=[w1 w2]; [b,a]=butter(N/2,wn); [z,p,k]=tf2zp(b,a); disp('Zeros are:');disp(z); disp('Poles are:');disp(p); disp('Gain constant is');disp(k); zplane(b,a); Poles are: -0.4963 + 0.7058i -0.4963 - 0.7058i 0.2419 + 0.8060i 0.2419 - 0.8060i -0.2651 + 0.5715i -0.2651 - 0.5715i 0.0305 + 0.6027i 0.0305 - 0.6027i Zeros are: 1.0001 1.0000 + 0.0001i 1.0000 - 0.0001i 0.9999 -1.0001 -1.0000 + 0.0001i -1.0000 - 0.0001i -0.9999 Gain constant is 0.0186 Adawy

49

Design of FIR filters

FIR filters are characterized by linear phase IIR filters cannot produce a linear phase Linear phase is necessary when it is desired to maintain the shape of the wave

50

Linear phase=15 samples
phase

X=sin(pi/30)*n

X1=sin(pi/30)*(n-15)

w

y=sin(pi/10)*n y1=sin(pi/10)*(n-15)

+

+

51

The result is the same wave shifted by 15 samples

Nonlinear phase
phase

X=sin(pi/30)*n

X1=sin(pi/30)*(n-15)

w

y=sin(pi/10)*n
y1=sin(pi/10)*(n-2)

+

+

52 The wave shape is not conserved because delay is not linear

Realization of linear phase FIR filter

It is always possible to design an FIR transfer function with an exact linear phase, but nearly impossible to design an IIR linear phase one. 1- Design of linear phase FIR filter using odd length symmetric impulse response h[n]=h[N-n], 0≤n≤N
For N=8

H(z)=h[8]+h[7]z-1+h[6]z-2+h[5]z-3+h[4]z-4+h[3]z-5+h[2]z-6 +h[1]z-7+h[0]z-8
53

The impulse response is symmetric around the h[4] sample, so h[0]=h[8], h[1]=h[7], h[2]=h[6], h[3]=h[5]. So we can write: H(z)=h[0](1+z-8)+h[1](z-1+z-7)+h[2](z-2+z-6)+h[3](z-3+z-5) +h[4]z-4 H(z)=h[0]z-4(z4+z-4)+h[1]z-4(z3+z-3)+ h[2]Z-4(z2+z-2)+ h[3]z-4(z1+z-1) +h[4]z-4 H(z)=z-4{h[0](z4+z-4)+h[1](z3+z-3)+ h[2](z2+z-2)+ h[3](z1+z-1) +h[4]} H(ejw)=e-jw4{2h[0]cos(4w)+2h[1]cos(3w)+2h[2]cos(2w)+ 2h[3]cos(w)+h[4]} Ө(w)=-4w, which is linear, GD=4 which is constant

N / 2  jw  jwN / 2 H (e )  e   a[n] cos(wn)  n 0 
54

Where a[0]=h[N/2],

a[n]=2h[N/2 – n],

1≤n≤N/2

2- Design of linear phase FIR filter using even length symmetric impulse response
For N=7

h[0]=h[7], h[1]=h[6], h[2]=h[5], h[3]=h[4]. H(z)=h[0](1+z-7)+h[1](z-1+z-6)+h[2](z-2+z-5)+h[3](z-3+z-4) H(z)=z-7/2{h[0](z7/2+z-72)+h[1](z5/2+z-5/2)+ h[2](z3/2+z-3/2)+ h[3](z1/2+z-1/2)} H(ejw)=e-jw7/2{2h[0]cos(7w/2)+2h[1]cos(5w/2)+ 2h[2]cos(3w/2)+ 2h[3]cos(w/2)} Ө(w)=-7w/2, which is linear, GD=7/2 which is constant

( N 1) / 2  jw  jwN / 2 H (e )  e   b[n] cos(w(n  0.5))  n1 
55

Where

b[n]=2h[(N+1)/2 – n],

1≤n≤(N+1)/2

Example

H(z)=(1/6){0.5+z-1+z-2+z-3+z-4+z-5+0.5z-6}
Impulse response

Frequency response

Phase response, it is linear in the pass band

56

3- Design of linear phase FIR filter using odd length antisymmetric impulse response

h[n]=-h[N-n],

0≤n≤N, h[N/2]=0

For N=8

H(z)=z-4{h[0](z4-z-4)+h[1](z3-z-3)+ h[2](z2-z-2)+ h[3](z1+z-1)} H(ejw)=e-j(4w-π/2){2h[0]sin(4w)+2h[1]sin(3w)+2h[2]sin(2w)+ 2h[3]sin(w)}

Ө(w)=-4w-π/2, which is linear,

GD=4 which is constant

N / 2  jw  j ( wN / 2 / 2 ) H (e )  e   c[n] sin(wn)  n1 
Where
57

c[n]=2h[N/2 – n],

1≤n≤N/2

4- Design of linear phase FIR filter using even length antisymmetric impulse response
For N=7

H(z)=z-7/2{h[0](z7/2-z-72)+h[1](z5/2-z-5/2)+ h[2](z3/2-z-3/2)+ h[3](z1/2-z-1/2)}
H(ejw)=e-(jw7/2-π/2){2h[0]sin(7w/2)+2h[1]sin(5w/2)+ 2h[2]sin(3w/2)+ 2h[3]sin(w/2)} Ө(w)=-7w/2 –π/2, which is linear, constant GD=7/2 which is

( N 1) / 2  jw  j ( wN / 2  / 2 ) H (e )  e   d [n] sin(w(n  0.5))  n 1 
Where
58

d[n]=2h[(N+1)/2 – n],

1≤n≤(N+1)/2

The four cases of symmetry required for linear phase realization
Even symmetry

even antisymmetry

Odd symmetry

Odd antisymmetry

Point of symmetry
59

Design of FIR filters using Windows
Simply truncate the response of an IIR filter to a desired length. Impulse response of the FIR filter
hd[n] h[n]= 0 Impulse response of the IIR filter N1≤n≤N2 otherwise
w[n]

h[n]=hd[n].w[n] 1 w[n]=
60

N1≤n≤N2 otherwise

0

From the properties of the Fourier transform
h[n]=hd[n].w[n] F.T


H(ejw)=Hd(eiw)*W(ejw)

1 jw H (e )  2
Hd(eiw)



H d (e jw )W (e j ( w ) )d 
W(ejw)

H(ejw)

wc

π

*
w

=
w w
4π/N

This convolution produced a smeared version of the ideal LPF frequency response. The wider the main lobe of W, the more spread is H. So, there is a trade off between large N to minimize smearing, and small N to get a reasonable implementation of the filter. 61 Adawy

Different window shapes
Much work has been done to adjust w[n] to satisfy certain main lobe width and relative side lobes levels requirements. These windows are: Rectangle, Bartlett, Hanning, Hamming, Blackman, Kaiser, and more.

62

Window Function
Rectangle Window

1 w[n]=
0

0≤n≤M
otherwise

Hanning Window

1  2n  w[n]  1  cos  2  2 M  1 
Hamming Window

-M ≤ n ≤ M

 2n  w[n]  0.54  0.46 cos   2M  1 

-M ≤ n ≤ M

Blackman Window

 2n   4n  w[n]  0.42  0.5 cos   0.08   2M  1   2M  1 
63

-M ≤ n ≤ M

Characteristics of different windows
Type of window Rectangle Hanning

Main lobe width
4π/(2M+1) 8π/(2M+1)

Relative side lobe level
13.3dB 31.5dB

Minimum stop band atten.
20.9dB 43.9dB

Transition band width 0.92π/M 3.11π/M

Hamming
Blackman

8π/(2M+1)

42.7dB

54.5dB
75.3dB

3.32π/M
5.56π/M

12π/(2M+1) 58.1dB

Transition band width (wp-ws)=c/M, where c is a constant depends on type of window as in the table. This relation can be used to estimate the order of the filter. 64 Adawy

Example
1

|H(eiw) 1 0 0≤w≤wc wc≤w≤π

Hd(ejw)=

Wc=π/2

π

w

The impulse response of the desired filter is obtained by calculating the inverse DTFT as follows:

1 h[n]  2

1 jw jwn  Hd (e )e dw  2 



wc

 wc

e

jwn

dw 

sin( wc n) -∞ ≤ n ≤∞ n

This response will be windowed for -M ≤n ≤M, and shifting by M samples to make it realizable (causal).

h' [ n] 

sin( wc (n  M ))  (n  M )

0≤ n ≤N-1, N=2M+1 Freq. response Adawy

65

N=51, wc=0.25pi, Rectangle window

66

Design of FIR filters based on frequency response sampling
Start by the desired frequency response Hd(ejw)

Uniformly sample it at N equally spaced points w k=2πk/N, k=0, 1, 2, …N-1

These samples represent N point DFT H[k]

Take IDFT of these N points to get the impulse response h[n]
67

FIR digital filter using MATLAB
MATLAB includes 5 window functions that can be used to design FIR filter using the window method, these windows are:

W=blackman(N), W=hamming(N), w=hanning(N), w=ch ebwin(N,Rs), w=kaiser(N,beta). N is the number of samples, beta, Rs are coefficients related transition bands. b=fir1(L,wn),
generates a length N=L+1 vector b that represent the impulse response coefficients of the LPF FIR filter. Hamming window is considered if no window is mentioned. Filter type (low, high, band stop or band stop). Also a type of the window can be included like this form;

b=fir1(L,wn, ‘filter type’, window)

b=fir1(L,wn, ‘stop’, kaiser)
68

Example
Design a LP FIR filter using kaiser, and hanning windows with the following specifications: wp=0.3pi, ws=0.4pi, As=50dB.

0.1102(As-8.7)

for As>50

β=

0.5842(As-21)0.4 + 0.07886(As-21) for 21 < As < 50

0

for As<21

As is the stop band attenuation=50, so, β=4.55126 Δf=(ws-wp)/2pi = 0.05

N

As  7.95  1 For As>21 14 .36 f

N=61

0.9222 N  1 For As<21 f 69

wp=0.3; kw=kaiser(61,4.55126); b=fir1(60,0.3,kw); w=0:pi/256:pi; h=freqz(b,1,w); mag=20*log10(abs(h)); subplot(2,1,1); plot(w/pi,mag); AXIS([0 1 -80 0]); ham=hamming(61); b2=fir1(60,0.3,ham); w=0:pi/256:pi; h2=freqz(b2,1,w); mag2=20*log10(abs(h2)); subplot(2,1,2); plot(w/pi,mag2); AXIS([0 1 -80 0]);
70

%design of FIR LPF N=61; beta=4.55126; wp=0.3; kw=kaiser(61,4.55126); b=fir1(60,0.3,kw); w=0:pi/256:pi; h=freqz(b,1,w); mag=20*log10(abs(h)); %fitering a sequence of data n=0:1:100; b=fir1(60,0.3,kw); x1=cos(0.1*pi*n); x2=cos(0.4*pi*n); y=filter(b,1,x1); subplot(3,1,1); plot(n,y,'r',n,x1,'k'); y=filter(b,1,x2); subplot(3,1,2); plot(n,y,'r',n,x2,'k'); y=filter(b,1,x1+x2); subplot(3,1,3); plot(n,y,'r',n,x1+x2,'k'); 71 ylabel('Amplitude'); xlabel('Time index');

Response of previous filter for different inputs