Chapter 16

INTRODUCTION

MATLAB stands for MATrix LABoratory. It is a technical computing environment for high performance 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 few dozen lines, to compute the solution with great accuracy in a few minutes on a computer, and to readily manipulate a three-dimensional display of 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 fraction of the time that it would take to write a program and execute in a language such as FORTRAN, BASIC, or C. It also features a family of application specific solutions, called toolboxes. Areas in which toolboxes are available include signal processing, image processing, control systems design, dynamic systems simulation, systems identification, 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 specific problems that are included in the previous chapters. All these MATLAB programs have been tested under version 7.1 of MATLAB and version 6.12 of the signal processing toolbox.

REPRESENTATION OF BASIC SIGNALS

MATLAB programs for the generation of unit impulse, unit step, ramp, exponential, sinusoidal and cosine sequences are as follows.

% Program for the generation of unit impulse signal

clc;clear all;close all; t522:1:2; y5[zeros(1,2),ones(1,1),zeros(1,2)];subplot(2,2,1);stem(t,y); ylabel('Amplitude --.'); xlabel('(a) n --.');

% Program for the generation of 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 for the generation of ramp sequence

n15input('enter the length of ramp sequence'); t50:n1; subplot(2,2,3);stem(t,t);ylabel('Amplitude --.'); xlabel('(c) n --.');

% Program for the generation of exponential sequence

n25input('enter the length of 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 for the generation of sine sequence t50:.01:pi; y5sin(2*pi*t);figure(2); subplot(2,1,1);plot(t,y);ylabel('Amplitude --.'); xlabel('(a) n --.'); % Program for the generation of 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 of ramp sequence 7 enter the length of exponential sequence 7 enter the a value 1

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

% Program for linear convolution of the sequence x5[1, 2] and h5[1, 2, 4]

clc; clear all; close all; x5input('enter the 1st sequence'); h5input('enter the 2nd sequence'); y5conv(x,h); figure;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(n)and h(n)and the convolved output signal y(n).

Circular Convolution

% Program for 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))]; for i = 1:N k = i; for j = 1:n2 H(i,j)=x(k)* b(j); k = k-1; if (k == 0) k = N; end end end y=zeros(1,N); M=H'; for j = 1:N for i = 1:n2 y(j)=M(i,j)+y(j); end end disp('The output sequence is y(n)= '); disp(y); if(j550) j5N1j; end y(n)5y(n)1g(i)*h(j); end end disp('The resultant signal is');y As an example, enter the first sequence [1 2 4] enter the 2nd sequence [1 2] The resultant signal is y51 4 8 8

Overlap Save Method and Overlap Add method

% Program for computing Block Convolution using Overlap Save Method

