Download Matlab programs in digital signal problem and more Exercises Digital Signal Processing in PDF only on Docsity!
MATLAB Programs
Chapter 16
16.1 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.
16.2 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; t 52 2:1:2; y 5 [zeros(1,2),ones(1,1),zeros(1,2)];subplot(2,2,1);stem(t,y);
816 Digital Signal Processing
ylabel(‘Amplitude --.’); xlabel(‘(a) n --.’);
% Program for the generation of unit step sequence [u(n) 2 u(n 2 N]
n 5 input(‘enter the N value’); t 5 0:1:n 2 1; y1 5 ones(1,n);subplot(2,2,2); stem(t,y1);ylabel(‘Amplitude --.’); xlabel(‘(b) n --.’);
% Program for the generation of ramp sequence
n1 5 input(‘enter the length of ramp sequence’); t 5 0:n1; subplot(2,2,3);stem(t,t);ylabel(‘Amplitude --.’); xlabel(‘(c) n --.’);
% Program for the generation of exponential sequence
n2 5 input(‘enter the length of exponential sequence’); t 5 0:n2; a 5 input(‘Enter the ‘a’ value’); y2 5 exp(a*t);subplot(2,2,4); stem(t,y2);ylabel(‘Amplitude --.’); xlabel(‘(d) n --.’);
% Program for the generation of sine sequence
t 5 0:.01:pi; y 5 sin(2pit);figure(2); subplot(2,1,1);plot(t,y);ylabel(‘Amplitude --.’); xlabel(‘(a) n --.’);
% Program for the generation of cosine sequence
t 5 0:.01:pi; y 5 cos(2pit); 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.
818 Digital Signal Processing
16.3 DISCRETE CONVOLUTION
16.3.1 Linear Convolution
Algorithm
- Get two signals x ( m )and h ( p )in matrix form
- The convolved signal is denoted as y ( n )
- y ( n )is given by the formula
y ( n ) 5 [ ( ) x k^^ h n (^^ k )] k
− =−∞
∞ ∑ where^ n^5 0 to^ m^^1 p^^2
- Stop
% Program for linear convolution of the sequence x 5 [1, 2] and h 5 [1, 2, 4]
clc; clear all; close all; x 5 input(‘enter the 1st sequence’); h 5 input(‘enter the 2nd sequence’); y 5 conv(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 y 5 1 4 8 8 Figure 16.2 shows the discrete input signals x ( n )and h ( n )and the convolved output signal y ( n ).
Fig. 16.2 ( Contd. )
Amplitude
Amplitude
n
n
0
0
1
1
1
2
(a)
(b)
2
3
1
2
3
2
4
MATLAB Programs 819
16.3.2 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:n 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:n y(j)=M(i,j)+y(j); end end disp(‘The output sequence is y(n)= ‘); disp(y);
Fig. 16.2 Discrete Linear Convolution
Amplitude
Amplitude
Amplitude
n
n
n
0
0
0
1
2
1
1
1
2
2
(a)
(b)
(c)
3
2
3
4
1
2
4
3
6
2
4
8
Amplitude
Amplitude
Amplitude
n
n
n
0
0
0
1
2
1
1
1
2
2
(a)
(b)
(c)
3
2
3
4
1
2
4
3
6
2
4
8
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=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:n 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)
822 Digital Signal Processing
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
- Get two signals x ( m )and h ( p )in matrix form
- The correlated signal is denoted as y ( n )
- y ( n )is given by the formula
y ( n ) 5 [ ( ) x k h k ( n )] k
=−∞
∞ ∑
where n 52 [max ( m , p ) 2 1] to [max ( m , p ) 2 1]
- Stop
% Program for computing cross-correlation of the sequences x 5 [1, 2, 3, 4] and h 5 [4, 3, 2, 1]
clc; clear all; close all; x 5 input(‘enter the 1st sequence’); h 5 input(‘enter the 2nd sequence’); y 5 xcorr(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 --.’);
824 Digital Signal Processing
% Program for computing autocorrelation function
x 5 input(‘enter the sequence’); y 5 xcorr(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 y 5 4 11 20 ↑
Figure 16.4 shows the discrete input signal x ( n )and its auto-correlated output signal y ( n ).
Amplitude
Amplitude
n
n
0
1
1
1
2
2
3
3 (a)
5
4 (b) y (n)
6
4
7
2
3
4
0
5
10
15
20
25
30
Fig. 16.4 Discrete Auto-correlation
16.5 STABILITY TEST
% Program for stability test
clc;clear all;close all; b 5 input(‘enter the denominator coefficients of the filter’); k 5 poly2rc(b); knew 5 fliplr(k); s 5 all(abs(knew)1); if(s 55 1) disp(‘“Stable system”’);
MATLAB Programs 825
else disp(‘“Non-stable system”’); end As an example, enter the denominator coefficients of the filter [1 2 1 .5] “Stable system”
16.6 SAMPLING THEOREM
The sampling theorem can be understood well with the following example.
Example 16.1 Frequency analysis of the amplitude modulated discrete-time signal x ( n ) 5 cos 2 p f 1 n 1 cos 2p f 2 n
where f 1
= (^) and f 2
= (^) modulates the amplitude-modulated signal is
x (^) c ( n )^5 cos 2p^ fc n where f (^) c 5 50/128. The resulting amplitude-modulated signal is x (^) am ( n ) 5 x ( n ) cos 2p fc n Using MATLAB program , (a) sketch the signals x ( n ), xc ( n ) and xam ( n ), 0 # n # 255 (b) compute and sketch the 128-point DFT of the signal xam ( n ), 0 # n # 127 (c) compute and sketch the 128-point DFT of the signal xam ( n ), 0 # n # 99
Solution
% Program
Solution for Section (a) clc;close all;clear all; f1 5 1/128;f2 5 5/128;n 5 0:255;fc 5 50/128; x 5 cos(2pif1n) 1 cos(2pif2n); xa 5 cos(2pifcn); xamp 5 x.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 computation 2 solution for Section (b) n 5 0:127;figure;n1 5 128; f1 5 1/128;f2 5 5/128;fc 5 50/128; x 5 cos(2pif1n) 1 cos(2pif2n); xc 5 cos(2pifcn); xa 5 cos(2pifcn); ( Contd. )
MATLAB Programs 827
20
Amplitude
40 60 80 100 120 n
0 140
− 10
− 5
0
5
10
15
20
25
Fig. 16.5(b) 128-point DFT of the Signal xam ( n ) , 0 # n # 127
20 40 60 80 100 120 140 n
0
− 5
0
5
10
15
20
25
30
35
Amplitude
Fig. 16.5(c) 128-point DFT of the Signal xam ( n ) , 0 # n # 99
828 Digital Signal Processing
xamp 5 x.xa;xam 5 fft(xamp,n1); stem(n,xam);title(‘xamp(n)’);xlabel(‘n --.’); ylabel(‘amplitude’); %128 point DFT computation 2 solution for Section (c) n 5 0:99;figure;n2 5 0:n1 2 1; f1 5 1/128;f2 5 5/128;fc 5 50/128; x 5 cos(2pif1n) 1 cos(2pif2n); xc 5 cos(2pifcn); xa 5 cos(2pifcn); xamp 5 x.xa; for i 5 1:100, xamp1(i) 5 xamp(i); end xam 5 fft(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 xa ( n ) and amplitude modulated signal xam ( n ) are shown in Fig. 16.5(a). Fig. 16.5 (b) shows the 128-point DFT of the signal xam ( n ) for 0 # n # 127 and Fig. 16.5 (c) shows the 128-point DFT of the signal xam ( n ), 0 # n # 99.
16.7 FAST FOURIER TRANSFORM
Algorithm
- Get the signal x ( n )of length N in matrix form
- Get the N value
- The transformed signal is denoted as
x k x n e k N
j N nk n
N ( ) = ( ) for ≤ ≤ −
−
− ∑
2
0
1 0 1
p
_ **\% Program for computing discrete Fourier transform_**
clc;close all;clear all; x 5 input(‘enter the sequence’); n 5 input(‘enter the length of fft’); X(k) 5 fft(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 of fft 8 X(k) 5 Columns 1 through 4 28.0000 2 4.0000 1 9.6569i 2 4.0000 1 4.0000i 2 4. 1 1.6569i Columns 5 through 8 2 4.0000 2 4.0000 2 1.6569i 2 4.0000 2 4.0000i 2 4. 2 9.6569i
830 Digital Signal Processing
w1 5 2wp/fs;w2 5 2ws/fs; [n,wn] 5 buttord(w1,w2,rp,rs,’s’); [z,p,k] 5 butter(n,wn); [b,a] 5 zp2tf(z,p,k); [b,a] 5 butter(n,wn,’s’); w 5 0:.01:pi; [h,om] 5 freqs(b,a,w); m 5 20*log10(abs(h)); an 5 angle(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. 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.
Fig. 16.7 Butterworth Low-pass Analog Filter ( a ) Amplitude Response and ( b ) Phase Response
− 250
− 200
− 150
− 4
− 2
2
4
0
− 100
− 50
50
Gain in dB
0
(a)
(b)
Normalised frequency
Normalised frequency
1
1
0
0
Phase in radians
− 250
− 200
− 150
− 4
− 2
2
4
0
− 100
− 50
50
Gain in dB
0
(a)
(b)
Normalised frequency
Normalised frequency
1
1
0
0
Phase in radians
MATLAB Programs 831
16.8.2 High-pass Filter
Algorithm
- Get the passband and stopband ripples
- Get the passband and stopband edge frequencies
- Get the sampling frequency
- Calculate the order of the filter using Eq. 8.
- Find the filter coefficients
- Draw the magnitude and phase responses.
% Program for the design of Butterworth analog high—pass filter
clc;
close all;clear all;
format long
rp 5 input(‘enter the passband ripple’);
rs 5 input(‘enter the stopband ripple’);
wp 5 input(‘enter the passband freq’);
ws 5 input(‘enter the stopband freq’);
fs 5 input(‘enter the sampling freq’);
w1 5 2wp/fs;w2 5 2ws/fs;
[n,wn] 5 buttord(w1,w2,rp,rs,’s’);
[b,a] 5 butter(n,wn,’high’,’s’);
w 5 0:.01:pi;
[h,om] 5 freqs(b,a,w);
m 5 20*log10(abs(h));
an 5 angle(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.
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.
MATLAB Programs 833
[n] 5 buttord(w1,w2,rp,rs); wn 5 [w1 w2]; [b,a] 5 butter(n,wn,’bandpass’,’s’); w 5 0:.01:pi; [h,om] 5 freqs(b,a,w); m 5 20*log10(abs(h)); an 5 angle(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. enter the stopband ripple... 36 enter the passband freq... 1500 enter the stopband freq... 2000 enter the sampling freq... 6000
The amplitude and phase responses of Butterworth bandpass analog filter are shown in Fig. 16.9.
Fig. 16.9 Butterworth Bandpass Analog Filter ( a ) Amplitude Response and ( b ) Phase Response
Gain in dB
Phase in radians
(a)
(b)
− 1000
− 4
4
− 2
2
0
− 800
− 600
− 200 − 400
0
200
Normalised frequency
Normalised frequency
1
1
0
0
834 Digital Signal Processing
16.8.4 Bandstop Filter
Algorithm
- Get the passband and stopband ripples
- Get the passband and stopband edge frequencies
- Get the sampling frequency
- Calculate the order of the filter using Eq. 8.
- Find the filter coefficients
- Draw the magnitude and phase responses.
% Program for the design of Butterworth analog Bandstop filter
clc;
close all;clear all; format long
rp 5 input(‘enter the passband ripple...’); rs 5 input(‘enter the stopband ripple...’);
wp 5 input(‘enter the passband freq...’); ws 5 input(‘enter the stopband freq...’);
fs 5 input(‘enter the sampling freq...’); w1 5 2wp/fs;w2 5 2ws/fs;
[n] 5 buttord(w1,w2,rp,rs,’s’); wn 5 [w1 w2];
[b,a] 5 butter(n,wn,’stop’,’s’); w 5 0:.01:pi;
[h,om] 5 freqs(b,a,w); m 5 20*log10(abs(h));
an 5 angle(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. 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.
