Digital Signal Processing-chapter 16: MATLAB Programs

MATLAB Programs

Chapter 16 16.1

INTRODUCTION

MATLAB stands or MATrix LABoratory. It is a technical computing environment or high perormance numeric computation and visualisation. It integrates numerical analysis, matrix computation, signal processing and graphics in an easy-to-use environment, where problems and solutions are expressed just as they are written mathematically, without traditional programming. MATLAB allows us to express the entire algorithm in a ew dozen lines, to compute the solution with great accuracy in a ew minutes on a computer, and to readily manipulate a three-dimensional display o the result in colour. MATLAB is an interactive system whose basic data element is a matrix that does not require dimensioning. It enables us to solve many numerical problems in a raction o the time that it would take to write a program and execute in a language such as FORTRAN, BASIC, or C. It also eatures a amily o application specifc solutions, called toolboxes . Areas in which toolboxes are available include signal processing, image processing, control systems design, dynamic systems simulation, systems identifcation, neural networks, wavelength communication and others. It can handle linear, non-linear, continuous-time, discrete-time, multivariable and multirate systems. This chapter gives simple programs to solve specifc problems that are included in the previous chapters. All these MATLAB programs have been tested under version 7.1 o MATLAB and version 6.12 o the signal processing toolbox.

16.2

REPRESENTATION OF BASIC SIGNALS

MATLAB programs or the generation o unit impulse, unit step, ramp, exponential, sinusoidal and cosine sequences are as ollows.

% Program or the generation o unit impulse signal clc;clear all;close all; t522:1:2; y5[zeros(1,2),ones(1,1),z [zeros(1,2),ones(1,1),zeros(1,2)];s eros(1,2)];subplot(2,2,1 ubplot(2,2,1);stem(t,y); );stem(t,y);

816

Digital Signal Processing

ylabel(‘Amplitude --.’); xlabel(‘(a) n --.’);