Overlap Save 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; h1= [h zeros(1,N-n1)]; n3=length(h1); y=zeros(1,N); x1= [zeros(1,n3-n2) x zeros(1,n3)]; H=fft(h1); for i=1:n2:N y1=x1(i:i+(2*(n3-n2))); y2=fft(y1); y3=y2.*H; y4=round(ifft(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 for 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=fft(h1); for i=1:n2:n1 if 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=fft(x1); x3=x2.*H; x4=round(ifft(x3)); if (i==1) 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

DISCRETE CORRELATION

Crosscorrelation

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

% Program for computing cross-correlation of 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); figure;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(fliplr(y));ylabel('Amplitude --.'); xlabel('(c) n --.'); disp('The resultant signal is');fliplr(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(n)and h(n)and the cross-correlated output signal y(n).

Autocorrelation

Algorithm 1. Get the signal x(n)of length N in matrix form 2. The correlated signal is denoted as y(n) 3. y(n)is given by the formula

% Program for computing autocorrelation function

x5input('enter the sequence'); y5xcorr(x,x); figure;subplot(2,1,1); stem(x);ylabel('Amplitude --.'); xlabel('(a) n --.'); subplot(2,1,2); stem(fliplr(y));ylabel('Amplitude --.'); xlabel('(a) n --.'); disp('The resultant signal is');fliplr(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(n)and its auto-correlated output signal y(n).

STABILITY TEST

% Program for stability test

clc;clear all;close all; b5input('enter the denominator coefficients of the filter'); k5poly2rc(b); knew5fliplr(k); s5all(abs(knew)1); if(s55 1) disp('"Stable system"'); else disp('"Non-stable system"'); end As an example, enter the denominator coefficients of the filter [1 21 .5] "Stable system"

SAMPLING THEOREM

The sampling theorem can be understood well with the following example. (2*pi*f1*n)1cos(2*pi*f2*n); xa5cos(2*pi*fc*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 for Section (b) n50:127;figure;n15128; f151/128;f255/128;fc550/128; x5cos(2*pi*f1*n)1cos(2*pi*f2*n); xc5cos(2*pi*fc*n); xa5cos(2*pi*fc*n); (Contd.) xamp5x.*xa;xam5fft(xamp,n1); stem(n,xam);title('xamp(n)');xlabel('n --.'); ylabel('amplitude'); %128 point DFT computation2solution for Section (c) n50:99;figure;n250:n121; f151/128;f255/128;fc550/128; x5cos(2*pi*f1*n)1cos(2*pi*f2*n); xc5cos(2*pi*fc*n); xa5cos(2*pi*fc*n); xamp5x.*xa; for i51:100, xamp1(i)5xamp(i); end xam5fft(xamp1,n1); s t e m ( n 2 , x a m ) ; t i t l e ( ' x a m p ( n ) ' ) ; x l a b e l ( ' n --.');ylabel('amplitude'); (a)Modulated signal x(n), carrier signal x a (n) and amplitude modulated signal x am (n) are shown in Fig. 16.5(a). Fig. 16.5 (b) shows the 128-point DFT of the signal x am (n) for 0 # n # 127 and Fig. 16.5 (c) shows the 128-point DFT of the signal x am (n), 0 # n # 99. The eight-point decimation-in-time fast Fourier transform of the sequence x(n)is computed using MATLAB program and the resultant output is plotted in Fig. 16 [h,om]5freqs(b,a,w); m520*log10(abs(h)); an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel('Gain in dB --.');xlabel('(a) Normalised frequency --.'); subplot(2,1,2);plot(om/pi,an); xlabel('(b) Normalised frequency --.'); ylabel('Phase in radians --.');

FAST FOURIER TRANSFORM

As an example, enter the passband ripple 0.15 enter the stopband ripple 60 enter the passband freq 1500 enter the stopband freq 3000 enter the stopband freq 7000 The amplitude and phase responses of the Butterworth low-pass analog filter are shown in Fig. 16.7.

High-pass Filter

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter using Eq. 8.46 5. Find the filter coefficients 6. Draw the magnitude and phase responses.

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter using Eq. 8.67 5. Find the filter coefficients 6. Draw the magnitude and phase responses. [h,om]5freqs(b,a,w); m520*log10(abs(h)); an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel('Gain in dB --.');xlabel('(a) Normalised frequency --.'); subplot(2,1,2);plot(om/pi,an); xlabel('(b) Normalised frequency --.'); ylabel('Phase in radians --.');

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter using Eq. 8.46 5. Find the filter coefficients 6. Draw the magnitude and phase responses.

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter using Eq. 8.57 5. Find the filter coefficients 6. Draw the magnitude and phase responses. [h,om]5freqz(b,a,w); m520*log10(abs(h)); an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel('Gain in dB --.');xlabel('(a) Normalised frequency --.'); subplot(2,1,2);plot(om/pi,an); xlabel('(b) Normalised frequency --.'); ylabel('Phase in radians --.');

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter using Eq. 8.67 5. Find the filter coefficients 6. Draw the magnitude and phase responses. (

% Program for the design of Butterworth analog high-pass filter

clc; close all;clear all; format long rp5input('enter the passband ripple'); rs5input('enter the stopband ripple'); wp5input('enter the passband freq'); ws5input('enter the stopband freq'); fs5input('enter the sampling freq'); w152*wp/fs;w252*ws/fs; [n,wn]5buttord(w1,w2,rp,rs,'s'); [b,a]5butter(n,wn,'high','s'); w50:.01:pi; [h,om]5freqs(b,a,w); m520*log10(abs(h)); an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel('Gain in dB --.');xlabel('(a) Normalised frequency --.'); subplot(2,1,2);plot(om/pi,an); xlabel('(b) Normalised frequency --.'); ylabel('Phase in radians --.');

As an example, enter the passband ripple 0.2 enter the stopband ripple 40 enter the passband freq 2000 enter the stopband freq 3500 enter the sampling freq 8000

The amplitude and phase responses of Butterworth high-pass analog filter are shown in Fig. 16.8.

Bandpass Filter

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter using Eq. 8.46 5. Find the filter coefficients 6. Draw the magnitude and phase responses.

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter using Eq. 8.67 5. Find the filter coefficients 6. Draw the magnitude and phase responses. [h,om]5freqs(b,a,w); m520*log10(abs(h)); an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel('Gain in dB --.');xlabel('(a) Normalised frequency --.'); subplot(2,1,2);plot(om/pi,an); xlabel('(b) Normalised frequency --.'); ylabel('Phase in radians --.');

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter using Eq. 8.57 5. Find the filter coefficients 6. Draw the magnitude and phase responses.

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequency 3. Get the sampling frequency 4. Calculate the order of the filter using Eq. 8.67 5. Find the filter coefficients 6. Draw the magnitude and phase responses.

% Program for the design of Butterworth analog Bandpass filter

clc; close all;clear all; format long rp5input('enter the passband ripple...'); rs5input('enter the stopband ripple...'); wp5input('enter the passband freq...'); ws5input('enter the stopband freq...'); fs5input('enter the sampling freq...'); w152*wp/fs;w252*ws/fs; The amplitude and phase responses of Butterworth bandpass analog filter are shown in Fig. 16.9. [h,om]5freqs(b,a,w); m520*log10(abs(h)); an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel('Gain in dB --.');xlabel('(a) Normalised frequency --.'); subplot(2,1,2);plot(om/pi,an); xlabel('(b) Normalised frequency --.'); ylabel('Phase in radians --.');

As an example, enter the passband ripple... 0.28 enter the stopband ripple... 28 enter the passband freq... 1000 enter the stopband freq... 1400 enter the sampling freq... 5000

The amplitude and phase responses of Butterworth bandstop analog filter are shown in Fig. 16.10.

16.9

CHEBYSHEV TYPE-1 ANALOG FILTERS

Low-pass Filter

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter using Eq. 8.57 5. Find the filter coefficients 6. Draw the magnitude and phase responses. The amplitude and phase responses of Chebyshev type -1 low-pass analog filter are shown in Fig. 16.11. The amplitude and phase responses of Chebyshev type -1 high-pass analog filter are shown in Fig. 16.12.

% Program for the design of Chebyshev Type-1 low-pass filter

Bandpass Filter Algorithm

1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter using Eq. 8.57 5. Find the filter coefficients 6. Draw the magnitude and phase responses.

% Program for the design of Chebyshev Type-1 Bandpass filter

clc; close all;clear all; format long rp5input('enter the passband ripple...'); rs5input('enter the stopband ripple...'); wp5input('enter the passband freq...'); ws5input('enter the stopband freq...'); fs5input('enter the sampling freq...'); w152*wp/fs;w252*ws/fs; [n]5cheb1ord(w1,w2,rp,rs,'s'); wn5[w1 w2]; [b,a]5cheby1(n,rp,wn,'bandpass','s'); w50:.01:pi; [h,om]5freqs(b,a,w); m520*log10(abs(h)); an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel('Gain in dB --.');xlabel('(a) Normalised frequency --.'); subplot(2,1,2);plot(om/pi,an); xlabel('(b) Normalised frequency --.'); ylabel('Phase in radians --.'); As an example, enter the passband ripple.. The amplitude and phase responses of Chebyshev type -1 bandpass analog filter are shown in Fig. 16.13.

Bandstop Filter

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequency 3. Get the sampling frequency 4. Calculate the order of the filter using Eq. 8.57 5. Find the filter coefficients 6. Draw the magnitude and phase responses.

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter using Eq. 8.67 5. Find the filter coefficients 6. Draw the magnitude and phase responses.

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter using Eq. 8.46 5. Find the filter coefficients 6. Draw the magnitude and phase responses.

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter using Eq. 8.67 5. Find the filter coefficients 6. Draw the magnitude and phase responses. (h) The amplitude and phase responses of Chebyshev type -2 bandstop digital filter are shown in Fig. 16.30.

% Program for the design of Chebyshev Type-1 Bandstop filter

clc; close all;clear all; format long rp5input('enter the passband ripple...'); rs5input('enter the stopband ripple...'); wp5input('enter the passband freq...'); ws5input('enter the stopband freq...'); The amplitude and phase responses of Chebyshev type -1 bandstop analog filter are shown in Fig. 16.14. The amplitude and phase responses of Chebyshev type -2 low-pass analog filter are shown in Fig. 16.15.

% Program for the design of Chebyshev Type-2 High pass analog filter

As an example, enter the passband ripple... 0.34 enter the stopband ripple... 34 enter the passband freq... 1400 enter the stopband freq... 1600 enter the sampling freq... 10000

The amplitude and phase responses of Chebyshev type -2 high-pass analog filter are shown in Fig. 16.16.

Gain in dB

Phase in radians

Phase in radians

Phase in radians

% Program for the design of Chebyshev

As an example, enter the passband ripple... 0.37 enter the stopband ripple... 37 enter the passband freq... 3000 enter the stopband freq... 4000 enter the sampling freq... 9000 The amplitude and phase responses of Chebyshev type -2 bandpass analog filter are shown in Fig. 16.17.

As an example, The amplitude and phase responses of Chebyshev type -1 high-pass digital filter are shown in Fig. 16.24.

% Program for the design of Chebyshev Type-2 Bandstop analog filter

clc; close all;clear all; format long rp5input('enter the passband ripple...'); rs5input('enter the stopband ripple...'); wp5input('enter the passband freq...'); ws5input('enter the stopband freq...'); fs5input('enter the sampling freq...'); w152*wp/fs;w252*ws/fs; [n]5cheb2ord(w1,w2,rp,rs,'s'); wn5[w1 w2]; [b,a]5cheby2(n,rs,wn,'stop','s'); w50:.01:pi; [h,om]5freqs(b,a,w); m520*log10(abs(h)); an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel('Gain in dB --.');xlabel('(a) Normalised frequency --.'); subplot(2,1,2);plot(om/pi,an); xlabel('(b) Normalised frequency --.'); ylabel('Phase in radians --.'); The amplitude and phase responses of Chebyshev type -2 bandstop analog filter are shown in Fig. 16.18. (h)); an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel('Gain in dB --.');xlabel('(a) Normalised frequency --.'); subplot(2,1,2);plot(om/pi,an); xlabel('(b) Normalised frequency --.'); ylabel('Phase in radians --.');

As an example, enter the passband ripple 0.5 enter the stopband ripple 50 enter the passband freq 1200 enter the stopband freq 2400 enter the sampling freq 10000

The amplitude and phase responses of Butterworth low-pass digital filter are shown in Fig. 16.19.

% Program for the design of Butterworth highpass digital filter

clc; close all;clear all; format long rp5input('enter the passband ripple'); rs5input('enter the stopband ripple'); wp5input('enter the passband freq'); ws5input('enter the stopband freq'); fs5input('enter the sampling freq'); w152*wp/fs;w252*ws/fs; [n,wn]5buttord(w1,w2,rp,rs); [b,a]5butter(n,wn,'high'); w50:.01:pi; [h,om]5freqz(b,a,w); m520*log10(abs(h)); an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel('Gain in dB --.');xlabel('(a) Normalised frequency --.'); subplot(2,1,2);plot(om/pi,an); xlabel('(b) Normalised frequency --.'); ylabel('Phase in radians --.'); As an example, enter the passband ripple 0.5 enter the stopband ripple 50 enter the passband freq 1200 The amplitude and phase responses of Butterworth high-pass digital filter are shown in Fig. 16.20.

Band-pass Filter

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter using Eq. 8.46 5. Find the filter coefficients 6. Draw the magnitude and phase responses.

% Program for the design of Butterworth Bandpass digital filter

clc; close all;clear all; format long rp5input('enter the passband ripple'); rs5input('enter the stopband ripple'); wp5input('enter the passband freq'); ws5input('enter the stopband freq'); fs5input('enter the sampling freq'); w152*wp/fs;w252*ws/fs; [n]5buttord(w1,w2,rp,rs); wn5[w1 w2]; [b,a]5butter(n,wn,'bandpass'); w50:.01:pi; [h,om]5freqz(b,a,w); m520*log10(abs(h)); an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel('Gain in dB --.');xlabel('(a) Normalised frequency --.'); subplot(2,1,2);plot(om/pi,an); xlabel('(b) Normalised frequency --.'); ylabel('Phase in radians --.'); As an example, enter the passband ripple 0.3 enter the stopband ripple 40 enter the passband freq 1500 enter the stopband freq 2000 enter the sampling freq 9000

The amplitude and phase responses of Butterworth band-pass digital filter are shown in Fig. 16.21.

% Program for the design of Butterworth Band stop digital filter

clc; close all;clear all; format long rp5input('enter the passband ripple'); rs5input('enter the stopband ripple'); wp5input('enter the passband freq'); ws5input('enter the stopband freq'); fs5input('enter the sampling freq'); w152*wp/fs;w252*ws/fs; [n]5buttord(w1,w2,rp,rs); wn5[w1 w2]; [b,a]5butter(n,wn,'stop'); w50:.01:pi; [h,om]5freqz(b,a,w); m520*log10(abs(h)); an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel('Gain in dB --.');xlabel('(a) Normalised frequency --.'); subplot(2,1,2);plot(om/pi,an); xlabel('(b) Normalised frequency --.'); ylabel('Phase in radians --.'); As an example, enter the passband ripple 0.4 enter the stopband ripple 46 enter the passband freq 1100 enter the stopband freq 2200 enter the sampling freq 6000

The amplitude and phase responses of the Butterworth bandstop digital filter are shown in Fig. 16.22. (h) The amplitude and phase responses of Chebyshev type -1 low-pass digital filter are shown in Fig. 16.23.

% Program for the design of Chebyshev Type-1 Bandpass digital filter

clc; close all;clear all; format long rp5input('enter the passband ripple...'); rs5input('enter the stopband ripple...'); wp5input('enter the passband freq...'); ws5input('enter the stopband freq...'); fs5input('enter the sampling freq...'); w152*wp/fs;w252*ws/fs; ( The amplitude and phase responses of Chebyshev type -1 bandstop digital filter are shown in Fig. 16.26. (

% Program for the design of Chebyshev Type-2 Bandpass digital filter

clc; close all;clear all; format long rp5input('enter the passband ripple...'); rs5input('enter the stopband ripple...'); wp5input('enter the passband freq...'); ws5input('enter the stopband freq...'); fs5input('enter the sampling freq...'); w152*wp/fs;w252*ws/fs; [n]5cheb2ord(w1,w2,rp,rs); wn5[w1 w2]; [b,a]5cheby2(n,rs,wn,'bandpass'); w50:.01/pi:pi; [h,om]5freqz(b,a,w); m520*log10(abs(h)); an5angle(h); subplot(2,1,1);plot(om/pi,m); ylabel('Gain in dB --.');xlabel('(a) Normalised frequency --.'); subplot(2,1,2);plot(om/pi,an); xlabel('(b) Normalised frequency --.'); ylabel('Phase in radians --.'); The amplitude and phase responses of Chebyshev type -2 bandpass digital filter are shown in Fig. 16.29.

FIR FILTER DESIGN USING WINDOW TECHNIQUES

In the design of FIR filters using any window technique, the order can be calculated using the formula given by where dp is the passband ripple, ds is the stopband ripple, f p is the passband frequency, f s is the stopband frequency and F s is the sampling frequency.

Rectangular Window

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter 5. Find the window coefficients using Eq. 7.37 6. Draw the magnitude and phase responses.

% Program for the design of FIR Low pass, High pass, Band pass and Bandstop filters using rectangular window

clc;clear all;close all; rp5input('enter the passband ripple'); rs5input('enter the stopband ripple'); fp5input('enter the passband freq'); fs5input('enter the stopband freq'); f5input('enter the sampling freq'); wp52*fp/f;ws52*fs/f; num5220*log10(sqrt(rp*rs))213; The gain responses of low-pass, high-pass, bandpass and bandstop filters using rectangular window are shown in Fig. 16.31.

Bartlett Window

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter 5. Find the filter coefficients 6. Draw the magnitude and phase responses.

% Program for the design of FIR Low pass, High pass, Band pass and Bandstop filters using Bartlett window

clc;clear all;close all; rp5input('enter the passband ripple'); rs5input('enter the stopband ripple'); fp5input('enter the passband freq'); fs5input('enter the stopband freq'); f5input('enter the sampling freq'); wp52*fp/f;ws52*fs/f; num5220*log10(sqrt(rp*rs))213; dem514.6*(fs2fp)/f; n5ceil(num/dem); n15n11; if (rem(n,2)˜50) n15n; n5n21; end y5bartlett(n1); The gain responses of low-pass, high-pass, bandpass and bandstop filters using Bartlett window are shown in Fig. 16.32.

Blackman window

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter 5. Find the window coefficients using Eq. 7.45 6. Draw the magnitude and phase responses.

% Program for the design of FIR Low pass, High pass, Band pass and Band stop digital filters using Blackman window

clc;clear all;close all; rp5input('enter the passband ripple'); rs5input('enter the stopband ripple'); fp5input('enter the passband freq'); fs5input('enter the stopband freq'); f5input('enter the sampling freq'); wp52*fp/f;ws52*fs/f; num5220*log10(sqrt(rp*rs))213; dem514.6*(fs2fp)/f; n5ceil(num/dem); n15n11; if (rem(n,2)˜50) n15n; n5n21; end y5blackman(n1); The gain responses of low-pass, high-pass, bandpass and bandstop filters using Blackman window are shown in Fig. 16.33.

Chebyshev Window

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter 5. Find the filter coefficients 6. Draw the magnitude and phase responses.

% Program for the design of FIR Lowpass, High pass, Band pass and Bandstop filters using Chebyshev window

clc;clear all;close all; rp5input('enter the passband ripple'); rs5input('enter the stopband ripple'); fp5input('enter the passband freq'); fs5input('enter the stopband freq'); f5input('enter the sampling freq'); Gain in dB Gain in dB r5input('enter the ripple value(in dBs)'); wp52*fp/f;ws52*fs/f; num5220*log10(sqrt(rp*rs))213; dem514.6*(fs2fp)/f; n5ceil(num/dem); if(rem(n,2)˜50) n5n11; end y5chebwin(n,r); The gain responses of low-pass, high-pass, bandpass and bandstop filters using Chebyshev window are shown in Fig. 16.34.

Hamming Window

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter 5. Find the window coefficients using Eq. 7.40 6. Draw the magnitude and phase responses.

% Program for the design of FIR Low pass, High pass, Band pass and Bandstop filters using Hamming window

clc;clear all;close all; rp5input('enter the passband ripple'); rs5input('enter the stopband ripple'); fp5input('enter the passband freq'); fs5input('enter the stopband freq'); f5input('enter the sampling freq'); wp52*fp/f;ws52*fs/f; num5220*log10(sqrt(rp*rs))213; dem514.6*(fs2fp)/f; n5ceil(num/dem); n15n11; if (rem(n,2)˜50) n15n; n5n21; end y5hamming(n1); % low-pass filter b5fir1(n,wp,y); [h,o]5freqz(b,1,256); m520*log10(abs(h)); subplot(2,2,1);plot(o/pi,m);ylabel('Gain in dB --.'); xlabel('(a) Normalised frequency --.'); % high-pass filter b5fir1(n,wp,'high',y); [h,o]5freqz(b,1,256); m520*log10(abs(h)); subplot(2,2,2);plot(o/pi,m);ylabel('Gain in dB --.'); xlabel('(b) Normalised frequency --.'); % band pass filter wn5[wp ws]; b5fir1(n,wn,y); [h,o]5freqz(b,1,256); m520*log10(abs(h)); subplot(2,2,3);plot(o/pi,m);ylabel('Gain in dB --.'); xlabel('(c) Normalised frequency --.'); % band stop filter b5fir1(n,wn,'stop',y); [h,o]5freqz(b,1,256); m520*log10(abs(h)); subplot(2,2,4);plot(o/pi,m);ylabel('Gain in dB --.'); xlabel('(d) Normalised frequency --.'); As an example, enter the passband ripple 0.02 enter the stopband ripple 0.01 enter the passband freq 1200 enter the stopband freq 1700 enter the sampling freq 9000 The gain responses of low-pass, high-pass, bandpass and bandstop filters using Hamming window are shown in Fig. 16.35.

Hanning Window

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter 5. Find the window coefficients using Eq. 7.44 6. Draw the magnitude and phase responses.

% Program for the design of FIR Low pass, High pass, Band pass and Band stop filters using Hanning window

clc;clear all;close all; rp5input('enter the passband ripple'); rs5input('enter the stopband ripple'); fp5input('enter the passband freq'); fs5input('enter the stopband freq'); f5input('enter the sampling freq'); wp52*fp/f;ws52*fs/f; num5220*log10(sqrt(rp*rs))213; dem514.6*(fs2fp)/f; n5ceil(num/dem); n15n11; if (rem(n,2)˜50) n15n; n5n21; end y5hamming(n1); The gain responses of low-pass, high-pass, bandpass and bandstop filters using Hanning window are shown in Fig. 16.36.

Kaiser Window

Algorithm 1. Get the passband and stopband ripples 2. Get the passband and stopband edge frequencies 3. Get the sampling frequency 4. Calculate the order of the filter 5. Find the window coefficients using Eqs 7.46 and 7.47 6. Draw the magnitude and phase responses.

% Program for the design of FIR Low pass, High pass, Band pass and Bandstop filters using Kaiser window

clc;clear all;close all; rp5input('enter the passband ripple'); rs5input('enter the stopband ripple'); fp5input('enter the passband freq'); fs5input('enter the stopband freq'); f5input('enter the sampling freq'); beta5input('enter the beta value'); wp52*fp/f;ws52*fs/f; num5220*log10(sqrt(rp*rs))213; dem514.6*(fs2fp)/f; n5ceil(num/dem); n15n11; if (rem(n,2)˜50) n15n; n5n21; end y5kaiser(n1,beta); % low-pass filter b5fir1(n,wp,y); [h,o]5freqz(b,1,256); m520*log10(abs(h)); subplot(2,2,1);plot(o/pi,m);ylabel('Gain in dB -->'); xlabel('(a) Normalised frequency -->'); % high-pass filter b5fir1(n,wp,'high',y); [h,o]5freqz(b,1,256); m520*log10(abs(h)); subplot(2,2,2);plot(o/pi,m);ylabel('Gain in dB -->'); xlabel('(b) Normalised frequency -->'); % band pass filter wn5[wp ws]; b5fir1(n,wn,y); [h,o]5freqz(b,1,256); m520*log10(abs(h)); subplot(2,2,3);plot(o/pi,m);ylabel('Gain in dB -->'); xlabel('(c) Normalised frequency -->'); % band stop filter b5fir1(n,wn,'stop',y); [h,o]5freqz(b,1,256); m520*log10(abs(h)); subplot(2,2,4);plot(o/pi,m);ylabel('Gain in dB -->'); xlabel('(d) Normalised frequency -->'); As an example, enter the passband ripple 0.02 enter the stopband ripple 0.01 enter the passband freq 1000 enter the stopband freq 1500 enter the sampling freq 10000 enter the beta value 5.8

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

UPSAMPLING A SINUSOIDAL SIGNAL

% Program for upsampling a sinusoidal signal by factor L

N5input('Input length of the sinusoidal sequence5'); L5input('Up Samping factor5'); fi5input('Input signal frequency5'); % Generate the sinusoidal sequence for the specified length N n50:N21; x5sin(2*pi*fi*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'); As an example, enter length of input sentence … 25 enter upsampling factor … 3 enter the value of a … 0.95

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

DOWN SAMPLING AN EXPONENTIAL SEQUENCE

% Program for downsampling an exponential sequence by a factor M

N5input('enter the length of the output sequence …'); M5input('enter the down sampling factor …'); As an example, enter the length of the output sentence … 25 enter the downsampling factor … 3 enter the value of a … 0.95

The input and output sequences of downsampling an exponential sequence a n are shown in Fig. 16.39.

ESTIMATION OF POWER SPECTRAL DENSITY (PSD)

% Program for estimating PSD of two sinusoids plus noise % Algorithm; % 1:Get the frequencies of the two sinusoidal waves % 2:Get the sampling frequency % 3:Get the length of the sequence to be considered % 4: Get the two FFT lengths for comparing the corresponding power spectral densities clc; close all; clear all; f15input('Enter the frequency of first sinusoid'); f25input('Enter the frequency of second sinusoid'); fs5input('Enter the sampling frequency'); N5input("Enter the length of the input sequence'); N15input("Enter the input FFT length 1'); N25input("Enter the input FFT length 2'); %Generation of input sequence t50:1/fs:1; x52*sin(2*pi*f1*1)13*sin(2*pi*f2*t)2randn(size(t));

%Generation of psd for two different FFT lengths

Pxx15abs(fft(x,N1)).^2/(N11); Pxx25abs(fft(x,N2)).^2/(N11); %Plot the psd; subplot(2,1,1); plot ((0:(N121))/N1*fs,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*fs,10*log10(Pxx2)); xlabel('Frequency in Hz'); ylabel('Power spectrum in dB'); title('[PSD with FFT length,num2str(N2)]');

Pxx15(abs(fft(x(1:256))).^21abs(fft(x(257:512))).^21 abs(fft(x(513:768))).^2/(256*3); %using nonoverlapping sections Pxx25(abs(fft(x(1:256))).^21abs(fft(x(129:384))).^21ab s(fft(x(257:512))).^21abs(fft(x(385:640))).^21abs(fft( x(513:768))).^21abs(fft(x(641:896))).^2/(256*6); %using overlapping sections % Plot the psd; subplot (2,1,1); plot ((0:255)/256*fs,10*log10 (Pxx1)

PSD ESTIMATOR

% Program for estimating PSD of a two sinusoids plus noise using

%(i)non-overlapping sections %(ii)overlapping sections and averaging the periodograms clc; close all; clear all; f15input('Enter the frequency of first sinusoid'); f25input('Enter the frequency of second sinusoid'); fs5input('Enter the sampling frequency'); N5input("Enter the length of the input sequence'); N15input("Enter the input FFT length 1'); N25input("Enter the input FFT length 2');

%Generation of input sequence

t50:1/fs:1; x52*sin(2*pi*f1*1)13*sin(2*pi*f2*t)2randn(size(t));

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 for the input signal of two sinusoids plus noise using Hanning window f15input('Enter the frequency of first sinusoid'); f25input('Enter the frequency of second sinusoid'); fs5input('Enter the sampling frequency'); t50:1/fs:1; w5hanning(256); x52*sin(2*pi*f1*t)13*sin(2*pi*f2*t)2randn(size(t)); Pxx5(abs(fft(w.*x(1:256))).^21abs(fft(w.*x(129:384))).2 1abs(fft(w.*x(257:512))).^21abs(fft(w.*x(385:640))).^2 1abs(fft(w.*x(513:768))).^21abs(fft(w.*x(641:896))).^2/ (norm(w)^2*6); Plot((0:255)/256*fs,10*log10(Pxx));

WELCH PSD ESTIMATOR

GROUP DELAY

% Program for computing group delay of a rational transfer function on a given frequency interval

function D5grpdly(b,a,K,theta); b5input ('enter the numerator polynomials5'); a5input ('enter the denominator polynomials5'); K5input ('enter the number of frequency response points5'); theta5input ('enter the theta value5'); a5reshape(a,1,length(a)); b5reshape(b,1,length(b)); if (length(a)551)%case of FIR bd52j*(0:length(b)21).*b; if(nargin553), B5frqresp(b,1,K); Bd5frqresp(bd,1,K); else, B5frqresp(b,1,K,theta); Bd5frqresp(bd,1,K,theta); end D5(real(Bd).*imag(B)2real(B).*imag(Bd))./abs(B

IIR FILTER DESIGN-BILINEAR TRANSFORMATION

% Program for transforming an analog filter into a digial filter using bilinear transformation

function [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 of the analog filter5'); 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

MULTIBAND FIR FILTER DESIGN

% Program for the design of multiband FIR filters

function h5firdes(N,spec,win); N5input('enter the length of the filter5'); spec5input('enter the low,high cutoff frequencies and gain5'); win5input('enter the window length5'); flag5rem(N,2); [K,m]5size(spec); n5(0:N)2N/2; if (˜flag),n(N/211)51; end,h5zeros(1,N11); for k51:K temp5(spec (k,3)/pi)*(sin(spec(k,2)*n)2sin(spec(k,1) *n))./n; if (˜flag);temp(N/211)5spec(k,3)*(spec(k,2)2 spec(k,1))/pi; end h5h1temp; end if (nargin553), h5h.*reshape(win,1,N11); end

ANALYSIS FILTER BANK

% Program for maximally decimated uniform DFT analysis filter bank

function u5dftanal(x,g,M); g5input('enter the filter coefficient5'); x5input('enter the input sequence5'); M5input('enter the decimation factor5'); 1g5length(g); 1p5floor( (1g21) Fig. 16.44 using LMS algorithm clc; close all; clear all; format short T5input('Enter the symbol interval'); br5input('Enter the bit rate value'); rf5input('Enter the roll off factor'); n5[210 10]; y55000*rcosfir(rf,n,br,T); %Transmit filter pulse shape is assumed as raised cosine ds5[5 2 5 2 5 2 5 2 5 5 5 5 2 2 2 5 5 5 5]; % data sequence lam5max(eig(r)); la5min(eig(r)); l5lam/la; w5inv(r)*p; % Initial value of filter coefficients e5rs2filter(w,l,ds); s51; mu51.5/lam; ni51; while (s 1 e210) w15w22*mu*(e.*ds)' ; % LMS algorithm adaptation rs y45filter(w1,1,ds); % Estimated echo signal using LMS algorithm e5y42rs; s50; e15xcorr(e); for i51:length(e1), s5s1e1(i); end s5s/length(e1); if (y455rs) break end ni5ni11; w5w1; end figure(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(fs); title('desired signal'); subplot(2,2,4); plot(fs1); title('echo added with desired signal'); figure(2); subplot(2,1,1); plot(y4); title('estimated echo signal using LMS algorithm'); subplot(2,1,2); plot(fs12y4); title('echo cancelled signal'); (c) Repeat parts (a) and (b) using the Hamming window (d) Repeat parts (a) and (b) using the Bartlett window. 16.21 A digital low-pass filter is required to meet the following specfications Passband ripple # 1 dB Passband edge 4 kHz Stopband attenuation  40 dB Stopband edge 6 kHz Sample rate 24 kHz The filter is to be designed by performing a bilinear transformation on an analog system function. Determine what order Butterworth, Chebyshev and elliptic analog design must be used to meet the specifications in the digital implementation. 16.22 An IIR digital low-pass filter is required to meet the following specfications 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 formulas to determine the filter order for (a) Digital Butterworth filter (b) Digital Chebyshev filter (c) Digital elliptic filter 16.23 An analog signal of the form x a (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. The signal x(n) is applied to a device (decimator) which reduces the rate by a factor of two. Determine the output spectrum. (c) Show that the spectrum is simply the Fourier transform of x(2n).

CANCELLATION OF ECHO PRODUCED ON THE TELEPHONE-PASS BAND CHANNEL

16.25 Design a digital type-I Chebyshev low-pass filter operating at a sampling rate of 44.1kHz with a passband frequency at 2kHz, a pass band ripple of 0.4dB, and a minimum stopband attenuation of 50dB at 12kHz using the impulse invariance method and the bilinear transformation method. Determine the order of analog filter prototype and design the analog prototype filter. Plot the gain and phase responses of the both designs using 16.26 Design a linear phase FIR high-pass filter with following specifications:

Stopband edge at 0.5p, passband edge at 0.7p, maximum passband attenuation of 0.15dB and a minimum stopband attenuation of 40dB. Use each of the following windows for the design. Hamming, Hanning, Blackman and Kaiser. Show the impulse response coefficients and plot the gain response of the designed filters for each case.

16.27 Design using the windowed Fourier series approach a linear phase FIR lowpass filter with the following specifications: pass band edge at 1 rad/s, stop band edge at 2rad/s, maximum passband attenuation of 0.2dB, minimum stopband attenuation of 50dB and a sampling frequency of 10rad/s. Use each of the following windows for the design: Hamming, Hanning, Blackman, Kaiser and Chebyshev. Show the impulse response coefficients and plot the gain response of designed filters for each case.

16.28 Design a two-channel crossover FIR low-pass and high-pass filter pair for digital audio applications. The low-pass and high-pass filters are of length 31 and have a crossover frequency of 2kHz operating at a sampling rate of 44.1KHz. Use the function 'fir1' with a Hamming window to design the lowpass filter while the high-pass filter is derived from the low-pass filter using the delay complementary property. Plot the gain responses of both filters on the same figure. What is the minimum number of delays and multipliers needed to implement the crossover network?

16.29 Design a digital network butterworth low-pass filter operating at sampling rate of 44.1kHz with a 0.5dB cutoff frequency at 2kHz and a minimum stopband attenuation of 45dB at 10kHz using the impulse invariance method and the bilinear transformation method. Assume the sampling interval for the impulse invariance design to be equal to 1. Determine the order of the analog filter prototype and then design the analog prototype filter. Plot the gain and phase responses of both designs. Compare the performances of the filters. Show all steps used in the design. Does the sampling interval have any effect on the design of the digital filter design based on the impulse invariance method? Hint The order of filter is