% Program or the generation o unit step sequence [u(n) 2 u(n 2 N] n5input(‘enter the N value’); t50:1:n21; y15ones(1,n);subplot(2,2,2); stem(t,y1);ylabel(‘Amplitude --.’); xlabel(‘(b) n --.’);

% Program or the generation o ramp sequence n15input(‘enter the length o ramp sequence’); t50:n1; subplot(2,2,3);stem(t,t);ylabel(‘Amplitude --.’); xlabel(‘(c) n --.’);

% Program or the generation o exponential sequence n25input(‘enter the length o exponential sequence’); t50:n2; a5input(‘Enter the ‘a’ value’); y25exp(a*t);subplot(2,2,4); stem(t,y2);ylabel(‘Amplitude --.’); xlabel(‘(d) n --.’);

% Program or the generation o sine sequence t50:.01:pi; y5sin(2*pi*t);gure(2); subplot(2,1,1);plot(t,y);ylabel(‘Amplitude --.’); xlabel(‘(a) n --.’);

% Program or the generation o cosine sequence t50:.01:pi; y5cos(2*pi*t); subplot(2,1,2);plot(t,y);ylabel(‘Amplitude --.’); xlabel(‘(b) n --.’);

As an example, enter the N value 7 enter the length o ramp sequence 7 enter the length o exponential sequence 7 enter the a value 1

Using the above MATLAB programs, we can obtain the waveorms o the unit impulse signal, unit step signal, ramp signal, exponential signal, sine wave signal and cosine wave signal as shown in Fig. 16.1.

817

MATLAB Programs

e d u t i l p m A

1

1

0.8

0.8

0.6

e d u t i l p m A

0.4 0.2

0.6 0.4 0.2 0

0

−2

−1

0

1

2

0

2

4

n

6

n

(a)

(b) 1

7 6

0.8

5 e d u t i l p m A

4

e d u t i l p m A

3 2 1 0

0.6 0.4 0.2 0

0

2

4

6

8

0

2

4

6

8

n

n

(c)

(d)

1 0.5 0 e d u t i l p m A

− 0.5 −1

0

0.5

1

1.5

2

2.5

3

3.5

2.5

3

3.5

n

(e) 1 0.5 e d u t i l p m A

0

− 0.5 −1

0

0.5

1

1.5 (f)

2 n

Fig. 16.1 Repr Repres esen enta tati tion on of Basi Basicc Sign Signal alss ( a a) Unit Impulse Signal (b) Unit-step Signal ( c c) Ramp Signal ( d d) Exponential Signal ( e e) Sinewave Signal ( f )Cosine Wave Signal

818


DISCRETE CONVOLUTION

16.3

16.3.1 Linear Convolution

Algorithm 1. Get two signals x signals x((m)and h( p)in p)in matrix orm 2. The convolved signal is denoted as y as y((n) 3. y( y(n)is given by the ormula ∞

y( y(n) 5

∑ [ x(k ) h(n − k )]

where n50 to m 1 p 2 1

k =−∞

4. Stop

% Program or linear convolution o the sequence x5 [1, [1, 2] and and h5 [1, [1, 2, 4] clc; clear all; close all; x5input(‘enter the 1st sequence’); h5input(‘enter the 2nd sequence’); y5conv(x,h); gure;subplot(3,1,1); stem(x);ylabel(‘Amplitude --.’); xlabel(‘(a) n --.’); subplot(3,1,2); stem(h);ylabel(‘Amplitude --.’); xlabel(‘(b) n --.’); subplot(3,1,3); stem(y);ylabel(‘Amplitude --.’); xlabel(‘(c) n --.’); disp(‘The resultant signal is’);y

As an example, enter the 1st sequence [1 2] enter the 2nd sequence [1 2 4] The resultant signal is y51 4 8 8

Figure 16.2 shows the discrete input signals x signals x((n)and h(n)and the convolved output signal y signal y((n). 2 e d u t i l p m A

1.5 1 0.5 0

1

1.1

1.2

1.3

1.4

1.5

1.6

(a)

Fig. 16.2

(Contd.)

1.7

1.8 n

1.9

2

819

MATLAB Programs

4 e d u t i l p m A

3 2 1 0

1

1.2

1.4

1.6

1.8

2

2.2

2.4

(b)

2.6

2.8

3

n

8 e d u t i l p m A

6 4 2 0 1

1.5

2

2.5

3

(c)

3.5 n

Fig. 16.2 Discrete Linear Convolution Convolution

16.3.2 Circular Convolution

% Program or Computing Circular Convolution clc; clear; a = input(‘enter the sequence x(n) = ’); b = input(‘enter the sequence h(n) = ’); n1=length(a); n2=length(b); N=max(n1,n2); x = [a zeros(1,(N-n1))]; or i = 1:N k = i; or j = 1:n2 H(i,j)=x(k)* b(j); k = k-1; i (k == 0) k = N; end end end y=zeros(1,N); M=H’; or j = 1:N or i = 1:n2 y(j)=M(i,j)+y(j); end end disp(‘The output sequence is y(n)= ‘); disp(y);

4

820


stem(y); title(‘Circular Convolution’); xlabel(‘n’); ylabel(‚y(n)‘);

As an Example, enter the sequence x(n) = [1 2 4] enter the sequence h(n) = [1 2] The output sequence is y(n)= 9 4 8

% Program or Computing Circular Convolution with zero padding clc; close all; clear all; g5input(‘enter the rst sequence’); h5input(‘enter the 2nd sequence’); N15length(g); N25length(h); N5 max(N1,N2); N35N12N2; %Loop or getting equal length sequence i(N350) h5[h,zeros(1,N3)]; else g5[g,zeros(1,2N3)]; end %computation o circular convolved sequence or n51:N, y(n)50; or i51:N, j5n2i11; i(j550) j5N1j; end y(n)5y(n)1g(i)*h(j); end end disp(‘The resultant signal is’);y

As an example, enter the rst sequence [1 2 4] enter the 2nd sequence [1 2] The resultant signal is y51 4 8 8

16.3.3 Overlap Save Method and Overlap Add method

% Program or computing Block Convolution using Overlap Save Method Overlap Save Method x=input(‘Enter the sequence x(n) = ’);

MATLAB Programs

821

h=input(‘Enter the sequence h(n) = ’); n1=length(x); n2=length(h); N=n1+n2-1; h1=[h zeros(1,N-n1)]; n3=length(h1); y=zeros(1,N); x1=[zeros(1,n3-n2) x zeros(1,n3)]; H=t(h1); or i=1:n2:N y1=x1(i:i+(2*(n3-n2))); y2=t(y1); y3=y2.*H; y4=round(it(y3)); y(i:(i+n3-n2))=y4(n2:n3); end disp(‘The output sequence y(n)=’); disp(y(1:N)); stem(y(1:N)); title(‘Overlap Save Method’); xlabel(‘n’); ylabel(‘y(n)’); Enter the sequence x(n) = [1 2 -1 2 3 -2 -3 -1 1 1 2 -1] Enter the sequence h(n) = [1 2 3 -1] The output sequence y(n) = 1 4 6 5 2 11 0 -16 -8 3 8 5 3 -5 1

%Program or computing Block Convolution using Overlap Add Method x=input(‘Enter the sequence x(n) = ’); h=input(‘Enter the sequence h(n) = ’); n1=length(x); n2=length(h); N=n1+n2-1; y=zeros(1,N); h1=[h zeros(1,n2-1)]; n3=length(h1); y=zeros(1,N+n3-n2); H=t(h1); or i=1:n2:n1 i i<=(n1+n2-1) x1=[x(i:i+n3-n2) zeros(1,n3-n2)]; else x1=[x(i:n1) zeros(1,n3-n2)]; end x2=t(x1); x3=x2.*H; x4=round(it(x3)); i (i==1)

822


y(1:n3)=x4(1:n3); else y(i:i+n3-1)=y(i:i+n3-1)+x4(1:n3); end end disp(‘The output sequence y(n)=’); disp(y(1:N)); stem((y(1:N)); title(‘Overlap Add Method’); xlabel(‘n’); ylabel(‘y(n)’);

As an Example, Enter the sequence x(n) = [1 2 -1 2 3 -2 -3 -1 1 1 2 -1] Enter the sequence h(n) = [1 2 3 -1] The output sequence y(n) = 1 4 6 5 2 11 0 -16 -8 3 8 5 3 -5 1

16.4

DISCRETE CORRELATION

16.4.1 Crosscorrelation

Algorithm 1. Get two signals x signals x((m)and h( p)in p)in matrix orm 2. The correlated signal is denoted as y as y((n) 3. y( y(n)is given by the ormula ∞

y( y(n) 5

∑ [ x(k ) h(k − n)] k =−∞

where n52 [max (m (m, p) p)2 1] to [max (m (m, p) p)2 1] 4. Stop

% Program or computing cross-correlation o the sequences x5 [1, 2, 3, 4] and h5 [4, 3, 2, 1] clc; clear all; close all; x5input(‘enter the 1st sequence’); h5input(‘enter the 2nd sequence’); y5xcorr(x,h); gure;subplot(3,1,1); stem(x);ylabel(‘Amplitude --.’); xlabel(‘(a) n --.’); subplot(3,1,2); stem(h);ylabel(‘Amplitude --.’); xlabel(‘(b) n --.’); subplot(3,1,3); stem(fiplr(y));ylabel(‘Amplitude --.’);

823

MATLAB Programs

xlabel(‘(c) n --.’); disp(‘The resultant signal is’);fiplr(y)

As an example, enter the 1st sequence [1 2 3 4] enter the 2nd sequence [4 3 2 1] The resultant signal is y51.0000 4.0000 10.0000 20.0000 25.0000 24.0000 16.0000 ↑

Figure 16.3 shows the discrete input signals x signals x((n)and h(n)and the cross-correlated output signal y signal y((n).

4 e d u t i l p m A

3 2 1 0

1

1.5

2

2.5

3

(a)

3.5

4

3.5

4

6

7

n

4 e d u t i l p m A

3 2 1 0

1

1.5

2

2.5

3

(b)

n

30 e d u t i l p m A

20 10 0

1

2

3

4

5

(c)

Fig. 16.3 Discrete Cross-correlation Cross-correlation

16.4.2 Autocorrelation

Algorithm 1. Get the signal x signal x((n)o length N length N in in matrix orm 2. The correlated signal is denoted as y as y((n) 3. y( y(n)is given by the ormula ∞

y( y(n) 5

∑ [ x(k ) x(k − n)] k =−∞

where n52( N N 2 1) to ( N 2 1)

n

824


% Program or computing autocorrelation unction x5input(‘enter the sequence’); y5xcorr(x,x); gure;subplot(2,1,1); stem(x);ylabel(‘Amplitude --.’); xlabel(‘(a) n --.’); subplot(2,1,2); stem(fiplr(y));ylabel(‘Amplitude --.’); xlabel(‘(a) n --.’); disp(‘The resultant signal is’);fiplr(y)

As an example, enter the sequence [1 2 3 4] The resultant signal is y54 11 20 30 20 11 4 ↑

Figure 16.4 shows the discrete input signal x( x(n)and its auto-correlated output signal y signal y((n). 4 e d u t i l p m A

3 2 1 0

1

1.5

2.5

2

1

2

3

3.5

( )

(a) 30 25 e 20 d 15 u t i l p 10 m 5 A 0

3

4

4 n

5

(b) (b) y (n) (n)

6

7 n

Fig. 16.4 Discrete Auto-correlation Auto-correlation

16.5

STABILITY STABILITY TEST

% Program or stability test clc;clear all;close all; b5input(‘enter the denominator coecients o the lter’); k5poly2rc(b); 5fiplr(k); knew s5all(abs(knew)1); i(s55 1) disp(‘“Stable system”’);

MATLAB Programs

825

else disp(‘“Non-stable system”’); end

As an example, enter the denominator coecients o the lter [1 21 .5] “Stable system”

SAMPLING THEOREM

16.6

The sampling theorem can be understood well with the ollowing example. discrete-time Example 16.1 Frequency analysis of the amplitude modulated discrete-time signal x( x(n)5cos 2 p f 1n 1 cos 2p f 2n where f 1

1 =

128

and f 2

5 =

128

modulates the amplitude-modulated signal is

xc(n)5cos 2p f c n where f where f c550/128. The resulting amplitude-modulated signal is xam(n)5 x( x(n) cos 2p f c n Using MATLAB program, program, (a) (b) (c)

sketch the signals x( x (n), x ), xc(n) and x xam(n), 0 # n # 255 compute and sketch the 128-point DFT of the signal x xam(n), 0 # n # 127 compute and sketch the 128-point DFT of the signal x am(n), 0 # n # 99

Solution

% Program Solution or Section (a) clc;close all;clear all; 151/128;255/128;n50:255;c550/128; x5cos(2*pi*1*n)1cos(2*pi*2*n); xa5cos(2*pi*c*n); xamp5x.*xa; subplot(2,2,1);plot(n,x);title(‘x(n)’); xlabel(‘n --.’);ylabel(‘amplitude’); subplot(2,2,2);plot(n,xc);title(‘xa(n)’); xlabel(‘n --.’);ylabel(‘amplitude’); subplot(2,2,3);plot(n,xamp); xlabel(‘n --.’);ylabel(‘amplitude’); %128 point DFT computation2solution or Section (b) n50:127;gure;n15128; 151/128;255/128;c550/128; x5cos(2*pi*1*n)1cos(2*pi*2*n); xc5cos(2*pi*c*n); xa5cos(2*pi*c*n);

(Contd.) Contd.)

826


2

1 e d u t i l p m A

0

−1

−2

0

100 (i)

Fig. 16.5(a)

200

300 n

(i) Modulating Signal x (n)

1

0.5 e d u t i l p m A

0

− 0.5

−1 0

100

200 (ii)

Fig. 16.5(a)

300

n

(ii) Carrier Signal and

2

1 e d u t i l p m A

0

−1

−2

0

100

200 (iii)

Fig. 16.5(a)

300

n

(iii) Amplitude Modulated Signal

(Contd.) Contd.)

MATLAB Programs

827

25

20

15

10

e d u t i l p m A

5

0

−5

− 10

0

20

40

60

80

100

120

140

n

Fig. 16.5(b)

128-point DFT of the Signal x am (n) , , 0 # n # 127

35

30

25

20

e d 15 u t i l p m A 10

5

0

−5

0

20

40

60

80

100

120 n

Fig. 16.5(c)

128-point DFT of the Signal x am (n) , , 0 # n # 99

140

828


5t(xamp,n1); xamp5x.*xa;xam stem(n,xam);title(‘xamp(n)’);xlabel(‘n --.’); ylabel(‘amplitude’); %128 point DFT computation2solution or Section (c) n50:99;gure;n250:n121; 151/128;255/128;c550/128; x5cos(2*pi*1*n)1cos(2*pi*2*n); xc5cos(2*pi*c*n); xa5cos(2*pi*c*n); xamp5x.*xa; or i51:100, xamp1(i)5xamp(i); end 5t(xamp1,n1); xam stem(n2,xam);title(‘xamp(n)’);xlabel(‘n --.’);ylabel(‘amplitude’);

(a)Modulated signal x( x(n), carrier signal xa(n) and amplitude modulated signal xam(n) are shown in Fig. 16.5(a). Fig. 16.5 (b) shows the 128-point DFT o the signal xam(n) or 0 # n # 127 and Fig. 16.5 (c) shows the 128-point DFT o the signal x signal xam(n), 0 # n # 99.

16.7

FAST FOURIER TRANSFORM

Algorithm 1. Get the signal x signal x((n)o length N length N in in matrix orm 2. Get the N the N value value 3. The transormed signal is denoted as N −1

x( k ) = ∑ x( n)e

− j

2p N

nk

for 0 ≤ k ≤ N −1

n= 0

\ % \

Program or computing discrete Fourier transorm

clc;close all;clear all; x5input(‘enter the sequence’); n5input(‘enter the length o t’); X(k)5t(x,n); stem(y);ylabel(‘Imaginary axis --.’); xlabel(‘Real axis --.’); X(k)

As an example, enter the sequence [0 1 2 3 4 5 6 7] enter the length o t 8 X(k)5 Columns 1 through 4 28.0000 24.000019.6569i 24.0000 14.0000i 24.0000 1 1.6569i Columns 5 through 8 24.0000 24.0000 21.6569i 24.0000 24.0000i 24.0000 29.6569i

MATLAB Programs

829

The eight-point decimation-in-time ast Fourier transorm o the sequence x( x(n)is computed using MATLAB program and the resultant output is plotted in Fig. 16.6. 10 8

6

4

2

s i x a y r a n i g a m I

0

−2 −4 −6 −8

− 10

−5

0

5

10

15

20

25

Real axis

Fig. 16.6

16.8

Fast Fourier Transform

BUTTERWORTH ANALOG FILTERS

16.8.1 Low-pass Filter

Algorithm 1. 2. 3. 4. 5. 6.

Get the passband and stopband ripples Get the passband and stopband edge requencies Get the sampling requency Calculate the order order o the flter flter using using Eq. 8.46 Find the flter coefcients Draw the magnitude and phase responses.

% Program or the design o Butterworth analog low pass flter clc; close all;clear ormat long rp5input(‘enter rs5input(‘enter wp5input(‘enter ws5input(‘enter s5input(‘enter

all; the the the the the

passband stopband passband stopband sampling

ripple’); ripple’); req’); req’); req’);

30

830


w152*wp/s;w252*ws/s; [n,wn]5buttord(w1,w2,rp,rs,’s’); [z,p,k]5butter(n,wn); [b,a]5zp2t(z,p,k); [b,a]5butter(n,wn,’s’); 50:.01:pi; w [h,om]5reqs(b,a,w); 520*log10(abs(h)); m an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel(‘Gain in dB --.’);xlabel(‘(a) Normalised requency --.’); subplot(2,1,2);plot(om/pi,an); xlabel(‘(b) Normalised requency --.’); ylabel(‘Phase in radians --.’);

As an example, enter enter enter enter enter

the the the the the

passband stopband passband stopband stopband

ripple ripple req req req

0.15 60 1500 3000 7000

The amplitude and phase responses o the Butterworth low-pass analog flter are shown in Fig. 16.7. 50 0

− 50 B d

− 100

ni ni

− 150 a G

− 200 − 250

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1

Normalised frequency (a) 4 2 s ai

n

0 ar

d

e

ni

−2 a

s h P

−4

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Normalised frequency (b)

Fig. 16.7 Butterworth Low-pass Low-pass Analog Filter ( a a) Amplitude Response and (b) Phase Response

MATLAB Programs

831

16.8.2 High-pass Filter

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge requencies 3. Get the sampling requency 4. Calculate the order order o the flter flter using using Eq. 8.46 5. Find the flter coefcients 6. Draw the magnitude and phase responses.

% Program or the design o Butterworth analog high—pass flter clc; close all;clear all; ormat long rp5input(‘enter the passband ripple’); rs5input(‘enter the stopband ripple’); wp5input(‘enter the passband req’); ws5input(‘enter the stopband req’); s5input(‘enter the sampling req’); w152*wp/s;w252*ws/s; [n,wn]5buttord(w1,w2,rp,rs,’s’); [b,a]5butter(n,wn,’high’,’s’);

50:.01:pi; w [h,om]5reqs(b,a,w);

520*log10(abs(h)); m an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel(‘Gain in dB --.’);xlabel(‘(a) Normalised requency --.’); subplot(2,1,2);plot(om/pi,an); xlabel(‘(b) Normalised requency --.’); ylabel(‘Phase in radians --.’);

As an example, enter the passband ripple

0.2

enter the stopband ripple

40

enter the passband req

2000

enter the stopband req

3500

enter the sampling req

8000

The amplitude and phase responses o Butterworth high-pass analog flter are shown in Fig. 16.8.

832


100 0 B

− 100 ni

− 200

in

d

a G

− 300 − 400 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1

Normalised frequency

(a) 4 2 s n ai

0 d ar ni

−2 e s a h P

−4 0

0.1

0.2

0.3

0.4

0.5

0.6

(b)

0.7

0.8


Fig. 16.8 Butterworth High-pass High-pass Analog Filter ( a) Amplitude Response and (b) Phase Response

16.8.3 Bandpass Filter

Algorithm 1. 2. 3. 4. 5. 6.


% Program or the design o Butterworth analog Bandpass flter clc; close all;clear all; ormat long rp5input(‘enter the passband ripple...’); rs5input(‘enter the stopband ripple...’); wp5input(‘enter the passband req...’); ws5input(‘enter the stopband req...’); s5input(‘enter the sampling req...’); w152*wp/s;w252*ws/s;

MATLAB Programs

833

[n]5buttord(w1,w2,rp,rs); wn5[w1 w2]; [b,a]5butter(n,wn,’bandpass’,’s’);

50:.01:pi; w [h,om]5reqs(b,a,w);


As an example, enter the passband ripple...

0.36

enter the stopband ripple...

36

enter the passband req...

1500

enter the stopband req...

2000

enter the sampling req...

6000

The amplitude and phase responses o Butterworth bandpass analog flter are shown in Fig. 16.9. 200 0

− 200 B d

− 400 ni ni

− 600 a G

− 800 − 1000 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1


(a) 4 2 s n ia

0 d a r ni

−2 e s a h P

−4

0

0.1

0.2

0.3

0.4

0.5 (b)

0.6

0.7

0.8


Fig. 16.9 Butterworth Bandpass Bandpass Analog Filter ( a) Amplitude Response and (b) Phase Response

834


16.8.4 Bandstop Filter


% Program or the design o Butterworth analog Bandstop flter clc; close all;clear all; ormat long rp5input(‘enter the passband ripple...’); rs5input(‘enter the stopband ripple...’); wp5input(‘enter the passband req...’); ws5input(‘enter the stopband req...’); s5input(‘enter the sampling req...’); w152*wp/s;w252*ws/s; [n]5buttord(w1,w2,rp,rs,’s’); wn5[w1 w2]; [b,a]5butter(n,wn,’stop’,’s’);

50:.01:pi; w [h,om]5reqs(b,a,w);



0.28


28


1000


1400


5000

The amplitude and phase responses o Butterworth bandstop analog flter are shown in Fig. 16.10.

835

MATLAB Programs

50 0

− 50 B d ni

−100 in a G

−150 −200

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1


(a) 4 2 s n ai

0 d ar in e

−2 s a h P

−4

0

0.1

0.2

0.3

0.4

0.5 (b)

0.6

0.7

0.8


Fig. 16.10 Butterworth Bandstop Bandstop Analog Filter ( a) Amplitude Response and (b) Phase Response

16.9

CHEBYSHEV TYPE-1 ANALOG FILTERS


Algorithm 1. 2. 3. 4. 5. 6.


% Program or the design o Chebyshev Type-1 low-pass flter clc; close all;clear ormat long rp5input(‘enter rs5input(‘enter wp5input(‘enter ws5input(‘enter s5input(‘enter



ripple...’); ripple...’); req...’); req...’); req...’);

836


w152*wp/s;w252*ws/s; [n,wn]5cheb1ord(w1,w2,rp,rs,’s’); [b,a]5cheby1(n,rp,wn,’s’); 50:.01:pi; w [h,om]5reqs(b,a,w); 520*log10(abs(h)); m an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel(‘Gain in dB --.’);xlabel(‘(a) Normalised requency --.’); subplot(2,1,2);plot(om/pi,an); xlabel(‘(b) Normalised requency --.’); ylabel(‘Phase in radians --.’);


the the the the the


ripple... ripple... req... req... req...

0.23 47 1300 1550 7800

The amplitude and phase responses o Chebyshev type - 1 low-pass analog flter are shown in Fig. 16.11.

0

− 20 − 40 B d ni

− 60 in a G

− 80 −100 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1


(a) 4 2 s n ia

0 d ar ni

−2 e s a h P

−4 0

0.1

0.2

0.3

0.4

0.5 (b)

Fig. 16.11

0.6

0.7

0.8


Chebyshev Type-I Low-pass Analog Filter ( a a) Amplitude Response Response and (b) Phase Response

MATLAB Programs

837

16.9.2 High-pass Filter Algorithm 1. 2. 3. 4. 5. 6.


%Program or the design o Chebyshev Type-1 high-pass flter clc; close all;clear all; ormat long rp5input(‘enter the passband ripple...’); rs5input(‘enter the stopband ripple...’); wp5input(‘enter the passband req...’); ws5input(‘enter the stopband req...’); s5input(‘enter the sampling req...’); w152*wp/s;w252*ws/s; [n,wn]5cheb1ord(w1,w2,rp,rs,’s’); [b,a]5cheby1(n,rp,wn,’high’,’s’); 50:.01:pi; w [h,om]5reqs(b,a,w); 520*log10(abs(h)); m an5angle(h); subplot(2,1,1);plot(om/pi,m); 0

− 50 B d

− 100 ni ni a G

− 150 − 200 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1


(a) 4 2 s n ia

0 d ar ni

−2 e s a h P

−4 0

0.1

0.2

0.3

0.4

0.5 (b)

Fig. 16.12

0.6

0.7

0.8


Chebyshev Type - 1 High-pass Analog Filter ( a a) Amplitude Response and (b) Phase Response

838


ylabel(‘Gain in dB --.’);xlabel(‘(a) Normalised requency --.’); subplot(2,1,2);plot(om/pi,an); xlabel(‘(b) Normalised requency --.’); ylabel(‘Phase in radians --.’);


the the the the the



0.29 29 900 1300 7500

The amplitude and phase responses o Chebyshev type - 1 high-pass analog flter are shown in Fig. 16.12.

16.9.3 Bandpass Filter Algorithm 1. 2. 3. 4. 5. 6.

Get the passband and stopband ripples Get the passband and stopband edge requencies Get the sampling requency Calculate the order o the flter using Eq. 8.57 Find the flter coefcients Draw the the magnitude magnitude and phase responses. responses.

% Program or the design o Chebyshev Type-1 Bandpass flter clc; close all;clear all; ormat long rp5input(‘enter the passband ripple...’); rs5input(‘enter the stopband ripple...’); wp5input(‘enter the passband req...’); ws5input(‘enter the stopband req...’); s5input(‘enter the sampling req...’); w152*wp/s;w252*ws/s; [n]5cheb1ord(w1,w2,rp,rs,’s’); wn5[w1 w2]; [b,a]5cheby1(n,rp,wn,’bandpass’,’s’); 50:.01:pi; w [h,om]5reqs(b,a,w); 520*log10(abs(h)); m an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel(‘Gain in dB --.’);xlabel(‘(a) Normalised requency --.’); subplot(2,1,2);plot(om/pi,an); xlabel(‘(b) Normalised requency --.’); ylabel(‘Phase in radians --.’);

As an example, enter the passband ripple... enter the stopband ripple... enter the passband req...

0.3 40 1400

MATLAB Programs

enter the stopband req... enter the sampling req...

839

2000 5000

The amplitude and phase responses o Chebyshev type - 1 bandpass analog flter are shown in Fig. 16.13. 0

− 100 B d n i n i a G

− 200 − 300 − 400 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1


a 3 2 1 s n a i d a r n i e s a h P

0

−1 −2 −3 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9


(b)

Fig. 16.13 Chebyshev Type-1 Bandpass Analog Filter ( a) Amplitude Response and (b) Phase Response


Algorithm 1. 2. 3. 4. 5. 6.

Get the passband and stopband ripples Get the passband and stopband edge requency Get the sampling requency Calculate the order order o the flter flter using using Eq. 8.57 Find the flter coefcients Draw the magnitude and phase responses.

% Program or the design o Chebyshev Type-1 Bandstop flter clc; close all;clear ormat long rp5input(‘enter rs5input(‘enter wp5input(‘enter ws5input(‘enter

all; the the the the

passband stopband passband stopband

ripple...’); ripple...’); req...’); req...’);

1

840


s5input(‘enter the sampling req...’); w152*wp/s;w252*ws/s; [n]5cheb1ord(w1,w2,rp,rs,’s’); wn5[w1 w2]; [b,a]5cheby1(n,rp,wn,’stop’,’s’); 50:.01:pi; w [h,om]5reqs(b,a,w); 520*log10(abs(h)); m an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel(‘Gain in dB --.’);xlabel(‘(a) Normalised requency --.’); subplot(2,1,2);plot(om/pi,an); xlabel(‘(b) Normalised requency --.’); ylabel(‘Phase in radians --.’);


the the the the the



0.15 30 2000 2400 7000

The amplitude and phase responses o Chebyshev type - 1 bandstop analog flter are shown in Fig. 16.14. 0

B

−50 −100 in

d

−150 G

a

ni

−200 −250 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1


(a) 4 2 s ai

n

0 ar

d

e

ni

−2 a

s P

h

−4 0

0.1

0.2

0.3

0.4

0.5 (b)

0.6

0.7

0.8


Fig. 16.14 Chebyshev Type - 1 Bandstop Analog Filter ( a) Amplitude Response and (b) Phase Response

MATLAB Programs

16.10

841

CHEBYSHEV TYPE-2 ANALOG FILTERS


Algorithm 1. 2. 3. 4. 5. 6.


% Program or the design o Chebyshev Type-2 low pass analog flter clc; close all;clear all; ormat long rp5input(‘enter the passband ripple...’); rs5input(‘enter the stopband ripple...’); wp5input(‘enter the passband req...’); ws5input(‘enter the stopband req...’); s5input(‘enter the sampling req...’); w152*wp/s;w252*ws/s; [n,wn]5cheb2ord(w1,w2,rp,rs,’s’); [b,a]5cheby2(n,rs,wn,’s’);

50:.01:pi; w [h,om]5reqs(b,a,w);



0.4


50


2000


2400


10000

The amplitude and phase responses o Chebyshev type - 2 low-pass analog flter are shown in Fig. 16.15.

842


0

− 20 − 40 B d ni

− 60 in a G

− 80 − 100 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1


(a) 4 2 s n ai

0 d ar ni

−2 e s a h P

−4 0

0.1

0.2

0.3

0.4

0.5 (b)

Fig. 16.15

0.6

0.7

0.8


Chebyshev Type - 2 Low-pass Analog Filter ( a a) Amplitude Response and (b) Phase Response


Algorithm 1. 2. 3. 4. 5. 6.


% Program or the design o Chebyshev Type-2 Type-2 High pass analog flter clc; close all;clear all; ormat long rp5input(‘enter the passband ripple...’); rs5input(‘enter the stopband ripple...’); wp5input(‘enter the passband req...’); ws5input(‘enter the stopband req...’); s5input(‘enter the sampling req...’); w152*wp/s;w252*ws/s; [n,wn]5cheb2ord(w1,w2,rp,rs,’s’); [b,a]5cheby2(n,rs,wn,’high’,’s’); 50:.01:pi; w

MATLAB Programs

843

[h,om]5reqs(b,a,w); 520*log10(abs(h)); m an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel(‘Gain in dB --.’);xlabel(‘(a) Normalised requency --.’); subplot(2,1,2);plot(om/pi,an); xlabel(‘(b) Normalised requency --.’); ylabel(‘Phase in radians --.’);


the the the the the



0.34 34 1400 1600 10000

The amplitude and phase responses o Chebyshev type - 2 high-pass analog flter are shown in Fig. 16.16. 0

− 20 B

− 40 d

ni in

− 60 a G

− 80 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1


(a) 4 2 0 s n ai d ar

−2 ni e s

−4 a h P

0

0.1

0.2

0.3

0.4

0.5 (b)

0.6

0.7

0.8


Fig. 16.16 Chebyshev Type - 2 High-pass Analog Filter ( a a) Amplitude Response and (b) Phase Response


Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge requencies

844

3. 4. 5. 6.


Get the sampling requency Calculate the order order o the flter flter using using Eq. 8.67 Find the flter coefcients Draw the magnitude and phase responses.

% Program or the design o Chebyshev Type-2 Bandpass analog flter clc; close all;clear all; ormat long rp5input(‘enter the passband ripple...’); rs5input(‘enter the stopband ripple...’); wp5input(‘enter the passband req...’); ws5input(‘enter the stopband req...’); s5input(‘enter the sampling req...’); w152*wp/s;w252*ws/s; [n]5cheb2ord(w1,w2,rp,rs,’s’); wn5[w1 w2]; [b,a]5cheby2(n,rs,wn,’bandpass’,’s’); 50:.01:pi; w [h,om]5reqs(b,a,w); 520*log10(abs(h)); m an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel(‘Gain in dB --.’);xlabel(‘(a) Normalised requency --.’); subplot(2,1,2);plot(om/pi,an); xlabel(‘(b) Normalised requency --.’); ylabel(‘Phase in radians --.’);

As an example, enter the passband ripple... enter enter enter enter

the the the the

stopband passband stopband sampling

ripple... req... req... req...

0.37 37 3000 4000 9000

The amplitude and phase responses o Chebyshev type - 2 bandpass analog flter are shown in Fig. 16.17. 20 0

− 20 − 40 ni

d

B

− 60 G

a

ni

− 80 − 100 0

0.1

0.2

0.3

0.4

0.5 (a)

Fig. 16.17

0.6

0.7

0.8


(Contd.)

0.9

1

MATLAB Programs

845

4 2 s n ia

0 d a r ni e

2 − s a h P

4 − 0

0.1

0.2

0.3

0.4

0.5 (b)

Fig. 16.17

16.10.4

0.6

0.7

0.8

0.9

1


Chebyshev Type - 2 Bandstop Analog Filter ( a a) Amplitude Response and (b) Phase Response

Bandstop Filter

Algorithm 1. 2. 3. 4. 5. 6.


% Program or the design o Chebyshev Type-2 Bandstop analog flter clc; close all;clear all; ormat long rp5input(‘enter the passband ripple...’); rs5input(‘enter the stopband ripple...’); wp5input(‘enter the passband req...’); ws5input(‘enter the stopband req...’); s5input(‘enter the sampling req...’); w152*wp/s;w252*ws/s; [n]5cheb2ord(w1,w2,rp,rs,’s’); wn5[w1 w2]; [b,a]5cheby2(n,rs,wn,’stop’,’s’); 50:.01:pi; w [h,om]5reqs(b,a,w); 520*log10(abs(h)); m an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel(‘Gain in dB --.’);xlabel(‘(a) Normalised requency --.’); subplot(2,1,2);plot(om/pi,an); xlabel(‘(b) Normalised requency --.’); ylabel(‘Phase in radians --.’);

846



the the the the the



0.25 30 1300 2000 8000

The amplitude and phase responses o Chebyshev type - 2 bandstop analog flter are shown in Fig. 16.18. 40 20 0 B d ni

− 20 a

− 40

in G

− 60 − 80 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1


(a) 4 2 s n ai

0 d ar in

−2 e s a h P

−4 0

0.1

0.2

0.3

0.4

0.5

0.6

(b)

0.7

0.8


Fig. 16.18 Chebyshev Type - 2 Bandstop Analog Filter ( a) Amplitude Response and (b) Phase Response

16.11

BUTTERWORTH DIGITAL IIR FILTERS


Algorithm 1. 2. 3. 4. 5. 6.


MATLAB Programs

847

% Program or the design o Butterworth low pass digital flter clc; close all;clear all; ormat long rp5input(‘enter the passband ripple’); rs5input(‘enter the stopband ripple’); wp5input(‘enter the passband req’); ws5input(‘enter the stopband req’); s5input(‘enter the sampling req’); w152*wp/s;w252*ws/s; [n,wn]5buttord(w1,w2,rp,rs); [b,a]5butter(n,wn);

50:.01:pi; w [h,om]5reqz(b,a,w);



0.5


50


1200


2400


10000

The amplitude and phase responses o Butterworth low-pass digital flter are shown in Fig. 16.19. 100 0

− 100 ni

d

B

− 200 ni G

a

− 300 − 400 0

0.1

0.2

0.3

0.4

0.5

0.6

0.8


(a)

Fig. 16.19

0.7

(Contd.)

0.9

1

848


4 2 s n ai

0 d a r in

−2 e s a h P

−4 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9


(b)

Fig. 16.19 Butterworth Low-pass Low-pass Digital Filter ( a a) Amplitude Response and (b) Phase Response


Algorithm 1. 2. 3. 4. 5. 6.


% Program or the design o Butterworth highpass digital flter clc; close all;clear all; ormat long rp5input(‘enter the passband ripple’); rs5input(‘enter the stopband ripple’); wp5input(‘enter the passband req’); ws5input(‘enter the stopband req’); s5input(‘enter the sampling req’); w152*wp/s;w252*ws/s; [n,wn]5buttord(w1,w2,rp,rs); [b,a]5butter(n,wn,’high’); 50:.01:pi; w [h,om]5reqz(b,a,w); 520*log10(abs(h)); m an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel(‘Gain in dB --.’);xlabel(‘(a) Normalised requency --.’); subplot(2,1,2);plot(om/pi,an); xlabel(‘(b) Normalised requency --.’); ylabel(‘Phase in radians --.’);

As an example, enter the passband ripple enter the stopband ripple enter the passband req

0.5 50 1200

1

MATLAB Programs

enter the stopband req enter the sampling req

849

2400 10000

The amplitude and phase responses o Butterworth high-pass digital flter are shown in Fig. 16.20. 50 0

− 50 B

− 100 d ni

− 150 in a G

− 200 − 250 − 300 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1


(a) 4 2 s n ia

0 d ar ni

−2 e s a h P

−4 0

0.1

0.2

0.3

0.4

0.5

0.6

(b)

0.7

0.8


Fig. 16.20 Butterworth High-pass High-pass Digital Filter ( a) Amplitude Response and (b) Phase Response

16.11.3 Band-pass Filter

Algorithm 1. 2. 3. 4. 5. 6.


% Program or the design o Butterworth Bandpass digital flter clc; close all;clear ormat long rp5input(‘enter rs5input(‘enter wp5input(‘enter

all; the passband ripple’); the stopband ripple’); the passband req’);

850


ws5input(‘enter the stopband req’); s5input(‘enter the sampling req’); w152*wp/s;w252*ws/s; [n]5buttord(w1,w2,rp,rs); wn5[w1 w2]; [b,a]5butter(n,wn,’bandpass’); 50:.01:pi; w [h,om]5reqz(b,a,w); 520*log10(abs(h)); m an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel(‘Gain in dB --.’);xlabel(‘(a) Normalised requency --.’); subplot(2,1,2);plot(om/pi,an); xlabel(‘(b) Normalised requency --.’); ylabel(‘Phase in radians --.’);


the the the the the



0.3 40 1500 2000 9000

The amplitude and phase responses o Butterworth band-pass digital flter are shown in Fig. 16.21. 0

− 100 − 200 − 300 B d

− 400 ni ni

− 500 a G

− 600 − 700

0

0.1

0.2

0.3

0.4

0.5 (a)

0.6

0.7

0.8

0.9

1


4 2 s n ai

0 d ar ni e

−2 s a h P

−4

0

0.1

0.2

0.3

0.4 (b)

0.5

0.6

0.7

0.8

0.9


Fig. 16.21 Butterworth Bandstop Bandstop Digital Filter ( a a) Amplitude Response and (b) Phase Response

1

MATLAB Programs

851



% Program or the design o Butterworth Band stop digital flter clc; close all;clear all; ormat long rp5input(‘enter the passband ripple’); rs5input(‘enter the stopband ripple’); wp5input(‘enter the passband req’); ws5input(‘enter the stopband req’); s5input(‘enter the sampling req’); w152*wp/s;w252*ws/s; [n]5buttord(w1,w2,rp,rs); wn5[w1 w2]; [b,a]5butter(n,wn,’stop’);

50:.01:pi; w [h,om]5reqz(b,a,w);



0.4


46


1100


2200


6000

The amplitude and phase responses o the Butterworth bandstop digital flter are shown in Fig. 16.22.

852


100 0

− 100 B d

− 200 ni in a G

− 300 − 400

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1


(a) 4 2 s n ai

0 d ar in

−2 e s a h P

−4 0

0.1

0.2

0.3

0.4

0.5 (b)

0.6

0.7

0.8


Fig. 16.22 Butterworth Bandstop Bandstop Digital Filter ( a) Amplitude Response and (b) Phase Response

16.12

CHEBYSHEV TYPE-1 DIGITAL FILTERS


Algorithm 1. 2. 3. 4. 5. 6.


% Program or the design o Chebyshev Type-1 lowpass digital flter clc; close all;clear ormat long rp5input(‘enter rs5input(‘enter wp5input(‘enter ws5input(‘enter s5input(‘enter




MATLAB Programs

853

w152*wp/s;w252*ws/s; [n,wn]5cheb1ord(w1,w2,rp,rs); [b,a]5cheby1(n,rp,wn); 50:.01:pi; w [h,om]5reqz(b,a,w); 520*log10(abs(h)); m an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel(‘Gain in dB --.’);xlabel(‘(a) Normalised requency --.’); subplot(2,1,2);plot(om/pi,an); xlabel(‘(b) Normalised requency --.’); ylabel(‘Phase in radians --.’);


the the the the the



0.2 45 1300 1500 10000

The amplitude and phase responses o Chebyshev type - 1 low-pass digital flter are shown in Fig. 16.23. 0

− 100 − 200 B d

− 300 ni ni a G

− 400 − 500

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1


(a) 4 2 s n ai

0 d ar ni

−2 e s a h P

−4 0

0.1

0.2

0.3

0.4

0.5 (b)

Fig. 16.23

0.6

0.7

0.8


Chebyshev Type - 1 Low-pass Digital Filter ( a a) Amplitude Response and (b) Phase Response

854



Algorithm 1. 2. 3. 4. 5. 6.


% Program or the design o Chebyshev Type-1 highpass digital ilter il ter clc; close all;clear all; ormat long rp5input(‘enter the passband ripple...’); rs5input(‘enter the stopband ripple...’); wp5input(‘enter the passband req...’); ws5input(‘enter the stopband req...’); s5input(‘enter the sampling req...’); w152*wp/s;w252*ws/s; [n,wn]5cheb1ord(w1,w2,rp,rs); [b,a]5cheby1(n,rp,wn,’high’);

50:.01/pi:pi; w [h,om]5reqz(b,a,w);



0.3


60


1500


2000


9000

The amplitude and phase responses o Chebyshev type - 1 high-pass digital flter are shown in Fig. 16.24.

MATLAB Programs

855

0

− 50 −100 −150 B d

− 200 ni ni

− 250 a G

− 300 − 350

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1


(a) 4 2 s n ia

0 d ar ni e

−2 s a h P

−4 0

0.1

0.2

0.3

0.4

0.5 (b)

0.6

0.7

0.8


Fig. 16.24 Chebyshev Type - 1 High-pass Digital Filter ( a) Amplitude Response and (b) Phase Response


Algorithm 1. 2. 3. 4. 5. 6.


% Program or the design o Chebyshev Type-1 Type-1 Bandpass digital flter clc; close all;clear all; ormat long rp5input(‘enter the passband rs5input(‘enter the stopband wp5input(‘enter the passband ws5input(‘enter the stopband s5input(‘enter the sampling w152*wp/s;w252*ws/s;


856


[n]5cheb1ord(w1,w2,rp,rs); wn5[w1 w2]; [b,a]5cheby1(n,rp,wn,’bandpass’);

50:.01:pi; w [h,om]5reqz(b,a,w);



the the the the the



0.4 35 2000 2500 10000

The amplitude and phase responses o Chebyshev type - 1 bandpass digital flter are shown in Fig. 16.25. 0

− 100 B d ni

ni

− 200 − 300

G

a

− 400 − 500

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1


(a) 4 2 s ia

n

0 ar

d

e

ni

−2 a

s P

h

−4 0

0.1

0.2

0.3

0.4

0.5 (b)

Fig. 16.25

0.6

0.7

0.8


Chebyshev Type - 1 Bandpass Digital Filter ( a) Amplitude Response and (b) Phase Response

MATLAB Programs

857


Algorithm 1. 2. 3. 4. 5. 6.


% Program or the design o Chebyshev Type-1 Bandstop digital flter clc; close all;clear all; ormat long rp5input(‘enter the passband ripple...’); rs5input(‘enter the stopband ripple...’); wp5input(‘enter the passband req...’); ws5input(‘enter the stopband req...’); s5input(‘enter the sampling req...’); w152*wp/s;w252*ws/s; [n]5cheb1ord(w1,w2,rp,rs); wn5[w1 w2]; [b,a]5cheby1(n,rp,wn,’stop’);

50:.1/pi:pi; w [h,om]5reqz(b,a,w);



0.25


40


2500


2750


7000

The amplitude and phase responses o Chebyshev type - 1 bandstop digital flter are shown in Fig. 16.26.

858


0

− 50 B

− 100 d ni ni a G

− 150 − 200 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1


(a) 4 3 2 s

1 a

d

0 ni

−1 a

s

−2 P

−3

ai

n r e h

0

0.1

0.2

0.3

0.4

0.5 (b)

0.6

0.7

0.8


Fig. 16.26 Chebyshev Type - 1 Bandstop Digital Filter ( a) Amplitude Response and (b) Phase Response

16.13

CHEBYSHEV TYPE-2 DIGITAL FILTERS

16.13.1 Low-pass Filter Algorithm 1. 2. 3. 4. 5. 6.


% Program or the design o Chebyshev Type-2 lowpass digital flter clc; close all;clear all; ormat long rp5input(‘enter the passband ripple...’); rs5input(‘enter the stopband ripple...’); wp5input(‘enter the passband req...’); ws5input(‘enter the stopband req...’); s5input(‘enter the sampling req...’); w152*wp/s;w252*ws/s; [n,wn]5cheb2ord(w1,w2,rp,rs); [b,a]5cheby2(n,rs,wn);

MATLAB Programs

859

50:.01:pi; w [h,om]5reqz(b,a,w); 520*log10(abs(h)); m an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel(‘Gain in dB --.’);xlabel(‘(a) Normalised requency --.’); subplot(2,1,2);plot(om/pi,an); xlabel(‘(b) Normalised requency --.’); ylabel(‘Phase in radians --.’);


the the the the the



0.35 35 1500 2000 8000

The amplitude and phase responses o Chebyshev type - 2 low-pass digital flter are shown in Fig. 16.27. 20 0

−20 B d ni

−40 a

−60

in G

−80 − 100 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1


(a) 4 2 s n ia

0 d ar ni e

−2 s a h P

−4 0

0.1

0.2

0.3

0.4

0.5 (b)

Fig. 16.27

0.6

0.7

0.8


Chebyshev Type - 2 Low-pass Digital Filter ( a a) Amplitude Response and (b) Phase Response


Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge requencies

860


3. Get the sampling requency 4. Calculate the order order o the flter flter using using Eq. 8.67 5. Find the flter coefcients 6. Draw the magnitude and phase responses.

% Program or the design o Chebyshev Type-2 high pass digital flter clc; close all;clear all; ormat long rp5input(‘enter the passband ripple...’); rs5input(‘enter the stopband ripple...’); wp5input(‘enter the passband req...’); ws5input(‘enter the stopband req...’); s5input(‘enter the sampling req...’); w152*wp/s;w252*ws/s; [n,wn]5cheb2ord(w1,w2,rp,rs); [b,a]5cheby2(n,rs,wn,’high’); 50:.01/pi:pi; w [h,om]5reqz(b,a,w); 520*log10(abs(h)); m an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel(‘Gain in dB --.’);xlabel(‘(a) Normalised requency --.’); subplot(2,1,2);plot(om/pi,an); xlabel(‘(b) Normalised requency --.’); ylabel(‘Phase in radians --.’);


the the the the the



0.25 40 1400 1800 7000

The amplitude and phase responses o Chebyshev type - 2 high-pass digital flter are shown in Fig. 16.28. 0

− 20 − 40 B ni

d

− 60 − 80

G

− 100

a

ni

− 120 0

0.1

0.2

0.3

0.4

0.5

0.6

0.8


(a)

Fig. 16.28

0.7

(Contd.)

0.9

1

MATLAB Programs

861

4 2 s n ai

0 d ar in

−2 e s a h P

−4 0

0.1

0.2

0.3

0.4

0.5 (b)

Fig. 16.28

0.6

0.7

0.8

0.9

1


Chebyshev Type - 2 High-pass Digital Filter ( a a) Amplitude Response and (b) Phase Response


Algorithm 1. 2. 3. 4. 5. 6.

Get the passband and stopband ripples Get the passband and stopband edge requency Get the sampling requency Calculate the order order o the flter flter using using Eq. 8.67 Find the flter coefcients Draw the magnitude and phase responses.

% Program or the design o Chebyshev Type-2 Type-2 Bandpass digital flter clc; close all;clear all; ormat long rp5input(‘enter the passband ripple...’); rs5input(‘enter the stopband ripple...’); wp5input(‘enter the passband req...’); ws5input(‘enter the stopband req...’); s5input(‘enter the sampling req...’); w152*wp/s;w252*ws/s; [n]5cheb2ord(w1,w2,rp,rs); wn5[w1 w2]; [b,a]5cheby2(n,rs,wn,’bandpass’); 50:.01/pi:pi; w [h,om]5reqz(b,a,w); 520*log10(abs(h)); m an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel(‘Gain in dB --.’);xlabel(‘(a) Normalised requency --.’); subplot(2,1,2);plot(om/pi,an); xlabel(‘(b) Normalised requency --.’); ylabel(‘Phase in radians --.’);

862



the the the the the



0.4 40 1400 2000 9000

The amplitude and phase responses o Chebyshev type - 2 bandpass digital flter are shown in Fig. 16.29. 100 0 B

−100 in

−200

in

d

a G

−300 −400 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1


(a) 4 2 s n ia

0 d a r ni

−2 e s a h P

−4 0

0.1

0.2

0.3

0.4

0.5 (b)

Fig. 16.29

0.6

0.7

0.8


Chebyshev Type - 2 Bandpass Digital Filter ( a a) Amplitude Response and (b) Phase Response



MATLAB Programs

863

% Program or the design o Chebyshev Type-2 Bandstop digital flter clc; close all;clear all; ormat long rp5input(‘enter the passband ripple...’); rs5input(‘enter the stopband ripple...’); wp5input(‘enter the passband req...’); ws5input(‘enter the stopband req...’); s5input(‘enter the sampling req...’); w152*wp/s;w252*ws/s; [n]5cheb2ord(w1,w2,rp,rs); wn5[w1 w2]; [b,a]5cheby2(n,rs,wn,’stop’); 50:.1/pi:pi; w [h,om]5reqz(b,a,w); 520*log10(abs(h)); m an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel(‘Gain in dB --.’);xlabel(‘(a) Normalised requency --.’); subplot(2,1,2);plot(om/pi,an); xlabel(‘(b) Normalised requency --.’); ylabel(‘Phase in radians --.’);


the the the the the



0.3 46 1400 2000 8000

The amplitude and phase responses o Chebyshev type - 2 bandstop digital flter are shown in Fig. 16.30.

20 0

− 20 ni

d

B

− 40 G

a

ni

− 60 − 80 0

0.1

0.2

0.3

0.4

0.5 (a)

Fig. 16.30

0.6

0.7

0.8


(Contd.)

0.9

1

864


3 2 1 s n

0 ia d

-1 ar ni

-2 e s a

-3 h P

-4 0

0.1

0.2

0.3

0.4

0.5

0.6

(b)

Fig. 16.30

0.7

0.8

0.9

1


Chebyshev Type - 2 Bandstop Digital Filter ( a a) Amplitude Response and (b) Phase Response

FIR FILTER DESIGN USING WINDOW TECHNIQUES

16.14

In the design o FIR flters using any window technique, the order can be calculated using the ormula given by N

20 log( d pd s ) 13

− =

−

14.6 ( f s

−

f p ) / F s

where d p is the passband ripple, d s is the stopband ripple, f ripple, f p is the passband requency, f s is the stopband requency and F and F s is the sampling requency.

16.14.1

Rectangular Rectangul ar Window

Algorithm 1. 2. 3. 4. 5. 6.

Get the passband and stopband ripples Get the passband and stopband edge requencies Get the sampling requency Calculate the order o the flter Find the window coefcients using Eq. 7.37 Draw the magnitude and phase responses.

% Program or the design o FIR Low pass, High pass, Band pass and Bandstop flters using rectangular window windo w clc;clear all;close all; rp5input(‘enter the passband ripple’); rs5input(‘enter the stopband ripple’); p5input(‘enter the passband req’); s5input(‘enter the stopband req’); 5input(‘enter the sampling req’); wp52*p/;ws52*s/; 52 5220*log10(sqrt(rp*rs))213; num

MATLAB Programs

865

514.6*(s2p)/; dem n5ceil(num/dem); n15n11; i (rem(n,2)˜50) n15n; n5n21; end y5boxcar(n1); %

low-pass filter

b5r1(n,wp,y); [h,o]5reqz(b,1,256); 520*log10(abs(h)); m subplot(2,2,1);plot(o/pi,m);ylabel(‘Gain in dB --.’); xlabel(‘(a) Normalised requency --.’); %

high-pass filter

b5r1(n,wp,’high’,y); [h,o]5reqz(b,1,256); 520*log10(abs(h)); m subplot(2,2,2);plot(o/pi,m);ylabel(‘Gain in dB --.’); xlabel(‘(b) Normalised requency --.’); %

band pass filter

wn5[wp ws]; b5r1(n,wn,y); [h,o]5reqz(b,1,256); 520*log10(abs(h)); m subplot(2,2,3);plot(o/pi,m);ylabel(‘Gain in dB -->’); xlabel(‘(c) Normalised requency -->’); %

band stop filter

b5r1(n,wn,’stop’,y); [h,o]5reqz(b,1,256); 520*log10(abs(h)); m subplot(2,2,4);plot(o/pi,m);ylabel(‘Gain in dB -->’); xlabel(‘(d) Normalised requency -->’);


the the the the the



0.05 0.04 1500 2000 9000

The gain responses o low-pass, high-pass, bandpass and bandstop flters using rectangular window are shown in Fig. 16.31.

866

B d n i n i a G


20

20

0

0

− 20

B d n i

− 40

n i a G

− 60 − 80

0

0.2

0.4

0.6

0.8

− 20 − 40 − 60 − 80

1

0


0.2

0.8

1

n i n i a G

0.2 0.4 0.6 0.8 Normalised frequency

1

(b)

20

5

0

0

− 20

B d n i

− 40

n i a G

− 60 − 80

0


1

−5 −10 −15 −20

0

(c)

Fig. 16.31

16.14.2

0.6


(a)

B d

0.4

(d)

Filters Using Rectangular Window ( a) Low-pass (b) High-pass ( c c) Bandpass and ( d) Bandstop

Bartlett Window

Algorithm 1. 2. 3. 4. 5. 6.

Get the passband and stopband ripples Get the passband and stopband edge requencies Get the sampling requency Calculate the order o the flter Find the flter coefcients Draw the magnitude and phase responses.

% Program or the design o FIR Low pass, High pass, Band pass and Bandstop flters using Bartlett window clc;clear all;close all; rp5input(‘enter the passband ripple’); rs5input(‘enter the stopband ripple’); p5input(‘enter the passband req’); s5input(‘enter the stopband req’); 5input(‘enter the sampling req’);

MATLAB Programs

867

wp52*p/;ws52*s/; 52 5220*log10(sqrt(rp*rs))213; num 514.6*(s2p)/; dem n5ceil(num/dem); n15n11; i (rem(n,2)˜50) n15n; n5n21; end y5bartlett(n1); %

low-pass filter


high-pass filter


band pass filter

wn5[wp ws]; b5r1(n,wn,y); [h,o]5reqz(b,1,256); 520*log10(abs(h)); m subplot(2,2,3);plot(o/pi,m);ylabel(‘Gain in dB --.’); xlabel(‘(c) Normalised requency --.’); %

band stop filter

b5r1(n,wn,’stop’,y); [h,o]5reqz(b,1,256); 520*log10(abs(h)); m subplot(2,2,4);plot(o/pi,m);ylabel(‘Gain in dB --.’); xlabel(‘(d) Normalised requency --.’);


the the the the the



0.04 0.02 1500 2000 8000

The gain responses o low-pass, high-pass, bandpass and bandstop flters using Bartlett window are shown in Fig. 16.32.

868


0

5

−5

0

−5

− 10 B d

− 15 in

− 20 G

− 25

ni a

B d

− 10 in

− 15 G

− 20

ni a

− 30

− 25

− 35

0

0.2

0.4

0.6

0.8

− 30

1

0


0.2

0.8

1


1

(b)

0

2

−0

− 10 B

B d

−2 in

−4

d in

− 20 in a G

a G

− 30 − 40

0

0.6


(a)

in

0.4

−6 0.2 0.4 0.6 0.8 Normalised frequency

1

−8

0

(c)

Fig. 16.32

(d)

Filters using Bartlett Window ( a) Low-pass (b) High-pass ( c c) Bandpass and ( d) Bandstop

16.14.3 Blackman window

Algorithm 1. 2. 3. 4. 5. 6.


% Program or the design o FIR Low pass, High pass, Band pass and Band stop digital flters using Blackman window wind ow clc;clear all;close all; rp5input(‘enter the passband ripple’); rs5input(‘enter the stopband ripple’); p5input(‘enter the passband req’); s5input(‘enter the stopband req’); 5input(‘enter the sampling req’);

MATLAB Programs

869

wp52*p/;ws52*s/; 52 5220*log10(sqrt(rp*rs))213; num 514.6*(s2p)/; dem n5ceil(num/dem); n15n11; i (rem(n,2)˜50) n15n; n5n21; end y5blackman(n1); %

low-pass filter


high-pass filter


band pass filter


band stop filter

b5r1(n,wn,’stop’,y); [h,o]5reqz(b,1,256); 520*log10(abs(h)); m subplot(2,2,4);plot(o/pi,m);;ylabel(‘Gain in dB --.’); xlabel(‘(d) Normalised requency --.’);


the the the the the



0.03 0.01 2000 2500 7000

The gain responses o low-pass, high-pass, bandpass and bandstop flters using Blackman window are shown in Fig. 16.33.

870


20

50

0 0

− 20 B

B

− 40 d ni

d in

− 50

− 60 ni a

in a

− 80 G

G

− 100 − 120

0

0.2

0.4

0.6

0.8

1

− 100 − 150

0


0.2

0.8

1


1

(b)

0

2

− 20

0

− 40 d ni

− 60 a

− 80

in G

B

0

d

−2 in

−4

ni a G

−6

−100 −120

0.6


(a)

B

0.4


1

−8

0

(c)

Fig. 16.33

(d)

Filters using Blackman Window ( a) Low-pass (b) High-pass ( c c) Bandpass and ( d) Bandstop

16.14.4 Chebyshev Window

Algorithm 1. 2. 3. 4. 5. 6.

Get the passband and stopband ripples Get the passband and stopband edge requencies Get the sampling requency Calculate the order o the flter Find the flter coefcients Draw the magnitude and phase responses.

% Program or the design o FIR Lowpass, High pass, Band pass and Bandstop flters using Chebyshev window clc;clear all;close all; rp5input(‘enter the passband ripple’); rs5input(‘enter the stopband ripple’); p5input(‘enter the passband req’); s5input(‘enter the stopband req’); 5input(‘enter the sampling req’);

MATLAB Programs

871

r5input(‘enter the ripple value(in dBs)’); wp52*p/;ws52*s/; 52 5220*log10(sqrt(rp*rs))213; num 514.6*(s2p)/; dem n5ceil(num/dem); i(rem(n,2)˜50) n5n11; end y5chebwin(n,r); %

low-pass filter

b5r1(n-1,wp,y); [h,o]5reqz(b,1,256); 520*log10(abs(h)); m subplot(2,2,1);plot(o/pi,m);ylabel(‘Gain in dB --.’); xlabel(‘(a) Normalised requency --.’); %

high-pass filter

b5r1(n21,wp,’high’,y); [h,o]5reqz(b,1,256); 520*log10(abs(h)); m subplot(2,2,2);plot(o/pi,m);ylabel(‘Gain in dB --.’); xlabel(‘(b) Normalised requency --.’); %

band-pass filter

wn5[wp ws]; b5r1(n21,wn,y); [h,o]5reqz(b,1,256); 520*log10(abs(h)); m subplot(2,2,3);plot(o/pi,m);ylabel(‘Gain in dB --.’); xlabel(‘(c) Normalised requency --.’); %

band-stop filter

b5r1(n21,wn,’stop’,y); [h,o]5reqz(b,1,256); 520*log10(abs(h)); m subplot(2,2,4);plot(o/pi,m);ylabel(‘Gain in dB --.’); xlabel(‘(d) Normalised requency --.’);

As an example, enter enter enter enter enter enter

the the the the the the

passband ripple 0.03 stopband ripple 0.02 passband req 1800 stopband req 2400 sampling req 10000 ripple value(in dBs)40

The gain responses o low-pass, high-pass, bandpass and bandstop flters using Chebyshev window are shown in Fig. 16.34.

872


20

20

0

0

− 20

− 20 B d in

− 40 a

−60

ni G

B d in

− 40 a

− 60

ni G

− 80

− 80 −100

0

0.2

0.4

0.6

0.8

1

− 100

0


0.2

0.8

1


1

(b)

0

2

− 20

0

−2

− 40 d

− 60 a

− 80

in

ni G

B

0

d

−4 in

−6 G

−8

in a

−100 −120


1

−12

0

(c)

Fig. 16.34

16.14.5

0.6


(a)

B

0.4

(d)

Filters using Chebyshev Window ( a) Low-pass (b) High-pass ( c c) Bandpass and ( d) Bandstop

Hamming Window

Algorithm 1. 2. 3. 4. 5. 6.


% Program or the design o FIR Low pass, High pass, Band pass and Bandstop flters using Hamming window clc;clear all;close all; rp5input(‘enter the passband ripple’); rs5input(‘enter the stopband ripple’); p5input(‘enter the passband req’); s5input(‘enter the stopband req’); 5input(‘enter the sampling req’);

MATLAB Programs

873

wp52*p/;ws52*s/; 52 5220*log10(sqrt(rp*rs))213; num 514.6*(s2p)/; dem n5ceil(num/dem); n15n11; i (rem(n,2)˜50) n15n; n5n21; end y5hamming(n1); %

low-pass filter


high-pass filter


band pass filter


band stop filter



the the the the the



0.02 0.01 1200 1700 9000

The gain responses o low-pass, high-pass, bandpass and bandstop flters using Hamming window are shown in Fig. 16.35.

874


20

20

0

0

− 20 B d

− 40 in

− 60 G

−80

ni a

d in

− 40 a

− 60

ni G

− 80

−100 −120

− 20 B

0

0.2

0.4

0.6

0.8

− 100

1

0


0.2

0.8

1


1

(b)

0

2

− 20

0

− 40 d

− 60 a

− 80

in

ni G

B d ni

0

−5 in a G

−100 −120

0.6


(a)

B

0.4


1

− 10 − 15

0

(c)

Fig. 16.35

(d)

Filters using Hamming Window ( a) Low-pass (b) High-pass ( c c) Bandpass and ( d) Bandstop

16.14.6 Hanning Window

Algorithm 1. 2. 3. 4. 5. 6.


% Program or the design o FIR Low pass, High pass, Band pass and Band stop flters using Hanning window wi ndow clc;clear all;close all; rp5input(‘enter the passband ripple’); rs5input(‘enter the stopband ripple’); p5input(‘enter the passband req’); s5input(‘enter the stopband req’); 5input(‘enter the sampling req’);

MATLAB Programs

875

wp52*p/;ws52*s/; 52 5220*log10(sqrt(rp*rs))213; num 514.6*(s2p)/; dem n5ceil(num/dem); n15n11; i (rem(n,2)˜50) n15n; n5n21; end y5hamming(n1); %

low-pass filter


high-pass filter


band pass filter


band stop filter



the the the the the



0.03 0.01 1400 2000 8000

The gain responses o low-pass, high-pass, bandpass and bandstop flters using Hanning window are shown in Fig. 16.36.

876


20

20

0

0

− 20 B d

− 40 ni

− 60 G

− 80

ni a

d

0

ni

− 40 a

− 60

ni G

− 80

− 100 − 120

− 20 B

0.2

0.4

0.6

0.8

− 100

1

0


0.2

0.8

1


1

(b)

0

2

− 20

0

−2

− 40 d in

B

−4 ni

−6 G

−8

in

d

− 60 ni a G

a

− 80 −100 −120

0


1

−12

0

(c)

Fig. 16.36

16.14.7

0.6


(a)

B

0.4

(d)

Filters using Hanning Window ( a a) Low-pass (b) High-pass ( c c) Bandpass and ( d) Bandstop

Kaiser Window

Algorithm 1. 2. 3. 4. 5. 6.

Get the passband and stopband ripples Get the passband and stopband edge requencies Get the sampling requency Calculate the order o the flter Find the window coefcients using Eqs 7.46 and 7.47 Draw the magnitude and phase responses.

% Program or the design o FIR Low pass, High pass, Band pass and Bandstop flters using Kaiser window wind ow clc;clear all;close all; rp5input(‘enter the passband ripple’); rs5input(‘enter the stopband ripple’); p5input(‘enter the passband req’); s5input(‘enter the stopband req’); 5input(‘enter the sampling req’); beta5input(‘enter the beta value’);

MATLAB Programs

877

wp52*p/;ws52*s/; 52 5220*log10(sqrt(rp*rs))213; num 514.6*(s2p)/; dem n5ceil(num/dem); n15n11; i (rem(n,2)˜50) n15n; n5n21; end y5kaiser(n1,beta); %

low-pass filter

b5r1(n,wp,y); [h,o]5reqz(b,1,256); 520*log10(abs(h)); m subplot(2,2,1);plot(o/pi,m);ylabel(‘Gain in dB -->’); xlabel(‘(a) Normalised requency -->’); %

high-pass filter

b5r1(n,wp,’high’,y); [h,o]5reqz(b,1,256); 520*log10(abs(h)); m subplot(2,2,2);plot(o/pi,m);ylabel(‘Gain in dB -->’); xlabel(‘(b) Normalised requency -->’); %

band pass filter

wn5[wp ws]; b5r1(n,wn,y); [h,o]5reqz(b,1,256); 520*log10(abs(h)); m subplot(2,2,3);plot(o/pi,m);ylabel(‘Gain in dB -->’); xlabel(‘(c) Normalised requency -->’); %

band stop filter

b5r1(n,wn,’stop’,y); [h,o]5reqz(b,1,256); 520*log10(abs(h)); m subplot(2,2,4);plot(o/pi,m);ylabel(‘Gain in dB -->’); xlabel(‘(d) Normalised requency -->’);

As an example, enter enter enter enter enter enter

the the the the the the

passband ripple stopband ripple passband req stopband req sampling req beta value

0.02 0.01 1000 1500 10000 5.8

The gain responses o low-pass, high-pass, bandpass and bandstop flters using Kaiser window are shown in Fig. 16.37.

878


20

20 0

0

− 20 B

− 20 B

− 40 d in

d in

− 60 ni a

a

−80 G

G

−100 −120

− 40 ni

0

0.2

0.4

0.6

0.8

− 60 − 80

1

0


0.2

0.4

0.6

0.8

1


1


(a)

(b)

0

5

− 20 0

B

− 40 d

− 60 a

− 80

in

in G

B d in a G

−100 −120

0


1

(c)

Fig. 16.37

16.15

−5 in

−10 −15

0

(d)

Filters using Kaiser Window Window ( a) Low-pass (b) High-pass ( c c) Bandpass and ( d) Bandstop

UPSAMPLING A SINUSOIDAL SIGNAL

% Program or upsampling a sinusoidal signal by actor L N5input(‘Input length o the sinusoidal sequence5’); L5input(‘Up Samping actor5’); 5input(‘Input signal requency5’);

% Generate the sinusoidal sequence or the specifed length N n50:N21; x5sin(2*pi**n);

% Generate the upsampled signal y5zeros (1,L*length(x)); y([1:L:length(y)])5x; %Plot the input sequence subplot (2,1,1); stem (n,x); title(‘Input Sequence’); xlabel(‘Time n’); ylabel(‘Amplitude’);

MATLAB Programs

879

%Plot the output sequence subplot (2,1,2); stem (n,y(1:length(x))); title(‘[output sequence,upsampling actor5‘,num2str(L)]); xlabel(‘Time n’); ylabel(‘Amplitude’);

16.16

UPSAMPLING AN EXPONENTIAL SEQUENCE

% Program or upsampling an exponential sequence by a actor M n5input(‘enter length o input sequence …’); l5input(‘enter up sampling actor …’);

% Generate the exponential sequence 50:n21; m a5input(‘enter the value o a …’); x5a.^m;

% Generate the upsampled signal y5zeros(1,l*length(x)); y([1:l:length(y)])5x; gure(1) stem(m,x); xlabel({‘Time n’;’(a)’}); ylabel(‘Amplitude’); gure(2) stem(m,y(1:length(x))); xlabel({‘Time n’;’(b)’}); ylabel(‘Amplitude’);

As an example, enter length o input sentence … 25 enter upsampling actor … 3 enter the value o a … 0.95

The input and output sequences o upsampling an exponential sequence a n are shown in Fig. 16.38.

Fig. 16.38

(Contd.)

880


Fig. 16.38 ( a a) Input Exponential Sequence (b) Output Sequence Upsampled Upsampled by a Factor of 3

16.17

DOWN SAMPLING A SINUSOIDAL SEQUENCE

% Program or down sampling a sinusoidal sequence by a actor M N5input(‘Input length o the sinusoidal signal5’); M5input(‘Down samping actor5’); 5input(‘Input signal requency5’);

%Generate the sinusoidal sequence n50:N21; 50:N*M21; m x5sin(2*pi**m);

%Generate the down sampled signal y5x([1:M:length(x)]); %Plot the input sequence subplot (2,1,1); stem(n,x(1:N)); title(‘Input Sequence’); xlabel(‘Time n’); ylabel(‘Amplitude’); %Plot the down sampled signal sequence subplot(2,1,2); stem(n,y); title([‘Output sequence down sampling actor’,num2str(M)]); xlabel(‘Time n’); ylabel(‘Amplitude’);

16.18

DOWN SAMPLING AN EXPONENTIAL SEQUENCE

% Program or downsampling an exponential sequence by a actor M N5input(‘enter the length o the output sequence …’); M5input(‘enter the down sampling actor …’);

MATLAB Programs

881

% Generate the exponential sequence n50:N21; 50:N*M21; m a5input(‘enter the value o a …’); x5a.^m;

% Generate the downsampled signal y5x([1:M:length(x)]); gure(1) stem(n,x(1:N)); xlabel({‘Time n’;’(a)’}); ylabel(‘Amplitude’); gure(2) stem(n,y); xlabel({‘Time n’;’(b)’}); ylabel(‘Amplitude’);

As an example, enter the length o the output sentence … 25 enter the downsampling actor … 3 enter the value o a … 0.95

The input and output sequences o downsampling an exponential sequence an are shown in Fig. 16.39.

Fig. 16.39 ( a) a) Input Exponential Exponential Sequence (b) Output Sequence Downsampled Downsampled by a Factor of 3

882


16.19

DECIMATOR

% Program or downsampling the sum o two sinusoids using MATLAB’ MA TLAB’ss inbuilt decimation unction by a actor a ctor M N5input(‘Length o the input signal5’); M5input(‘Down samping actor5’); 15input(‘Frequency o rst sinusoid5’); 25input(‘Frequency o second sinusoid5’); n50:N21;

% Generate the input sequence x52*sin(2*pi*1*n)13*sin(2*pi*2*n);

%Generate the decimated signal % FIR low pass decimation is used y5decimate(x,M,‘r’); %Plot the input sequence subplot (2,1,1); stem (n,x(1:N)); title(‘Input Sequence’); xlabel(‘Time n’); ylabel(‘Amplitude’); %Plot the output sequence subplot (2,1,2); 50:N/M21; m stem (m,y(1:N/M)); title([‘Output sequence down sampling actor’,num2str(M)]); xlabel(‘Time n’); ylabel(‘Amplitude’);

16.20

DECIMATOR AND INTERPOLA INTE RPOLATOR TOR

% Program or downsampling and upsampling the sum o two sinusoids using MATLAB’s inbuilt decimation and interpolation unction by a actor o 20. %Generate the input sequence or F s5200 Hz,  1550 Hz and  25100 Hz t50:1/200:10; y53.*cos(2*pi*50.*t/200)11.*cos(2*pi*100.*t/200); gure(1) stem(y); xlabel({‘Times in Seconds’;’(a)}); ylabel(‘Amplitude’);

883

MATLAB Programs

%Generate the decimated and interpolated signals gure(2) stem(decimate(y,20)); xlabel({‘Times in Seconds’;’(b)}); ylabel(‘Amplitude’); gure(3) stem(interp(decimate(y,20),2)); xlabel({‘Times in Seconds’;’(c)}); ylabel(‘Amplitude’);

5 e d u t i l p m A

0

−5

0

500

1000

1500

2000

2500

Time in Seconds (a) 5 e d u t i l p m A

0

−5

0

20

40

60

80

100

120

Time in Seconds (b) 5 e d u t i l p m A

0

−5

0

50

100

150

200

250

Time in Seconds (c)

Fig. 16.40

16.21

( a) Input Sequence, (b) Decimated Sequence and ( c) Interpolated s equence

ESTIMATION OF POWER SPECTRAL DENSITY (PSD)

% Program or estimating PSD o two sinusoids plus noise % % % %

Algorithm; 1:Get the requencies o the two sinusoidal waves 2:Get the sampling requency 3:Get the length o the sequence to be considered

884


% 4:Get 4:Get the two FFT lengths or comparing the corresponding power spectral densities clc; close all; clear all; 15input(‘Enter the requency o rst sinusoid’); 25input(‘Enter the requency o second sinusoid’); s5input(‘Enter the sampling requency’); N5input(“Enter the length o the input sequence’); N15input(“Enter the input FFT length 1’); N25input(“Enter the input FFT length 2’);

%Generation o input sequence t50:1/s:1; x52*sin(2*pi*1*1)13*sin(2*pi*2*t)2randn(size(t));

%Generation o psd or two dierent FFT lengths Pxx15abs(t(x,N1)).^2/(N11); Pxx25abs(t(x,N2)).^2/(N11); %Plot the psd; subplot(2,1,1); plot ((0:(N121))/N1*s,10*log10(Pxx1)); xlabel(‘Frequency in Hz’); ylabel(‘Power spectrum in dB’); title(‘[PSD with FFT length,num2str(N1)]’); subplot (2,1,2); plot ((0:(N221))/N2*s,10*log10(Pxx2)); xlabel(‘Frequency in Hz’); ylabel(‘Power spectrum in dB’); title(‘[PSD with FFT length,num2str(N2)]’);

16.22

PSD ESTIMA ESTI MATOR TOR

% Program or estimating PSD o a two sinusoids plus noise using %(i)non-overlapping sections %(ii)overlapping sections and averaging the periodograms clc; close all; clear all; 15input(‘Enter the requency o rst sinusoid’); 25input(‘Enter the requency o second sinusoid’); s5input(‘Enter the sampling requency’); N5input(“Enter the length o the input sequence’); N15input(“Enter the input FFT length 1’); N25input(“Enter the input FFT length 2’);

%Generation o input sequence t50:1/s:1; x52*sin(2*pi*1*1)13*sin(2*pi*2*t)2randn(size(t));

MATLAB Programs

885

%Generation o psd or two dierent FFT lengths Pxx15(abs(t(x(1:256))).^21abs(t(x(257:512))).^21 abs(t(x(513:768))).^2/(256*3); %using nonoverlapping sections Pxx25(abs(t(x(1:256))).^21abs(t(x(129:384))).^21ab s(t(x(257:512))).^21abs(t(x(385:640))).^21abs(t( x(513:768))).^21abs(t(x(641:896))).^2/(256*6); %using overlapping sections % Plot the psd; subplot (2,1,1); plot ((0:255)/256*s,10*log10(Pxx1)); xlabel(‘Frequency in Hz’); ylabel(‘Power spectrum in dB’); title(‘[PSD with FFT length,num2str(N1)]’); subplot (2,1,2); plot ((0:255)/256*s,10*log10(Pxx2)); xlabel(‘Frequency in Hz’); ylabel(‘Power spectrum in dB’); title(‘[PSD with FFT length,num2str(N2)]’);

16.23

PERIODOGRAM ESTIMATION

% Periodogram estimate cold be done by applying a nonrectangular data windows to the sections prior to computing the periodogram % This program estimates PSD or the input signal o two sinusoids plus noise using Hanning window 15input(‘Enter the requency o rst sinusoid’); 25input(‘Enter the requency o second sinusoid’); s5input(‘Enter the sampling requency’); t50:1/s:1; 5hanning(256); w x52*sin(2*pi*1*t)13*sin(2*pi*2*t)2randn(size(t)); Pxx5(abs(t(w.*x(1:256))).^21abs(t(w.*x(129:384))).^ 21abs(t(w.*x(257:512))).^21abs(t(w.*x(385:640))).^2 1abs(t(w.*x(513:768))).^21abs(t(w.*x(641:896))).^2/ (norm(w)^2*6); Plot((0:255)/256*s,10*log10(Pxx));

16.24

WELCH PSD ESTIMATOR

% Program or estimating the PSD o sum o two sinusoids plus noise using Welch Welch method n50.01:0.001:.1; x5sin(.25*pi*n)13*sin(.45*pi*n)1rand(size(n)); pwelch(x)

886


) 5 e l p m a s / 0 d a r / B d ( y −5 t i s n e D −10 m u r t c e p −15 S r e w o P −20

Welch PSD Estimate

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

Normalised Frequency (x pi rad/sample)

Fig. 16.41

Welch PSD Estimate

xlabel(‘Normalised Frequency (x pi rad/sample)’); ylabel(‘Power Spectrum Density (dB/rad/sample)’)

16.25

WELCH PSD ESTIMATOR USING WINDOWS

% Program or estimating the PSD o sum o two sinusoids using Welch method with an overlap over lap o 50 percent and with Hanning, Hamming, Bartlett, Blackman and rectangular windows. s51000; t50:1/s:3; x5sin(2*pi*200*t)1sin(2*pi*400*t); gure(1) subplot(211) pwelch(x,[],[],[],s); title(‘Overlay plot o 50 Welch estimates’) subplot(212) pwelch(x,hanning(512),0,512,s) title(‘N5512 Overlap550% Hanning’) gure(2) subplot(211) pwelch(x,[],[],[],s); title(‘Overlay plot o 50 Welch estimates’) subplot(212) pwelch(x,hamming(512),0,512,s) title(‘N5512 Overlap550% Hamming’) gure(3) subplot(211) pwelch(x,[],[],[],s); title(‘Overlay plot o 50 Welch estimates’) subplot(212) pwelch(x,bartlett(512),0,512,s) title(‘N5512 Overlap550% Bartlett’) gure(4) subplot(211) pwelch(x,[],[],[],s);

MATLAB Programs

title(‘Overlay plot o 50 Welch estimates’) subplot(212) pwelch(x,blackman(512),0,512,s) title(‘N5512 Overlap550% Blackman’) gure(5) subplot(211) pwelch(x,[],[],[],s); title(‘Overlay plot o 50 Welch estimates’) subplot(212) pwelch(x,boxcar(512),0,512,s) title(‘N5512 Overlap550% Rectangular’) ) z H / B d ( y t i s n e D l a r t c e p S r e w o P

Overlay plot of 50 Welch estimates

0 −20 −40 −60 −80

0

) z 0 H / B d ( y t i s −50 n e D l a r t c −100 e p S r e w o −150 P 0

50

100

150

200 250 300 Frequency (Hz)

350

400

450

500

400

450

500

N=512 Overlap = 50% Hanning

50

100

150


350

Fig. 16.42 ( a a) Welch Estimate with N 5512 , 50% Overlap Overlap Hanning Hanning ) z H / B d ( y t i s n e D l a r t c e p S r e w o P

) z H / B d ( y t i s n e D l a r t c e p S r e w o P


0 −20 −40 −60 −80

0

50

100

150


350

400

450

500

400

450

500

N=512 Overlap = 50% Hamming

0 −20 −40 −60 −80

0

50

100

150


350

Fig. 16.42 (b) Welch Estimate with N 5512 , 50% Overlap Overlap Hamming Hamming

887

888


) z H / B d ( y t i s n e D l a r t c e p S r e w o P ) z H / B d ( y t i s n e D l a r t c e p S r e w o P


0 −20 −40 −60 −80

0

50

100

150


350

400

450

500

400

450

500

N=512 Overlap = 50% Bartlett

0 −20 −40 −60 −80

0

50

Fig. 16.42

100

150


350

( c c) Welch Estimate with N 5512 , 50% Overlap Overlap Bartlett Bartlett

) z 0 H / B d ( y−20 t i s n e D−40 l a r t c e p−60 S r e w o−80 P 0

50

100

150

) z 0 H / B d ( y −50 t i s n e D l −100 a r t c e p−150 S r e w o−200 P 0

50

100

150


Fig. 16.42

200 250 300 350 Frequency (Hz) N=512 Overlap = 50% Blackman


350

400

450

500

400

450

500

d) Welch Estimate with N 5512 , 50% Overlap Overlap Blackman Blackman ( d

MATLAB Programs

) z 0 H / B d ( y−20 t i s n e D−40 l a r t c e p−60 S r e w o−80 P 0

889


50

) z 0 H / B d ( −10 y t i s−20 n e D−30 l a r t c−40 e p S−50 r e w o−60 P 0

100

150 N =

50

100

150


350

400

450

500

400

450

500

512 Overlap = 50% Rectangular

200 250 3 00 Frequency Frequency (Hz)

350

Fig. 16.42 ( e e) Welch Estimate with N 5512 , 50% Overlap Overlap Rectangu Rectangular lar

16.26

WELCH PSD ESTIMATOR USING WINDOWS

% Program or estimating the PSD o sum o two sinusoids plus noise using Welch Welch method with an overlap o 50 percent and with w ith Hanning, Hamming, Bartlett, Blackman and rectangular windows w indows s51000; t50:1/s:3; x52*sin(2*pi*200*t)15*sin(2*pi*400*t); y5x1randn(size(t)); gure(1) subplot(211); pwelch(y,[],[],[],s); title(‘Overlay plot o 50 Welch estimates’); subplot(212); pwelch(y,hanning(512),0,512,s); title(‘N5512 Overlap550% Hanning’); gure(2) subplot(211); pwelch(y,[],[],[],s); title(‘Overlay plot o 50 Welch estimates’); subplot(212); pwelch(y,hamming(512),0,512,s); title(‘N5512 Overlap550% Hamming’);

890


gure(3) subplot(211); pwelch(y,[],[],[],s); title(‘Overlay plot o 50 Welch estimates’); subplot(212); pwelch(y,bartlett(512),0,512,s); title(‘N5512 Overlap550% Bartlett’); gure(4) subplot(211); pwelch(y,[],[],[],s); title(‘Overlay plot o 50 Welch estimates’); subplot(212); pwelch(y,blackman(512),0,512,s); title(‘N5512 Overlap550% Blackman’); gure(5) subplot(211); pwelch(y,[],[],[],s); title(‘Overlay plot o 50 Welch estimates’); subplot(212); pwelch(y,boxcar(512),0,512,s); title(‘N5512 Overlap550% Rectangular’);

) z H / 10 B d ( 0 y t i s n e −10 D l a r t −20 c e p S −30 r e w o −40 P 0 ) z H / 10 B d ( 0 y t i s n e −10 D l a r t −20 c e p S −30 r e w o −40 P 0


50

10 0

1 50


3 50

400

450

500

400

45 0

500


50

100

150


350

Fig. 16.43 ( a a) Welch Estimate with N 5512, 50% Overlap Hanning

MATLAB Programs

) z 10 H / B d ( 0 y t i s n e −10 D l a r −20 t c e p S −30 r e w o −40 P 0


50

) z H / 10 B d ( 0 y t i s n e −10 D l a r t −20 c e p S −30 r e w o −40 P 0

Fig. 16.43

) z H / 10 B d ( 0 y t i s n e −10 D l a r t −20 c e p S −30 r e w o −40 P 0 ) z H / 10 B d ( 0 y t i s n e −10 D l a r t −20 c e p S −30 r e w o −40 P 0

100

150

200 250 30 0 Frequency (Hz)

350

400

450

500

400

450

500


50

100

150

2 00 250 300 Frequency (Hz)

350

(b) Welch Estimate with N 5512, 50% Overlap Hamming


50

100

150


350

400

450

500

400

450

500

N=512 Overlap = 50% Bartlett

50

100

150


350

Fig. 16.43 ( c c) Welch Estimate with N 5512, 50% Overlap Bartlett

891

892


) z H / 10 B d ( 0 y t i s n e −10 D l a r −20 t c e p S −30 r e w −40 o 0 P

50

100

150

400

450

500

) z H / 10 B d ( 0 y t i s n e −10 D l a r −20 t c e p S −30 r e w −40 o 0 P

200 250 300 350 Frequency (Hz) N=512 Overlap = 50% Blackman

50

100

150

400

450

500



350

( d d) Welch Estimate with N 5512, 50% Overlap Blackman

Fig. 16.43




50

100

150


350

400

450

500

400

450

500


50

100

150

200 250 30 0 Frequency Frequency (Hz)

350

Fig. 16.43 ( e e) Welch Estimate with N 5512, 50% Overlap Rectangular

MATLAB Programs

16.27

893

STA STATE-SPACE TE-S PACE REPRES RE PRESENTA ENTATION TION

% Program or computing the state-space matrices rom the given transer unction unction [A,B,C,D]5t2ss(b,a); a5input (‘enter the denominator polynomials5’); b5input (‘enter the numerator polynomials5’); p5length(a)21;q5length(b)21;N5 max(p,q); i(Np),a5[a,zeros(1,N2p)];end i(Nq),b5[b,zeros(1,N2q)];end A5[2a(2:N11);[eye(N21),zeros(N21,1)]]; B5[1;zeros(N21,1)]; C5b(2:N11)2b(1)*(2:N11); D5b(1);

16.2 16.28 8

PARTI ARTIAL AL FRAC FRACTI TION ON DECO DECOMP MPOS OSIT ITIO ION N

% Program or partial raction decomposition o a rational transer unction unction[c,A,alpha]5t2p(b,a); a5input (‘enter the denominator polynomials5’); b5input (‘enter the numerator polynomials5’); p5length(a)21; q5length(b)21; a5(1/a(1))*reshape(a,1,p11); b5(1/a(1))*reshape(b,1,q11); i(q5p),%case o nonempty c(z) temp5toeplitz([a,zeros(1,q2p)]’,[a(1),zeros(1,q2p)]); temp5[temp,[eye(p);zeros(q2p11,p)]); temp5temp/b’; c5temp(1:;q2p11); d5temp(q2p12:q11)’; else c5[]; d5[b,zeros(1,p2q21)]; end alpha5cplxpair (roots(a));’; A5zeros(1,p); or k51 :p temp5prod(alpha(k)2alpha(nd(1:p 5k))); i(temp550),error(‘repeated roots in TF2PF’); else,A(k)5polyval(d,alpha(k))/temp; end end 

894


16.2 16.29 9

INVE NVERS RSE E z-T z-TRANS RANSF FORM ORM

% Program or computing inverse z-transorm o a rational transer unction unction x5invz(b,a,N); b5input (‘enter the numerator polynomials5’); a5input (‘enter the denominator polynomials5’); N5input (‘enter the number o points to be computed5’); [c,A,alpha]5t2p(b,a); x5zeros (1,N); x(1:length(c))5c; or k51:length(A), x5x1A(k)*(alpha(k)).^(0:N21); end x5real(x);

16.30

GROUP DELAY

% Program or computing group delay o a rational transer unction on a given requency interval interva l unction D5grpdly(b,a,K,theta); b5input (‘enter the numerator polynomials5’); a5input (‘enter the denominator polynomials5’); K5input (‘enter the number o requency response points5’); theta5input (‘enter the theta value5’); a5reshape(a,1,length(a)); b5reshape(b,1,length(b)); i (length(a)551)%case o FIR bd52j*(0:length(b)21).*b; i(nargin553), B5rqresp(b,1,K); Bd5rqresp(bd,1,K); else, B5rqresp(b,1,K,theta); Bd5rqresp(bd,1,K,theta); end D5(real(Bd).*imag(B)2real(B).*imag(Bd))./abs(B).^2; else %case o IIR i(nargin553), D5grpdly (b,1,K)2grpdly(a,1,K); else, D5grpdly(b,1,K,theta)2grpdly(a,1,K,theta); end end

MATLAB Programs

16.31

895

IIR FILTER DESIGN-IMPULSE INVARIANT METHOD

% Program or transorming an analog flter into a digital flter using impulse invariant technique unction [bout,aout]5impinv(bin,ain,T); bin5input(‘enter the numerator polynomials5’); ain5input(‘enter the denominator polynomials5’); T5input(‘enter the sampling interval5’); i(length(bin)5length(ain)), error(‘Anlog lter in IMPINV is not strictly proper’); end [r,p,k]5residue(bin,ain); [bout,aout]5p2t([],T*r,exp(T*p));

16.32

IIR FILTER DESIGN-BILINEAR TRANSFORMATION

% Program or transorming an analog flter into a digial flter using bilinear transormation unction [b,a,vout,uout,Cout]5bilin(vin,uin,Cin,T); pin5input(‘enter the poles5’); zin5input(‘enter the zero5’); T5input(‘enter the sampling interval5’); Cin5input(‘enter the gain o the analog lter5’); p5length(pin); q5length(zin); Cout5Cin*(0.5*T)^(p2q)*prod(120.5*T*zin)/ prod(120.5*T*pin); zout5[(110.5*T*zin)./(120.5*T*pin),2ones(1,p2q)]; pout5(110.5*T*pin)./(120.5*T*pin); a51; b51; or k51 :length(pout),a5conv(a,[1,2pout(k)]); end or k51 :length(zout),b5conv(b,[1,2zout(k)]); end a5real(a); b5real(Cout*b); Cout5real(Cout);

16.33

DIRECT REALISATION OF IIR DIGITAL FILTERS

% Program or computing direct realisation values o IIR digital flter unction y5direct(typ,b,a,x); x5input(‘enter the input sequence5’);

896


b5input(‘enter the numerator polynomials5’); a5input(‘enter the denominator polynomials5’); typ5input(‘type o realisation5’); p5length(a)21; q5length(b)21; pq5 max(p,q); a5a(2:p11);u5zeros(1,pq);%u is the internal state; i(typ551) or i51:length(x),

5x(i)2sum(u(1:p).*a); unew u5[unew,u]; y(i)5sum(u(1:q11).*b); u5u(1:pq); end elsei(typ552) or i51:length(x) y(i)5b(1)*x(i)1u(1); u5u[(2:pq),0]; u(1:q)5u(1:q)1b(2:q11)*x(i); u(1:p)5u(1:p)2a*y(i); end end

16.34

PARALLEL REALISATION OF IIR DIGITAL FILTERS

% Program or computing parallel realisation values o IIR digital flter unction y5parallel(c,nsec,dsec,x); x5input(‘enter the input sequence5’); b5input(‘enter the numerator polynomials5’); a5input(‘enter the denominator polynomials5’); c5input(‘enter the gain o the lter5’); [n,m]5size(a);a5a(:,2:3); u5zeros(n,2); or i51:length(x), y(i)5c*x(i); or k51:n,

5x(i)2sum(u(k,:).*a(k,:));u(k,:)5[unew,u(k,1)]; unew y(i)5y(i)1sum(u(k,:).*b(k,:)); end end

MATLAB Programs

16.35

897

CASCADE REALISATION OF IIR DIGITAL FILTERS

% Program or computing cascade realisation values o digital IIR flter unction y5cascade(c,nsec,dsec,x); x5input(‘enter the input sequence5’); b5input(‘enter the numerator polynomials5’); a5input(‘enter the denomiator polynomials5’); c5input(‘enter the gain o the lter5’); [n,m]5size(b); a5a(:,2:3);b5b(:,2,:3); u5zeros(n,2); or i51 :length(x), or k51 :n, 5x(i)2sum(u(k,:).*a(k,:)); unew 1sum(u(k,:).*b(k,:)) x(i)52unew u(k,:)5[unew,u(k,1)]; end y(i)5c*x(i); end

16.36

DECIMATION BY POLYPHASE POLYPHASE DECOMPOSITION

% Program or computing convolution and m-old decimation by polyphase decomposition unction y5ppdec(x,h,M); x5input(‘enter the input sequence5’); h5input(‘enter the FIR lter coecients5’); M5input(‘enter the decimation actor5’); 1h5length(h); 1p5foor((1h21)/M)11; p5reshape([reshape(h,1,1h),zeros(1,1p*M21h)],M,1p); lx5length(x); ly5foor ((1x11h22)/M)11; 1u5oor((1x1M22)/M)11; %length o decimated sequences u5[zeros(1,M21),reshape(x,1,1x),zeros(1,M*lu2lx2M11)]; y5zeros(1,1u11p21); 51:M,y5y1conv(u(m,: ),p(m,: )); end or m y5y(1:1y);

16.3 16.37 7

MUL MULTIBA TIBAND ND FIR FIR FIL FILTER TER DESI DESIGN GN

% Program or the design o multiband FIR flters unction h5rdes(N,spec,win); N5input(‘enter the length o the lter5’);

898


spec5input(‘enter the low,high cuto requencies and gain5’); win5input(‘enter the window length5’); fag5rem(N,2); [K,m]5size(spec); n5(0:N)2N/2; i (˜fag),n(N/211)51; end,h5zeros(1,N11); or k51:K temp5(spec (k,3)/pi)*(sin(spec(k,2)*n)2sin(spec(k,1) *n))./n; i (˜fag);temp(N/211)5spec(k,3)*(spec(k,2)2 spec(k,1))/pi; end h5h1temp; end i (nargin553), h5h.*reshape(win,1,N11); end

16.3 16.38 8

ANAL ANALYS YSIS IS FIL FILTER TER BANK BANK

% Program or maximally decimated uniorm DFT analysis flter bank unction u5dtanal(x,g,M); g5input(‘enter the lter coecient5’); x5input(‘enter the input sequence5’); M5input(‘enter the decimation actor5’); 1g5length(g); 1p5foor((1g21)/M)11; p5reshape([reshape(g,1,1g),zeros(1,1p*M21g)],M,1p); lx5length(x); lu5foor ((1x1M22)/M)11; x5[zeros(1,M21),reshape(x,1,1x),zeros(1,M*lu2lx2M11)]; x5fipud(reshape(x,M,1u)); %the decimated sequences u5[];

51:M,u5[u; cov(x(m,:),p(m,:))]; end or m u5it(u);

16.3 16.39 9

SYNT SYNTHE HESI SIS S FIL FILTER TER BANK BANK

% Program or maximally decimated uniorm DFT synthesis flter bank unction y5udtsynth(v,h,M); 1h5length(h); ‘1q5foor((1h21)/M)11; q5fipud(reshape([reshape(h,1,1h),zeros(1,1q*M21h)],M,1q));

MATLAB Programs

899

v5t(v); y5[ ];

51:M,y5[conv(v(m,:),q(m,:));y]; end or m y5y(:).’;

16.4 16.40 0

LEVI LEVINS NSON ON-D -DUR URBI BIN N ALGO ALGORI RITH THM M

% Program or the solution o normal equations using Levinson Durbin algorithm unction [a,rho,s]5levdur(kappa); % Input; % kappa: covariance sequence values rom 0 to p % Output parameters: % a: AR polynomial,with leading entry 1 % rho set o p refection coecients % s: innovation variance kappa5input(‘enter the covariance sequence5’); p5length(kappa)21; kappa5reshape(kappa,p11,1); a51; s5kappa(1); rho5[]; or i51:p, rhoi5(a*kappa(i11:21:2))/s; rho5[rho,rhoi]; s5s*(12rhoi^2); a5[a,0]; a5a2rhoi*fiplr(a); end

16.4 16.41 1

WIEN WIENER ER EQUA EQUATI TION ON’S ’S SOLU SOLUTI TION ON

% Program unction b5 wiener(kappax,kappayx); kappax5input(‘enter the covariance sequence5’); kappyx5input(‘enter the joint covariance sequence5’); q5length(kappax)21; kappax5reshape(kappax,q11,1); kappayx5reshape(kappayx,q11,1); b5(toeplitz(kappax)/(kappayx)’;

16.42 16.42

SHOR SHORT-TIME -TIME SPECT SPECTRAL RAL ANAL ANALYSIS YSIS

% Program unction X5stsa(x,N,K,L,w,opt,M,theta0,dtheta); x5input(‘enter the input signal5’); L5input(‘enter the number consecutive DFTs to average5’); N5input(‘enter the segment length5’); K5input(‘enter 5input(‘enter the number o overlapping points5’); w

900


the window coecients5’); opt5input(‘opt5’); M5input(‘enter the length o DFT5’); theta05input(‘theta05’); dtheta5input(‘dtheta5’); 1x5length(x); nsec5ceil((1x2N)/(N2K)11; x5[reshape(x,1,1x),zeros(1,N1(nsec21)*(N2K))2lx)]; nout5N; i (nargin 5),nout5M; else,opt5‘n’; end X5zeros(nsec,nout); or n51: nsec, temp5 w.*x((n21)*(N2K)11:(n21)*(N2K)1N); i (opt(1) 55 ‘z’),temp5[temp,zeros(1,M2N)]; end i (opt(1)55‘c’),temp5chirp (temp,theta0,dtheta,M); else,temp5tshit(t(temp)); end X(n,: )5abs(temp).^2; end i(L1); nsecL5foor(nsec/L); or n51:nsecL,X(n,:)5 mean (X((n21)*L11:n*L,:)); end i (nsec55nsecL*L11), X(nsecL11,:)5X(nsecL*L11,:); X5X(1:nsecL11),: ); elsei(nsec nsecL*L), X(nsecL11,:)5 mean(x(nsecL*L11:nsec,:)); X5X(1:nsecL11,:); else,X5X(1:nsecL,:); end end LKh

16.43

CANCELLATION OF ECHO PRODUCED ON THE TELEPHONE-BASE BAND CHANNEL Base band transmit filter

Desired sequence

Echo Canceler

Echo path

Estimated sequence

Fig. 16.44 Baseband Channel Channel Echo Canceler Canceler % Simulation program or baseband echo cancellation shown in Fig. 16.44 using LMS algorithm clc; close all; clear all; ormat short T5input(‘Enter the symbol interval’); br5input(‘Enter the bit rate value’); r5input(‘Enter the roll o actor’); n5[210 10]; y55000*rcosr(r,n,br,T); %Transmit lter pulse shape is assumed as raised cosine

MATLAB Programs

901

ds5[5 2 5 2 5 2 5 2 5 5 5 5 2 2 2 5 5 5 5]; % data sequence 5length(ds); m nl5length(y); i51; z5conv(ds(i),y); while(i) z15[z, zeros(1,1.75*br)]; z5conv(ds(i11),y); z25[zeros(1,i*1.75*br),z]; z5z11z2; i5i11; end %plot(z); %near end signal h5randn(1,length(ds)); %echo path impulse response rs15lter(h,1,z); or i51; length(ds); rs(i)5rs 1(i)/15; end or i51: round(x3/3), rs(i)5randn(1); % rs2echo produced in the hybrid end s5[5 5 2 2 2 2 2 5 2 2 2 5 5 5 2 5 2 5 2]; % Desired data signal 5length(ds); m nl5length(y); i51; z5conv(s(i),y); while(i) z15[z,zeros(1,1.75*br)]; z5conv(s(i11),y); z25[zeros(1,i*1.75*br),z]; z5z11z2; i5i11; end s15rs1s; % echo added with desired signal ar5xcorr(ds,ds); crd5xcorr(rs,ds); ll5length(ar); j51; or i5round(11/2): 11, ar1(j)5ar(i); j5j11; end r5toeplitz(ar1); l25length(crd); j51; or i5round(l2/2):12, crdl(j)5crd(i); j5j11; end p5crd1’;

902


5 lam max(eig(r)); la5 min(eig(r)); l5lam/la; 5inv(r)*p; % Initial value o lter coecients w e5rs2lter(w,l,ds); s51; mu51.5/lam; ni51; while (s 1 e210) 22*mu*(e.*ds)’ ; % LMS algorithm adaptation w15 w rs y45lter(w1,1,ds); % Estimated echo signal using LMS algorithm e5y42rs; s50; e15xcorr(e); or i51:length(e1), s5s1e1(i); end s5s/length(e1); i (y455rs) break end ni5ni11; 5 w w1; end gure(1); subplot(2,2,1); plot(z); title(‘near end signal’); subplot(2,2,2); plot(rs); title(‘echo produced in the hybrid’); subplot(2,2,3); plot(s); title(‘desired signal’); subplot(2,2,4); plot(s1); title(‘echo added with desired signal’); gure(2); subplot(2,1,1); plot(y4); title(‘estimated echo signal using LMS algorithm’); subplot(2,1,2); plot(s12y4); title(‘echo cancelled signal’);

16.44

CANCELLATION OF ECHO PRODUCED ON THE TELEPHONE—PASS BAND CHANNEL

Pass band transmit filter

Desired sequence

Echo Canceler

Echo path

Passband receive filter

Estimated sequence

Fig. 16.45

Pass Band Channel Echo Canceler

MATLAB Programs

903

% Simulation program or passband echo cancellation shown in Fig. 16.45 using LMS algorithm clc; close all; clear all; ormat long d58000; s516000; c58000; 54000; t50:.01:1; %d5sin(2*pi**t/d); % Near end signal ns5[5 2 5 2 5 5 2 2 2 5 5 2 2 2 2 2 2 5 5 2 5 2 5 5 5 5 5 5 5 5 5 5 5 5 5 5]; % Near end input signal is digitally modulated and plotted y5dmod(ns,c,d,s,‘psk’); subplot(2,2,1); plot(y); title(‘input signal’); xlabel(‘Time ———’); ylabel(‘Amplitude ———’); % Echo is generated due to mismatch o hybrid impedances h55*randn(1,length(ns));——— rsl5lter(h,1,y); % or i51; length(ns); % rsl(i)5rs6(i); % end or i51; length(ns); rs(i)5rs1(i); end subplot(2,2,2); plot(rs); title(‘noise signal’); xlabel(‘Time ———’); ylabel(‘Amplitude ———’); % Far end signal s15[5 5 2 5 2 5 2 5 2 5 2 5 2 5 5 5 5 5 5 5 5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 5]; % rs5sign(rs2); % Far end signal is digitally modulated and plotted z15dmod(s1,c,d,s,‘psk’); or i51:length(ns), z(i)5z1(i); end subplot(2,2,3); plot(z); title(‘ar-end signal’); xlabel(‘Time ———’); ylabel(‘Amplitude ———’); % Echo and the ar end modulated signal is added in the hybrid q15z11rs1; or i51; length(ns); q(i)5q1(i);; end subplot(2,2,4); plot(q); title(‘received signal’); xlabel(‘Time ———’); ylabel(‘Amplitude ———’); q25xcorr(q); % Auto correlation is taken or the near end signal ar5xcorr(ns); % cross correction is taken or the near end and ar end signal crd5xcorr(rs,ns);

904


l15length(ar); j51; or i5round(ll/2): l1, ar1(j)5ar(i) j5j11; end % Toeplitz matrix is taken or the auto correlated signal r5toeplitz(ar1); l25length(crd); j51; or i5round(l2/2):l2, crd1(j)5crd(i); j5j11; end p5crd1’; % Maximum and minimum eigen values are calculated rom the toeplitz matrix 5 lam max(eig(r)); la5 min(eig(r)); l5lam/la; % initial lter taps are ound using the below relation 5inv(r)*p; w % The step size actor is calculated 5length(ns)22.5; m 2.95367)/.274274; a5(m mu5a/lam; % The initial error is calculated s51; e5rs2lter(w,1,ns); ni51; gure(2); subplot(2,2,1); % Filter taps are iterated until the mean squared error becomes E225 while (‘s25’ ! s0’) 22*mu*(e.*ns)’; w15 w i (ni5100) break; end rs y45lter(w1,1,ns) e5y42rs; s50; el5e.*e; or i51: length(e1), s5s1e1(i); end s5s/length(e1); ni5ni11; 5 w w1; plot (ni,e); hold on; title(‘ MSE vs no. o iterations’); end end subplot(2,2,2); plot(y4); title(‘estimated noise signal’); xlabel(‘Time ———’); ylabel(‘Amplitude ———’); subplot(2,2,3); plot(q-y4); title(‘noise cancelled signal’); xlabel(‘Time ———’); ylabel(‘Amplitude ———’);

MATLAB Programs

905

Review Questions 16.1

(a) Generate unit impulse unction (b) Generate signal x signal x((n)5u(n) 2 u(n 2 N ) (c) Plot the sequence, x sequence, x((n)5 A cos ((2p f n f n)/ f f s ), where n50 to 100 f s5100 Hz, f Hz, f 51 Hz, A Hz, A50.5 (d) Generate x Generate x (n)5exp (25n), where n50 to 10.

16.2

Consider a system with impulse response (1 / 2) n ,   h( n) =    0,

16.3

n=0

to 4 elsewhere.

Consider the fgure.

( )

h 1 (n )

y ( y (n )

h 2 (n ) (−)

h 3 (n )

h 4 (n )

Fig. Q16.3

(a) Express the overall impulse impulse response in terms o h o h1(n), h2(n), h3(n), h4(n) (b) Determine h(n) when h1(n)5{1/2, 1/4, 1/2} h2(n)5h3(n)5(n 1 1) u(n) h4(n)5d(n 2 2) (c) Determine the response o the system in part (b) i x i x((n)5 d(n 1 2) 1 3d(n 2 1) 2 4d(n 2 3). 16.4

Compute the overall impulse response o the system ( )

H 1

3 (n )

H 2

H 3

H 4

y ( y (n )

y 3 (n )

Fig. Q16.4

or 0 or 0 # n # 99. 99. The system H system H 1 , H 2 , H 3 , H 4 , are specifed by H 1 : h1[n]5{1, 1/2, 1/4, 1/8, 1/16, 1/32} H 2 : h2[n]5{1, 1, 1, 1, 1} H 3 : y3[n]5(1/4) x (1/4) x((n) 1 (1/2) x (1/2) x((n 2 1) 1 (1/4) x (1/4) x((n 2 2) H 4 : y[ y[n]50.9 y 0.9 y((n 2 1) 2 0.81 y 0.81 y((n 2 2) 1 V (n) 1 V (n 2 1) Plot h(V ) or 0 or 0 # n # 99. 16.5

Consider the system with h(n)5an u(n), 21 , a , 1. Determine the response. x( x(n)5u(n 1 5) 2 u(n 2 10)

906


y ( y (n )

h ( h (n ) (−)

z −2

h (n )

Fig. Q16.5

16.6

(i) Determine the range o values o the parameter ‘a ‘ a’ or which the linear time invariant system with impulse response n   a , n ≥ 0, n even h( n) =    0, otherwise

is stable. stable. (ii) Determine the response o the system with impulse response h(n)5an u(n) to the input signal x(n)5u(n) 2 u(n 2 10) 16.7 Consider the system described by the dierence equation y(n)5a ( y y 1 1) 1 bx (n). Determine ‘b ‘ b’ in terms o a so that S h(n)51.

16.8

16.9

(a) Compute the the zero-state zero-state step response s response s((n) o the system and choose ‘b ‘ b’ so that s that s((`)51. (b) Compare the values o ‘ b’ obtained in parts (a) and (b). What did you observe? Compute and sketch the convolution y( y(n) and correlation r xh(n) sequences or the ollowing pair o signals and comment on the results obtained.       1, 2, 4  1, 1, 1, 1, 1  x1 ( n) =  h n = ( )    1         ↑  ↑        0, 1, − 2, 3, − 4  1 / 2, 1, 2, 1, 1 / 2  x2 ( n) =  h2 ( n) =    ↑ ↑                 1, 2, 3, 4   4, 3, 2, 1   x3 ( n) =  h3 ( n) =            ↑  ↑      1 , 2 , 3 , 4 1 , 2 , 3 , 4         x4 ( n) =  h4 ( n) =            ↑  ↑  Consider the recursive discrete-time system described by the dierence equation. y(n)5a1 y( y(n 2 1) 2 a2(n 2 2) 1 b0 x( x(n) where a1520.8, a250.64, b050.866 (a) Write a program to compute compute and plot the impulse response response h(n) o the system or 0 # n # 49. 49. (b) Write a program to compute compute and plot the zero-state zero-state step response s( s(n) o the system or 0 # n # 100. (c) Defne an FIR system with impulse response response h FIR(n) given by  h( n),  h FIR ( n) = 

0 ≤ n ≤ 19 elsewhere   0, where h(n) is the impulse response computed in part (a). Write a program to compute and plot its step response. response. (d) Compare the results obtained obtained in part (b) and part part (c) and explain their similarities and dierences.

MATLAB Programs

16.10

16.11

Determine and plot the real and imaginary imaginary parts parts and the magnitude magnitude and phase spectra o the ollowing DTFT or various values o r and u G( z) z)51/(1 2 2r (cos (cos u) z21 1 r 2 z22) or 0 , r , 1. Using MATLAB program compute the circular convolution o two length-N sequences via the DFT based approach. Using this problem determine the circular convolution o the ollowing pairs o sequences: sequences : (a) x(n)5{1, 2 3, 4, 2, 0, 2 2}

h(n)5{3, 0, 1, 2 1, 2, 1}

(b) x(n)5{3 1 j2, j2, 22 1 j, j, j3, j3,

h(n)5{ 1 2 j3, j3, 4 1 j2, j2, 2 2 j2, j2,

1 1 j4, j4, 23 1 j3} j3} , 16.12

907

23 1 j5, j5, 2 1 j }

(c) x(n)5cos (pn/ 2) 2) h(n)53n 0 # n # 5 Determine the actored orm o the ollowing z ollowing z-transorms -transorms (a) H 1( z) z)5(2 z4 1 16 z3 1 44 z2 1 56 z 1 32) / / (3 z3 1 3 z3 2 15 z2 1 18 z 2 12) (b) H 2( z) z)5(4 z4 2 8.68 z3 2 17.98 z2 1 26.74 z 2 8.04)/

16.13

( z z4 2 2 z3 1 10 z2 1 6 z 1 65) and show their pole-zero plots. Determine all regions o convergence o each o the above z above z-transorms, -transorms, and describe the type o their inverse z-transorm z-transorm (let-sided, right-sided, two-sided sequences) associated with each o the ROCs. Determine the z the z-transorm -transorm as a ratio o two polynomials in z21 rom each o the partial-raction expansions listed below: (a)

(b)

H 2 ( z ) = 3 +

(c)

H 3 ( z ) =

(d) 16.14 16.15

16.16

H1 ( z ) = 3 +

12 (2 − z −1 )

−

16 ( 4 − z −1 )

−1

20

H 4 ( z ) = 8 +

( 4 − z 1 )

−

(1 + 0.5 − z )

(5 + 2 − z )

z > 0.5

−

3

−1

,

2

(5 + 2 z −1 )

(1 + 0.25 z ) 10

−

10

−2

−1

+

(5 + 2 z )

+

z

,

z > 0.5

4 (1 + 0.9 z −2 )

,

z > 0.4

−1

(6 + 5 z −1 + z −2 )

,

z > 0.4

Determine the inverse z inverse z-transorm -transorm o each z each z-transorm -transorm given in Q16.13. Q16.13 . Consider the system − − − (1− 2 z 1 + 2 z 2 − z 3 ) H ( z ) = , ROC 0.5 < z < 1 (1− z −1 )(1 − 0.5 z −1 )(1 − 0.2 z −1 ) (a) Sketch the pole-zero pattern. Is the system stable? (b) Determine the impulse response o the system. Determine the impulse impulse response and the step response o the ollowing causal systems. Plot the pole-zero patterns and determine which o the systems are stable. stable . 3 1 (a) y( n) = y( n −1) − y( n − 2) + x( n) 4 8

908


(b) y( y(n)5 y (n 2 1) 2 0.5y( 0.5y(n 2 2) 1 x(n) 1 x (n 2 1) (c)

H ( z ) =

z

−1

(1 + z −1 )

(1− z )3 (d) y( y(n)50.6 y( y(n 2 1) 2 0.08 y( y(n 2 2) 1 x( x(n) (e) y( y(n)50.7 y( y(n 2 1) 2 0.1 y( y(n 2 2) 1 2 x( x(n) 2 x( x(n 2 2) Ans:

16.17

(a), (b), (d) and (e) are stable, (c) is unstable

The requency requency analysis o an amplitude-modulated amplitude-modulated discrete-time signal x( x(n)5sin 2p f 1n 1 sin 2p f 2 n where f 1

1 =

128

and f 2

5 =

128

modulates the amplitude-modulated siganl is

xc(n) 5 sin 2p f c n where f where f c 550/128. The resulting amplitude-modulated signal is xam(n)5 x( x(n) sin 2p f c n (a) Sketch the signals x signals x((n), x ), xc(n) and x and xam(n), 0 # n # 255 (b) Compute and sketch the 128-point DFT o the signal xam(n), 0 # n # 127 (c) Compute and sketch the 128-point DFT o the signal signal x xam(n), 0 # n # 99 (d) Compute and sketch the 256-point DFT o the signal xam(n), 0 # n # 179 (e) Explain the results obtained obtained in parts (b) through (d) by deriving deriving the spectrum o the amplitude modulated signal and comparing it with the experimental results. 16.18 A continuous time signal xa(t ) consists o a linear combination o sinusoidal signals o requencies 300Hz , 300Hz , 400Hz, 1.3kHz, 3.6KHz and 4.3KHz. The xa(t ) is sampled at 4kHz rate and the sampled sequence is passed through an ideal low-pass flter with cut o requency o 1kHz, generating a continuous time signal y signal ya(t ). ). What are the requency components present in the reconstructed signal y signal ya(t )? )? 16.19 Design an FIR linear phase, phase, digital digital flter flter approximating approximating the ideal ideal requency requency response 1,  for ␻ ≤ ␲ / 6  H d ( ␻) =   0, for ␲ / 6 < ␻ ≤ ␲  

16.20

(a) Determine the coefcient coefcient o a 25 tap flter flter based on the window method with a rectangular window. (b) Determine and plot plot the magnitude and phase phase response o the flter. flter. (c) Repeat parts (a) and (b) using the Hamming Hamming window (d) Repeat parts (a) and (b) (b) using the Bartlett window. Design an FIR FIR Linear Phase, bandstop flter having the ideal ideal requency requency response   1, for ␻ ≤ ␲ / 6    H d ( ␻) = 0, for ␲ / 6 < ␻ ≤ ␲ / 3   1, for ␲ / 3 ≤ ␻ ≤ ␲    (a) Determine the coefcient o a 25 tap flter based on the window method with a rectangular window. (b) Determine and plot plot the magnitude and phase phase response o the flter. flter.

MATLAB Programs

16.21

16.22

16.23

909

(c) Repeat parts (a) and (b) using the Hamming Hamming window (d) Repeat parts (a) and (b) using the Bartlett Bartlett window. A digital low-pass low-pass flter is required to meet the ollowing ollowing specfcations specfcations Passband ripple # 1 dB Passband edge 4 kHz Stopband attenuation  40 dB Stopband edge 6 kHz Sample rate 24 kHz The flter is to be designed by perorming a bilinear transormation on an analog system unction. Determine what order Butterworth, Chebyshev and elliptic analog design must be used to meet the specifcations in the digital implementation. An IIR digital low-pass low-pass flter is required to meet the ollowing ollowing specfcations specfcations Passband ripple # 0.5 dB Passband edge 1.2 kHz Stopband attenuation  40 dB Stopband edge 2 kHz Sample rate 8 kHz Use the design ormulas to determine the flter order or (a) Digital Butterworth flter (b) Digital Chebyshev flter (c) Digital elliptic flter An analog signal o the orm x orm xa(t )5a(t ) cos(2000 pt ) is bandlimited to the range 900 # F # 1100Hz. It is used as an input to the system shown in Fig. Q16.23. (t )

a

A /D R

= 2500

( )

ϖ

H ( ω)

(n ) D /A

a (n )

cos (0.8 πn )

Fig. Q16.23

16.24

(a) Determine and and sketch the spectra or the signal x signal x((n) and w(n). (b) Use Hamming window window o length M 531 to design a low-pass linear phase FIR flter H H () that passes {a { a(n)}. (c) Determine the sampling sampling rate o A/D converter converter that would allow us to eliminate the requency conversion in the above fgure. Consider the signal x signal x((n)5an u(m), |a |a| , 1 (a) Determine the spectrum spectrum X X () (b) The signal x signal x((n) is applied to a device (decimator) which reduces the rate by a actor o two. Determine the output spectrum. (c) Show that the spectrum spectrum is simply the Fourier transorm o x o x(2 (2n n).

16.25

Design a digital type-I Chebyshev low-pass flter operating at a sampling rate o 44.1kHz with a passband requency at 2kHz, a pass band ripple o 0.4dB, and a minimum stopband attenuation o 50dB at 12kHz using the impulse invariance method and the bilinear transormation method. Determine the order o analog flter prototype and design the analog prototype flter. Plot the gain and phase responses o the both designs using

910


MATLAB. Compare the perormances o the two flters. Show all steps used in the design. Hint 1. The order o flter cosh 1 ( ( A2 −

N

1) / ␧

−

=

cosh 1 ( ⍀ s / ⍀p) 2. Use the unction cheblap. cheblap. −

16.26

Design a linear linear phase phase FIR high-pass flter with ollowing specifcations: Stopband edge at 0.5p, passband edge at 0.7p, maximum passband attenuation o 0.15dB and a minimum stopband attenuation o 40dB. Use each o the ollowing windows or the design. Hamming, Hanning, Blackman and Kaiser. Show the impulse response coefcients and plot the gain response o the designed flters or each case.

16.27

Design using the windowed Fourier series approach a linear phase FIR low pass flter with the ollowing specifcations: pass band edge at 1 rad/s, stop band edge at 2rad/s, maximum passband attenuation o 0.2dB, minimum stopband attenuation o 50dB and a sampling requency o 10rad/s. Use each o the ollowing windows or the design: Hamming, Hanning, Blackman, Kaiser and Chebyshev. Show the impulse response coefcients and plot the gain response o designed flters or each case.

16.28

Design a two-channel crossover FIR low-pass low-pass and high-pass flter pair or digital audio applications. The low-pass and high-pass flters are o length 31 and have a crossover requency o 2kHz 2 kHz operating at a sampling rate o 44.1 KHz. KHz. Use the unction ‘fr1’ with a Hamming window to design the low pass flter while the high-pass flter is derived derived rom the low-pass flter using the delay complementary property. Plot the gain responses o both flters on the same fgure. What is the minimum number o delays and multipliers needed to implement the crossover network?

16.29

Design a digital network network butterworth butterworth low-pass flter operating operating at sampling rate o 44.1kHz with a 0.5dB cuto requency at 2kHz and a minimum stop band attenuation o 45dB at 10kHz using the impulse invariance method and the bilinear transormation method. Assume the sampling interval or the impulse invariance design to be equal to 1. Determine the order o the analog flter prototype and then design the analog prototype flter. Plot the gain and phase responses o both designs. Compare the perormances o the flters. Show all steps used in the design. Does the sampling interval have any eect on the design o the digital flter design based on the impulse invariance method? Hint The order o flter is log10 (1 / k 1 ) N

=

log10 (1 / k )

and use the unction ‘ buttap’.

Digital Signal Processing-chapter 16: MATLAB Programs

Recommend Documents