Digital Signal Processing

Scilab Textbook Companion for
Digital Signal Processing

by R. Babu
1
Created by
Mohammad Faisal Siddiqui
B.Tech (pursuing)
Electronics Engineering
Jamia Milia Islamia
College Teacher
Dr. Sajad A. Loan, JMI, New D
Cross-Checked by
Santosh Kumar, IITB
August 10, 2013
1
Funded by a grant from the National Mission on Education through ICT,
http://spoken-tutorial.org/NMEICT-Intro. This Textbook Companion and Scilab
codes written in it can be downloaded from the Textbook Companion Project
section at the website http://scilab.in
Book Description
Title: Digital Signal Processing
Author: R. Babu
Publisher: Scitech Publications
Edition: 4
Year: 2010
ISBN: 81-8371-081-7
1
Scilab numbering policy used in this document and the relation to the
above book.
Exa Example (Solved example)
Eqn Equation (Particular equation of the above book)
AP Appendix to Example(Scilab Code that is an Appednix to a particular
Example of the above book)
For example, Exa 3.51 means solved example 3.51 of this book. Sec 2.3 means
a scilab code whose theory is explained in Section 2.3 of the book.
2
Contents
List of Scilab Codes 4
1 DISCRETE TIME SIGNALS AND LINEAR SYSTEMS 11
2 THE Z TRANSFORM 34
3 THE DISCRETE FOURIER TRANSFORM 59
4 THE FAST FOURIER TRANSFORM 79
5 INFINITE IMPULSE RESPONSE FILTERS 90
6 FINITE IMPULSE RESPONSE FILTERS 102
7 FINITE WORD LENGTH EFFECTS IN DIGITAL FIL-
TERS 139
8 MULTIRATE SIGNAL PROCESSING 141
9 STATISTICAL DIGITAL SIGNAL PROCESSING 143
11 DIGITAL SIGNAL PROCESSORS 146
3
List of Scilab Codes
Exa 1.1 Continuous Time Plot and Discrete Time Plot . . . . 11
Exa 1.2 Continuous Time Plot and Discrete Time Plot . . . . 13
Exa 1.3.a Evaluate the Summations . . . . . . . . . . . . . . . . 14
Exa 1.3.b Evaluate the Summations . . . . . . . . . . . . . . . . 14
Exa 1.4.a Check for Energy or Power Signals . . . . . . . . . . . 15
Exa 1.4.d Check for Energy or Power Signals . . . . . . . . . . . 15
Exa 1.5.a Determining Periodicity of Signal . . . . . . . . . . . . 16
Exa 1.5.c Determining Periodicity of Signal . . . . . . . . . . . . 18
Exa 1.5.d Determining Periodicity of Signal . . . . . . . . . . . . 19
Exa 1.11 Stability of the System . . . . . . . . . . . . . . . . . 21
Exa 1.12 Convolution Sum of Two Sequences . . . . . . . . . . 21
Exa 1.13 Convolution of Two Signals . . . . . . . . . . . . . . . 22
Exa 1.18 Cross Correlation of Two Sequences . . . . . . . . . . 22
Exa 1.19 Determination of Input Sequence . . . . . . . . . . . . 23
Exa 1.32.a Plot Magnitude and Phase Response . . . . . . . . . . 24
Exa 1.37 Sketch Magnitude and Phase Response . . . . . . . . . 24
Exa 1.38 Plot Magnitude and Phase Response . . . . . . . . . . 26
Exa 1.45 Filter to Eliminate High Frequency Component . . . . 27
Exa 1.57.a Discrete Convolution of Sequences . . . . . . . . . . . 29
Exa 1.61 Fourier Transform . . . . . . . . . . . . . . . . . . . . 29
Exa 1.62 Fourier Transform . . . . . . . . . . . . . . . . . . . . 30
Exa 1.64.a Frequency Response of LTI System . . . . . . . . . . . 31
Exa 1.64.c Frequency Response of LTI System . . . . . . . . . . . 31
Exa 2.1 z Transform and ROC of Causal Sequence . . . . . . . 34
Exa 2.2 z Transform and ROC of Anticausal Sequence . . . . . 34
Exa 2.3 z Transform of the Sequence . . . . . . . . . . . . . . 35
Exa 2.4 z Transform and ROC of the Signal . . . . . . . . . . 36
Exa 2.5 z Transform and ROC of the Signal . . . . . . . . . . 36
4
Exa 2.6 Stability of the System . . . . . . . . . . . . . . . . . 37
Exa 2.7 z Transform of the Signal . . . . . . . . . . . . . . . . 37
Exa 2.8.a z Transform of the Signal . . . . . . . . . . . . . . . . 37
Exa 2.10 z Transform Computation . . . . . . . . . . . . . . . . 38
Exa 2.13.a z Transform of Discrete Time Signals . . . . . . . . . . 39
Exa 2.13.b z Transform of Discrete Time Signals . . . . . . . . . . 40
Exa 2.13.c z Transform of Discrete Time Signals . . . . . . . . . . 40
Exa 2.13.d z Transform of Discrete Time Signals . . . . . . . . . . 41
Exa 2.16 Impulse Response of the System . . . . . . . . . . . . 41
Exa 2.17 Pole Zero Plot of the Dierence Equation . . . . . . . 42
Exa 2.19 Frequency Response of the System . . . . . . . . . . . 44
Exa 2.20.a Inverse z Transform Computation . . . . . . . . . . . 44
Exa 2.22 Inverse z Transform Computation . . . . . . . . . . . 45
Exa 2.23 Causal Sequence Determination . . . . . . . . . . . . . 45
Exa 2.35.a Pole Zero Plot of the System . . . . . . . . . . . . . . 47
Exa 2.35.b Unit Sample Response of the System . . . . . . . . . . 49
Exa 2.38 Determine Output Response . . . . . . . . . . . . . . 50
Exa 2.40 Input Sequence Computation . . . . . . . . . . . . . . 52
Exa 2.41.a z Transform of the Signal . . . . . . . . . . . . . . . . 53
Exa 2.41.b z Transform of the Signal . . . . . . . . . . . . . . . . 53
Exa 2.41.c z Transform of the Signal . . . . . . . . . . . . . . . . 54
Exa 2.45 Pole Zero Pattern of the System . . . . . . . . . . . . 55
Exa 2.53.a z Transform of the Sequence . . . . . . . . . . . . . . 55
Exa 2.53.b z Transform of the Signal . . . . . . . . . . . . . . . . 55
Exa 2.53.c z Transform of the Signal . . . . . . . . . . . . . . . . 56
Exa 2.53.d z Transform of the Signal . . . . . . . . . . . . . . . . 56
Exa 2.54 z Transform of Cosine Signal . . . . . . . . . . . . . . 56
Exa 3.1 DFT and IDFT . . . . . . . . . . . . . . . . . . . . . 59
Exa 3.2 DFT of the Sequence . . . . . . . . . . . . . . . . . . 60
Exa 3.3 8 Point DFT . . . . . . . . . . . . . . . . . . . . . . . 62
Exa 3.4 IDFT of the given Sequence . . . . . . . . . . . . . . . 62
Exa 3.7 Plot the Sequence . . . . . . . . . . . . . . . . . . . . 63
Exa 3.9 Remaining Samples . . . . . . . . . . . . . . . . . . . 64
Exa 3.11 DFT Computation . . . . . . . . . . . . . . . . . . . . 65
5
Exa 3.13 Circular Convolution . . . . . . . . . . . . . . . . . . . 65
Exa 3.15 Determine Sequence x3 . . . . . . . . . . . . . . . . . 67
Exa 3.18 Output Response . . . . . . . . . . . . . . . . . . . . . 69
Exa 3.20 Output Response . . . . . . . . . . . . . . . . . . . . . 70
Exa 3.21 Linear Convolution . . . . . . . . . . . . . . . . . . . . 70
Exa 3.23.a N Point DFT Computation . . . . . . . . . . . . . . . 71
Exa 3.23.b N Point DFT Computation . . . . . . . . . . . . . . . 71
Exa 3.23.c N Point DFT Computation . . . . . . . . . . . . . . . 72
Exa 3.23.d N Point DFT Computation . . . . . . . . . . . . . . . 72
Exa 3.23.e N Point DFT Computation . . . . . . . . . . . . . . . 72
Exa 3.23.f N Point DFT Computation . . . . . . . . . . . . . . . 73
Exa 3.24 DFT of the Sequence . . . . . . . . . . . . . . . . . . 73
Exa 3.25 8 Point Circular Convolution . . . . . . . . . . . . . . 73
Exa 3.26 Linear Convolution using DFT . . . . . . . . . . . . . 74
Exa 3.27.a Circular Convolution Computation . . . . . . . . . . . 75
Exa 3.27.b Circular Convolution Computation . . . . . . . . . . . 75
Exa 3.30 Calculate value of N . . . . . . . . . . . . . . . . . . . 76
Exa 3.32 Sketch Sequence . . . . . . . . . . . . . . . . . . . . . 76
Exa 3.36 Determine IDFT . . . . . . . . . . . . . . . . . . . . . 78
Exa 4.3 Shortest Sequence N Computation . . . . . . . . . . . 79
Exa 4.4 Twiddle Factor Exponents Calculation . . . . . . . . . 80
Exa 4.6 DFT using DIT Algorithm . . . . . . . . . . . . . . . 80
Exa 4.8 DFT using DIF Algorithm . . . . . . . . . . . . . . . 81
Exa 4.9 8 Point DFT of the Sequence . . . . . . . . . . . . . . 81
Exa 4.11 IDFT of the Sequence using DIT Algorithm . . . . . . 82
Exa 4.14 DFT using DIT Algorithm . . . . . . . . . . . . . . . 83
Exa 4.15 DFT using DIF Algorithm . . . . . . . . . . . . . . . 84
Exa 4.16.a 8 Point DFT using DIT FFT . . . . . . . . . . . . . . 84
Exa 4.16.b 8 Point DFT using DIT FFT . . . . . . . . . . . . . . 85
Exa 4.17 IDFT using DIF Algorithm . . . . . . . . . . . . . . . 85
Exa 4.18 IDFT using DIT Algorithm . . . . . . . . . . . . . . . 86
Exa 4.19 FFT Computation of the Sequence . . . . . . . . . . . 86
6
Exa 4.20 8 Point DFT by Radix 2 DIT FFT . . . . . . . . . . . 86
Exa 4.21 DFT using DIT FFT Algorithm . . . . . . . . . . . . 87
Exa 4.22 Compute X using DIT FFT . . . . . . . . . . . . . . . 87
Exa 4.23 DFT using DIF FFT Algorithm . . . . . . . . . . . . 88
Exa 5.1 Order of the Filter Determination . . . . . . . . . . . 90
Exa 5.2 Order of Low Pass Butterworth Filter . . . . . . . . . 90
Exa 5.4 Analog Butterworth Filter Design . . . . . . . . . . . 91
Exa 5.5 Analog Butterworth Filter Design . . . . . . . . . . . 92
Exa 5.6 Order of Chebyshev Filter . . . . . . . . . . . . . . . . 92
Exa 5.7 Chebyshev Filter Design . . . . . . . . . . . . . . . . . 93
Exa 5.8 Order of Type 1 Low Pass Chebyshev Filter . . . . . . 93
Exa 5.9 Chebyshev Filter Design . . . . . . . . . . . . . . . . . 94
Exa 5.10 HPF Filter Design with given Specications . . . . . . 94
Exa 5.11 Impulse Invariant Method Filter Design . . . . . . . . 95
Exa 5.16 Bilinear Transformation Method Filter Design . . . . . 97
Exa 5.17 HPF Design using Bilinear Transform . . . . . . . . . 98
Exa 5.18 Bilinear Transformation Method Filter Design . . . . . 99
Exa 5.19 Single Pole LPF into BPF Conversion . . . . . . . . . 99
Exa 5.29 Pole Zero IIR Filter into Lattice Ladder Structure . . 100
Exa 6.1 Group Delay and Phase Delay . . . . . . . . . . . . . 102
Exa 6.5 LPF Magnitude Response . . . . . . . . . . . . . . . . 104
Exa 6.6 HPF Magnitude Response . . . . . . . . . . . . . . . . 104
Exa 6.7 BPF Magnitude Response . . . . . . . . . . . . . . . . 106
Exa 6.8 BRF Magnitude Response . . . . . . . . . . . . . . . . 108
Exa 6.9.a HPF Magnitude Response using Hanning Window . . 110
Exa 6.9.b HPF Magnitude Response using Hamming Window . . 112
Exa 6.10 Hanning Window Filter Design . . . . . . . . . . . . . 114
Exa 6.11 LPF Filter Design using Kaiser Window . . . . . . . . 116
Exa 6.12 BPF Filter Design using Kaiser Window . . . . . . . . 118
Exa 6.13.a Digital Dierentiator using Rectangular Window . . . 122
Exa 6.13.b Digital Dierentiator using Hamming Window . . . . 124
Exa 6.14.a Hilbert Transformer using Rectangular Window . . . . 124
Exa 6.14.b Hilbert Transformer using Blackman Window . . . . . 126
Exa 6.15 Filter Coecients obtained by Sampling . . . . . . . . 127
7
Exa 6.16 Coecients of Linear phase FIR Filter . . . . . . . . . 128
Exa 6.17 BPF Filter Design using Sampling Method . . . . . . 128
Exa 6.18.a Frequency Sampling Method FIR LPF Filter . . . . . 129
Exa 6.18.b Frequency Sampling Method FIR LPF Filter . . . . . 130
Exa 6.19 Filter Coecients Determination . . . . . . . . . . . . 131
Exa 6.20 Filter Coecients using Hamming Window . . . . . . 133
Exa 6.21 LPF Filter using Rectangular Window . . . . . . . . . 135
Exa 6.28 Filter Coecients for Direct Form Structure . . . . . . 137
Exa 6.29 Lattice Filter Coecients Determination . . . . . . . . 138
Exa 7.2 Subtraction Computation . . . . . . . . . . . . . . . . 139
Exa 7.14 Variance of Output due to AD Conversion Process . . 139
Exa 8.9 Two Component Decomposition . . . . . . . . . . . . 141
Exa 8.10 Two Band Polyphase Decomposition . . . . . . . . . . 142
Exa 9.7.a Frequency Resolution Determination . . . . . . . . . . 143
Exa 9.7.b Record Length Determination . . . . . . . . . . . . . . 144
Exa 9.8.a Smallest Record Length Computation . . . . . . . . . 144
Exa 9.8.b Quality Factor Computation . . . . . . . . . . . . . . 145
Exa 11.3 Program for Integer Multiplication . . . . . . . . . . . 146
Exa 11.5 Function Value Calculation . . . . . . . . . . . . . . . 146
8
List of Figures
1.1 Continuous Time Plot and Discrete Time Plot . . . . . . . . 12
1.2 Continuous Time Plot and Discrete Time Plot . . . . . . . . 13
1.3 Determining Periodicity of Signal . . . . . . . . . . . . . . . 17
1.6 Plot Magnitude and Phase Response . . . . . . . . . . . . . 23
1.7 Sketch Magnitude and Phase Response . . . . . . . . . . . . 25
1.8 Plot Magnitude and Phase Response . . . . . . . . . . . . . 26
1.9 Filter to Eliminate High Frequency Component . . . . . . . 28
1.10 Frequency Response of LTI System . . . . . . . . . . . . . . 30
1.11 Frequency Response of LTI System . . . . . . . . . . . . . . 32
2.1 Pole Zero Plot of the Dierence Equation . . . . . . . . . . . 42
2.2 Frequency Response of the System . . . . . . . . . . . . . . 43
2.3 Impulse Response of the System . . . . . . . . . . . . . . . . 46
2.4 Pole Zero Plot of the System . . . . . . . . . . . . . . . . . . 48
2.5 Unit Sample Response of the System . . . . . . . . . . . . . 49
2.6 Determine Output Response . . . . . . . . . . . . . . . . . . 51
2.7 Pole Zero Pattern of the System . . . . . . . . . . . . . . . . 54
2.8 Impulse Response of the System . . . . . . . . . . . . . . . . 57
3.1 DFT of the Sequence . . . . . . . . . . . . . . . . . . . . . . 60
3.2 Plot the Sequence . . . . . . . . . . . . . . . . . . . . . . . . 63
3.3 Sketch Sequence . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.1 LPF Magnitude Response . . . . . . . . . . . . . . . . . . . 103
6.2 HPF Magnitude Response . . . . . . . . . . . . . . . . . . . 105
6.3 BPF Magnitude Response . . . . . . . . . . . . . . . . . . . 107
6.4 BRF Magnitude Response . . . . . . . . . . . . . . . . . . . 109
9
6.5 HPF Magnitude Response using Hanning Window . . . . . . 111
6.6 HPF Magnitude Response using Hamming Window . . . . . 113
6.7 Hanning Window Filter Design . . . . . . . . . . . . . . . . 115
6.8 LPF Filter Design using Kaiser Window . . . . . . . . . . . 117
6.9 BPF Filter Design using Kaiser Window . . . . . . . . . . . 119
6.10 Digital Dierentiator using Rectangular Window . . . . . . . 121
6.11 Digital Dierentiator using Hamming Window . . . . . . . . 123
6.12 Hilbert Transformer using Rectangular Window . . . . . . . 125
6.13 Hilbert Transformer using Blackman Window . . . . . . . . 126
6.14 Frequency Sampling Method FIR LPF Filter . . . . . . . . . 129
6.15 Frequency Sampling Method FIR LPF Filter . . . . . . . . . 131
6.16 Filter Coecients Determination . . . . . . . . . . . . . . . 132
6.17 Filter Coecients using Hamming Window . . . . . . . . . . 134
6.18 LPF Filter using Rectangular Window . . . . . . . . . . . . 136
10
Chapter 1
DISCRETE TIME SIGNALS
AND LINEAR SYSTEMS
Scilab code Exa 1.1 Continuous Time Plot and Discrete Time Plot
1 // Example 1 . 1
2 // Sket ch t he c ont i nuous ti me s i g n a l x ( t ) =2exp(2 t )
and a l s o i t s d i s c r e t e ti me e q ui v a l e nt s i g n a l wi th
a s ampl i ng pe r i o d T = 0 . 2 s e c
3 clear all;
4 clc ;
5 close ;
6 t=0:0.01:2;
7 x1=2*exp(-2*t);
8 subplot (1,2,1);
9 plot(t,x1);
10 xlabel( t );
11 ylabel( x ( t ) );
12 title( CONTINUOUS TIME PLOT );
13 n=0:0.2:2;
14 x2=2*exp(-2*n);
15 subplot (1,2,2);
11
Figure 1.1: Continuous Time Plot and Discrete Time Plot
12
Figure 1.2: Continuous Time Plot and Discrete Time Plot
16 plot2d3(n,x2);
17 xlabel( n );
18 ylabel( x ( n) );
19 title( DISCRETE TIME PLOT );
Scilab code Exa 1.2 Continuous Time Plot and Discrete Time Plot
1 // Example 1 . 2
2 // Sket ch t he c ont i nuous ti me s i g n a l x=s i n ( 7 t )+s i n
( 10 t ) and a l s o i t s d i s c r e t e ti me e q ui v a l e nt
s i g n a l wi th a s ampl i ng pe r i o d T = 0 . 2 s e c
13
3 clear all;
4 clc ;
5 close ;
6 t=0:0.01:2;
7 x1=sin(7*t)+sin (10*t);
8 subplot (1,2,1);
9 plot(t,x1);
10 xlabel( t );
11 ylabel( x ( t ) );
13 n=0:0.2:2;
14 x2=sin(7*n)+sin (10*n);
15 subplot (1,2,2);
16 plot2d3(n,x2);
17 xlabel( n );
Scilab code Exa 1.3.a Evaluate the Summations
1 // Example 1 . 3 ( a )
2 //MAXIMA SCILAB TOOLBOX REQUIRED FOR THIS PROGRAM
3 // Ca l c ul a t e Fol l owi ng Summations
4 clear all;
5 clc ;
6 close ;
7 syms n;
8 X= symsum (sin (2*n),n ,2, 2);
9 // Di s pl ay t he r e s u l t i n command window
10 disp (X,The Val ue o f summati on comes out t o be : );
Scilab code Exa 1.3.b Evaluate the Summations
14
1 // Example 1 . 3 ( b)
3 // Ca l c ul a t e Fol l owi ng Summations
4 clear all;
5 clc ;
6 close ;
7 syms n;
8 X= symsum (%e^(2*n),n ,0, 0);
10 disp (X,The Val ue o f summati on comes out t o be : );
Scilab code Exa 1.4.a Check for Energy or Power Signals
1 // Example 1 . 4 ( a )
3 // Fi nd Energy and Power o f Gi ven S i g n a l s
4 clear all;
5 clc ;
6 close ;
7 syms n N;
8 x=(1/3)^n;
9 E= symsum (x^2,n ,0, %inf);
11 disp (E, Energy : );
12 p=(1/(2*N+1))*symsum (x^2,n ,0, N);
13 P=limit(p,N,%inf);
14 disp (P,Power : );
15 //The Energy i s Fi n i t e and Power i s 0 . The r e f or e t he
gi ve n s i g n a l i s an Energy Si g na l
Scilab code Exa 1.4.d Check for Energy or Power Signals
1 // Example 1 . 4 ( d)
15
3 // Fi nd Energy and Power o f Gi ven S i g n a l s
4 clear all;
5 clc ;
6 close ;
7 syms n N;
8 x=%e^(2*n);
9 E= symsum (x^2,n ,0, %inf);
11 disp (E, Energy : );
12 p=(1/(2*N+1))*symsum (x^2,n ,0, N);
13 P=limit(p,N,%inf);
14 disp (P,Power : );
15 //The Energy andPower i s i n f i n i t e . The r e f or e t he
gi ve n s i g n a l i s an n e i t h e r Energy Si g na l nor
Power Si g na l
Scilab code Exa 1.5.a Determining Periodicity of Signal
1 // Example 1 . 5 ( a )
2 //To Determi ne Whether Gi ven Si g na l i s Pe r i o di c or
not
3 clear all;
4 clc ;
5 close ;
6 t=0:0.01:2;
7 x1=exp(%i*6*%pi*t);
8 subplot (1,2,1);
9 plot(t,x1);
10 xlabel( t );
11 ylabel( x ( t ) );
13 n=0:0.2:2;
16
Figure 1.3: Determining Periodicity of Signal
17
14 x2=exp(%i*6*%pi*n);
15 subplot (1,2,2);
16 plot2d3(n,x2);
17 xlabel( n );
20 // Hence Gi ven Si g na l i s Pe r i o di c wi th N=1
Scilab code Exa 1.5.c Determining Periodicity of Signal
1 // Example 1 . 5 ( c )
18
not
3 clear all;
4 clc ;
5 close ;
6 t=0:0.01:10;
7 x1=cos(2* %pi*t/3);
8 subplot (1,2,1);
9 plot(t,x1);
10 xlabel( t );
11 ylabel( x ( t ) );
13 n=0:0.2:10;
14 x2=cos(2* %pi*n/3);
15 subplot (1,2,2);
16 plot2d3(n,x2);
17 xlabel( n );
Scilab code Exa 1.5.d Determining Periodicity of Signal
1 // Example 1 . 5 ( d)
not
3 clear all;
4 clc ;
5 close ;
6 t=0:0.01:50;
7 x1=cos(%pi*t/3)+cos (3*%pi*t/4);
8 subplot (1,2,1);
9 plot(t,x1);
19
20
10 xlabel( t );
11 ylabel( x ( t ) );
13 n=0:0.2:50;
14 x2=cos(%pi*n/3)+cos (3*%pi*n/4);
15 subplot (1,2,2);
16 plot2d3(n,x2);
17 xlabel( n );
Scilab code Exa 1.11 Stability of the System
1 // Example 1. 11
3 // Te s t i ng S t a b i l i t y o f Gi ven System
4 clear all;
5 clc ;
6 close ;
7 syms n;
8 x =(1/2)^n
9 X= symsum (x,n ,0, %inf );
11 disp (X,Summation i s : );
12 disp( Hence Summation < i n f i n i t y . Gi ven System i s
St a bl e );
Scilab code Exa 1.12 Convolution Sum of Two Sequences
1 // Example 1. 12
2 // Program t o Compute c o nvo l ut i o n o f gi ve n s e que nc e s
3 //x ( n) =[ 3 2 1 2 ] , h( n) =[ 1 2 1 2 ] ;
21
4 clear all;
5 clc ;
6 close ;
7 x=[3 2 1 2];
8 h=[1 2 1 2];
9 y=convol(x,h);
10 disp(y);
Scilab code Exa 1.13 Convolution of Two Signals
1 // Example 1. 13
2 // Program t o Compute c o nvo l ut i o n o f gi ve n s e que nc e s
3 //x ( n) =[ 1 2 1 1 ] , h( n) =[ 1 1 1 1] ;
4 clear all;
5 clc ;
6 close ;
7 x=[1 2 1 1];
8 h=[1 -1 1 -1];
9 y=convol(x,h);
10 disp(round(y));
Scilab code Exa 1.18 Cross Correlation of Two Sequences
1 // Example 1. 18
2 // Program t o Compute Cross c o r r e l a t i o n o f gi ve n
s e que nc e s
3 //x ( n) =[ 1 2 1 1 ] , h( n) =[ 1 1 2 1 ] ;
4 clear all;
5 clc ;
6 close ;
7 x=[1 2 1 1];
8 h=[1 1 2 1];
9 h1=[1 2 1 1];
22
Figure 1.6: Plot Magnitude and Phase Response
10 y=convol(x,h1);
11 disp(round(y));
Scilab code Exa 1.19 Determination of Input Sequence
1 // Example 1. 19
2 //To f i nd i nput x ( n)
3 //h( n) =[ 1 2 1 ] , y ( n) =[ 1 5 10 11 8 4 1 ]
4 clear all;
5 clc ;
6 close ;
7 z=%z;
8 a=z^6+5*(z^(5))+10*(z^(4))+11*(z^(3))+8*(z^(2))+4*(z
^(1))+1;
9 b=z^6+2*z^(5) +1*z^(4);
10 x =ldiv(a,b,5);
11 disp (x,x ( n)=);
23
Scilab code Exa 1.32.a Plot Magnitude and Phase Response
1 // Example 1. 32
2 // Program t o Pl ot Magni tude and Phase Responce
3 clear all;
4 clc ;
5 close ;
6 w=-%pi :0.01: %pi;
7 H=1+2* cos(w)+2*cos(2*w);
8 // c a l u c u l a t i o n o f Phase and Magni tude o f H
9 [phase_H ,m]= phasemag(H);
10 Hm=abs(H);
11 a=gca();
12 subplot (2,1,1);
13 a.y_location= o r i g i n ;
14 plot2d(w/%pi ,Hm);
15 xlabel( Frequency i n Radi ans )
16 ylabel( abs (Hm) );
17 title( MAGNITUDE RESPONSE );
18 subplot (2,1,2);
19 a=gca();
20 a.x_location= o r i g i n ;
22 plot2d(w/(2* %pi),phase_H);
23 xlabel( Frequency i n Radi ans );
24 ylabel( <(H) );
25 title( PHASE RESPONSE );
Scilab code Exa 1.37 Sketch Magnitude and Phase Response
24
Figure 1.7: Sketch Magnitude and Phase Response
1 // Example 1. 37
3 //y ( n) =1/2[ x ( n)+x ( n2) ]
4 clear all;
5 clc ;
6 close ;
7 w=0:0.01: %pi;
8 H=(1+ cos(2*w)-%i*sin(2*w))/2;
11 Hm=abs(H);
12 a=gca();
13 subplot (2,1,1);
19 subplot (2,1,2);
20 a=gca();
25
Figure 1.8: Plot Magnitude and Phase Response
25 ylabel( <(H) );
Scilab code Exa 1.38 Plot Magnitude and Phase Response
1 // Example 1. 38
3 // 0 . 5 de l t a ( n)+de l t a ( n1) +0. 5 de l t a ( n2)
4 clear all;
5 clc ;
6 close ;
7 w=-%pi :0.01: %pi;
8 H=0.5+ exp(-%i*w)+0.5* exp(-%i*w);
11 Hm=abs(H);
12 a=gca();
13 subplot (2,1,1);
26
19 subplot (2,1,2);
20 a=gca();
25 ylabel( <(H) );
Scilab code Exa 1.45 Filter to Eliminate High Frequency Component
1 // Example 1. 45
2 //
3 clear all;
4 clc ;
5 close ;
6 t=0:0.01:10;
7 x=2*cos(5*t)+cos (300*t);
8 x1=2*cos (5*t);
9 b=[0.05 0.05];
10 a=[1 -0.9];
11 y=filter(b,a,x);
12 subplot (2,1,1);
13 plot(t,x);
14 xlabel( Time i n Sec );
15 ylabel( Ampl i tude );
16 subplot (2,1,2);
17 plot(t,y);
27
Figure 1.9: Filter to Eliminate High Frequency Component
28
18 subplot (2,1,2);
19 plot(t,x1, : );
20 title( x : SIGNAL WITHOUT NOISE y : SIGNAL WITH NOISE )
;
21 xlabel( Time i n Sec );
22 ylabel( Ampl i tude );
Scilab code Exa 1.57.a Discrete Convolution of Sequences
1 // Example 1. 57 ( a )
2 // Program t o Compute d i s c r e t e c o nvo l ut i o n o f gi ve n
s e que nc e s
3 //x ( n) =[ 1 2 1 1 ] , h( n) =[ 1 0 1 1 ] ;
4 clear all;
5 clc ;
6 close ;
7 x=[1 2 -1 1];
8 h=[1 0 1 1];
9 y=convol(x,h);
10 disp(round(y));
Scilab code Exa 1.61 Fourier Transform
1 // Example 1. 61
3 // Four i e r t r ans f or m o f ( 3 ) n u( n)
4 clear all;
5 clc ;
6 close ;
7 syms n;
8 x =(3) ^n;
9 X= symsum (x,n ,0, %inf )
29
Figure 1.10: Frequency Response of LTI System
11 disp (X,The Four i e r Transf orm does not e x i t as x ( n)
i s not a b s o l u t e l y summable and appr oac he s
i n f i n i t y i . e . );
Scilab code Exa 1.62 Fourier Transform
1 // Example 1. 62
3 // Four i e r t r ans f or m o f ( 0 . 8 ) | n | u( n)
4 clear all;
5 clc ;
6 close ;
7 syms w n;
8 X= symsum ((0.8)^n*%e^(%i*w*n),n ,1, %inf )+symsum
((0.8)^n*%e^(-%i*w*n),n ,0, %inf )
10 disp (X,The Four i e r Transf orm comes out t o be : );
30
Scilab code Exa 1.64.a Frequency Response of LTI System
1 // Example 1. 64 ( a )
2 // Program t o Ca l c ul a t e Pl ot Magni tude and Phase
Responce
3 clear all;
4 clc ;
5 close ;
6 w=0:0.01: %pi;
7 H=1/(1 -0.5*%e^(-%i*w));
10 Hm=abs(H);
11 a=gca();
12 subplot (2,1,1);
18 subplot (2,1,2);
19 a=gca();
24 ylabel( <(H) );
Scilab code Exa 1.64.c Frequency Response of LTI System
31
Figure 1.11: Frequency Response of LTI System
32
1 // Example 1. 64 ( c )
2 // Program t o Ca l c ul a t e Pl ot Magni tude and Phase
Responce
3 clear all;
4 clc ;
5 close ;
6 w=0:0.01: %pi;
7 H=1/(1 -0.9*%i*%e^(-%i*w));
10 Hm=abs(H);
11 a=gca();
12 subplot (2,1,1);
18 subplot (2,1,2);
19 a=gca();
24 ylabel( <(H) );
33
Chapter 2
THE Z TRANSFORM
Scilab code Exa 2.1 z Transform and ROC of Causal Sequence
1 // Example 2 . 1
2 //Z t r ans f or m o f [ 1 0 3 1 2 ]
3 clear all;
4 clc ;
5 close ;
6 function[za]= ztransfer(sequence ,n)
7 z=poly(0, z , r )
8 za=sequence *(1/z)^n
9 endfunction
10 x1=[1 0 3 -1 2];
11 n=0: length(x1) -1;
12 zz=ztransfer(x1,n);
14 disp (zz,Zt r ans f or m o f s e que nc e i s : );
15 disp( ROC i s t he e n t i r e pl ane e xc e pt z = 0 );
Scilab code Exa 2.2 z Transform and ROC of Anticausal Sequence
34
1 // Example 2 . 2
2 //Z t r ans f or m o f [ 3 2 1 0 1 ]
3 clear all;
4 clc ;
5 close ;
7 z=poly(0, z , r )
9 endfunction
10 x1=[-3 -2 -1 0 1];
11 n=-(length(x1) -1):0;
15 disp( ROC i s t he e n t i r e pl ane e xc e pt z = %i nf );
Scilab code Exa 2.3 z Transform of the Sequence
1 // Example 2 . 3
2 //Z t r ans f or m o f [ 2 1 3 2 1 0 2 3 1]
3 clear all;
4 clc ;
5 close ;
7 z=poly(0, z , r )
9 endfunction
10 x1=[2 -1 3 2 1 0 2 3 -1];
11 n=-4:4;
15 disp( ROC i s t he e n t i r e pl ane e xc e pt z = 0 and z =
%i nf );
35
Scilab code Exa 2.4 z Transform and ROC of the Signal
1 // Example 2 . 4
3 //Z t r ans f or m o f an u( n)
4 clear all;
5 clc ;
6 close ;
7 syms a n z;
8 x =a^n
9 X= symsum (x*(z^(-n)),n ,0, %inf );
11 disp (X,Zt r ans f or m o f an u( n) wi th i s : );
12 disp( ROC i s t he Regi on mod( z ) > a )
Scilab code Exa 2.5 z Transform and ROC of the Signal
1 // Example 2 . 5
3 //Z t r ans f or m o f bn u(n1)
4 clear all;
5 clc ;
6 close ;
7 syms b n z;
8 x =b^n
9 X= symsum (x*(z^(-n)),n ,0, %inf );
11 disp (X,Zt r ans f or m o f bn u( n) wi th i s : );
12 disp( ROC i s t he Regi on mod( z ) < b )
36
Scilab code Exa 2.6 Stability of the System
1 // Example 2 . 6
3 //Z t r ans f or m o f 2n u( n)
4 clear all;
5 clc ;
6 close ;
7 syms n z;
8 x =(2) ^n
9 X= symsum (x*(z^(-n)),n ,0, %inf );
11 disp (X,Zt r ans f or m o f 2n u( n) i s : );
12 disp( ROC i s t he Regi on mod( z ) > 2 );
Scilab code Exa 2.7 z Transform of the Signal
1 // Example 2 . 7
3 //Z t r ans f or m o f [ 3 ( 3 n) 4(2) n ] u( n)
4 clear all;
5 clc ;
6 close ;
7 syms n z;
8 x1 =(3) ^(n);
9 X1= symsum (3* x1 *(z^(-n)),n ,0, %inf );
10 x2 =(4) ^(n);
11 X2= symsum (4* x2 *(z^(-n)),n ,0, %inf );
12 X = (X1 -X2);
14 disp (X,Zt r ans f or m o f [ 3 ( 3 n) 4(2) n ] u( n) i s : );
Scilab code Exa 2.8.a z Transform of the Signal
37
1 // Example 2 . 8 ( a )
3 //Z t r ans f or m o f c os (Won)
4 clc;
5 syms Wo n z;
6 x1=exp(sqrt(-1)*Wo*n);
7 X1=symsum(x1*(z^-n),n,0,%inf);
8 x2=exp(-sqrt(-1)*Wo*n);
10 X=(X1+X2)/2;
11 disp(X, X( z )= );
1 // Example 2 . 9
3 //Z t r ans f or m o f ( 1/3) n u( n1)
4 clear all;
5 clc ;
6 close ;
7 syms n z;
8 x =(1/3)^n;
9 X= (1/z)*symsum (x*(z^(-n)),n ,0, %inf );
11 disp (X,Zt r ans f or m o f ( 1/3) n u( n1) i s : );
Scilab code Exa 2.10 z Transform Computation
1 // Example 2. 10
3 //Z t r ans f or m o f r n . c os (Won)
4 clc;
5 syms r Wo n z;
38
6 x1=(r^n)*exp(sqrt(-1)*Wo*n);
8 x2=(r^n)*exp(-sqrt(-1)*Wo*n);
10 X=(X1+X2)/2;
11 disp(X, X( z )= );
1 // Example 2. 11
3 //Z t r ans f or m o f n . an u( n)
4 clear all;
5 clc ;
6 close ;
7 syms a n z;
8 x =(a) ^n;
9 X= symsum (x*(z^(-n)),n ,0, %inf )
10 Y = diff (X,z);
12 disp (Y,Zt r ans f or m o f n . an u( n) i s : );
Scilab code Exa 2.13.a z Transform of Discrete Time Signals
1 // Example 2. 13 ( a )
3 //Z t r ans f or m o f ( 1/5) n u( n) +5(1/2) (n) u(n1)
4 clear all;
5 clc ;
6 close ;
7 syms n z;
8 x1 =( -1/5)^n ;
9 X1= symsum (x1 *(z^(-n)),n ,0, %inf );
39
10 x2 =(1/2)^(-n);
11 X2= symsum (5* x2 *(z^(-n)),n ,0, %inf );
12 X = (X1 -X2);
14 disp (X,Zt r ans f or m o f [ 3 ( 3 n) 4(2) n ] u( n) i s : );
15 disp( ROC i s t he Regi on 1/5 < mod( z ) < 2 );
Scilab code Exa 2.13.b z Transform of Discrete Time Signals
1 // Example 2. 13 ( b)
3 //Z t r ans f or m
4 clc;
5 syms n z k;
6 x1=1;
7 X1=symsum(x1*z^(-n),n,0,0);
8 x2=1;
9 X2=symsum(x2*z^(-n),n,1,1);
10 x3=1;
11 X3=symsum(x3*z^(-n),n,2,2);
12 X=0.5* X1+X2 -1/3* X3;
13 disp(X, X( z )= );
Scilab code Exa 2.13.c z Transform of Discrete Time Signals
1 // Example 2. 13 ( c )
3 //Z t r ans f or m o f u( n2)
4 clear all;
5 clc ;
6 close ;
7 syms n z;
8 x =1;
40
9 X= (1/(z^2))*symsum (x*(z^(-n)),n ,0, %inf );
11 disp (X,Zt r ans f or m o f u( n2) i s : );
Scilab code Exa 2.13.d z Transform of Discrete Time Signals
1 // Example 2. 13 ( d)
3 //Z t r ans f or m o f ( n+0. 5) ( ( 1 / 3 ) n) u( n)
4 clear all;
5 clc ;
6 close ;
7 syms n z;
8 x1 =(1/3)^n;
9 X11= symsum (x1*(z^(-n)),n ,0, %inf )
10 X1 = diff (X11 ,z);
11 x2 =(1/3) ^(n);
12 X2= symsum (0.5* x2 *(z^(-n)),n ,0, %inf );
13 X = (X1+X2);
15 disp (X,Zt r ans f or m o f ( n+0. 5) ( ( 1 / 3 ) n) u( n) i s : );
Scilab code Exa 2.16 Impulse Response of the System
1 // Example 2. 16
2 //To f i nd i nput h( n)
3 // a=[ 1 2 4 1 ] , b =[ 1]
4 clear all;
5 clc ;
6 close ;
7 z=%z;
8 a=z^3+2*(z^(2)) -4*(z)+1;
9 b=z^3;
41
Figure 2.1: Pole Zero Plot of the Dierence Equation
10 h =ldiv(a,b,4);
11 disp (h,h( n)=);
Scilab code Exa 2.17 Pole Zero Plot of the Dierence Equation
1 // Example 2. 17
2 //To draw t he pol e z e r o pl o t
3 clear all;
4 clc ;
5 close ;
6 z=%z
7 H1Z =((z)*(z-1))/((z -0.25) *(z-0.5));
8 xset( window ,1);
9 plzr(H1Z);
42
Figure 2.2: Frequency Response of the System
43
Scilab code Exa 2.19 Frequency Response of the System
1 // Example 2. 19
3 clear all;
4 clc ;
5 close ;
6 w=-%pi :0.01: %pi;
7 H=1/(1 -0.5*( cos(w)-%i*sin(w)));
10 Hm=abs(H);
11 a=gca();
12 subplot (2,1,1);
18 subplot (2,1,2);
19 a=gca();
24 ylabel( <(H) );
Scilab code Exa 2.20.a Inverse z Transform Computation
1 // Example 2. 10 ( a )
2 //To f i nd i nput h( n)
3 //X( z ) =( z +0. 2) / ( ( z +0. 5) ( z 1) ;
4 clear all;
5 clc ;
44
6 close ;
7 z=%z;
8 a=(z+0.5)*(z-1);
9 b=z+0.2;
10 h =ldiv(b,a,4);
11 disp (h,h( n)=);
Scilab code Exa 2.22 Inverse z Transform Computation
1 // Example 2. 22
3 //X( z ) =1/(2 z ( 2)+2z ( 1) +1) ;
4 clear all;
5 clc ;
6 close ;
7 z=%z;
8 a=(2+2*z+z^2);
9 b=z^2;
10 h =ldiv(b,a,6);
11 disp (h, Fi r s t s i x va l ue s o f h( n)=);
Scilab code Exa 2.23 Causal Sequence Determination
1 // Example 2. 23
3 //X( z ) =1/(12z ( 1) ) (1z ( 1) ) 2 ;
4 clear all;
5 clc ;
6 close ;
7 z=%z;
8 a=(z-2)*(z-1) ^2;
9 b=z^3;
10 h =ldiv(b,a,6);
45
Figure 2.3: Impulse Response of the System
11 disp (h, Fi r s t s i x va l ue s o f h( n)=);
1 // Example 2. 34
2 //To pl o t t he i mpul s e r e s ponc e o f t he system
a n a l y i c a l l y and us i ng s c i l a b
3 clear all;
4 clc ;
5 close ;
6 n=0:1:50;
46
7 x=[1, zeros (1,50)];
8 b=[1 2];
9 a=[1 -3 -4];
10 yanaly =6/5*4.^n-1/5*( -1).^n;// Ana l y t i c a l So l ut i o n
11 ymat=filter(b,a,x);
12 subplot (3,1,1);
13 plot2d3(n,x);
14 xlabel( n );
16 title( INPUT SEQUENCE (IMPULSE FUNCTION) );
17 subplot (3,1,2);
18 plot2d3(n,yanaly);
19 xlabel( n );
20 ylabel( y ( n) );
21 title( OUTPUT SEQUENCE yanal y );
22 subplot (3,1,3);
23 plot2d3(n,ymat);
24 xlabel( n );
26 title( OUTPUT SEQUENCE ymat );
27 //As t he An a l t i c a l Pl ot matches t he S c i l a b Pl ot
hence i t i s t he Responce o f t he system
Scilab code Exa 2.35.a Pole Zero Plot of the System
1 // Example 2. 35 ( a )
3 clear all;
4 clc ;
5 close ;
6 z=%z
7 H1Z=(z)/(z^2-z-1);
8 xset( window ,1);
47
Figure 2.4: Pole Zero Plot of the System
48
Figure 2.5: Unit Sample Response of the System
9 plzr(H1Z);
Scilab code Exa 2.35.b Unit Sample Response of the System
1 // Example 2. 35 ( b)
2 //To pl o t t he r e s ponc e o f t he system a n a l y i c a l l y and
us i ng s c i l a b
3 clear all;
4 clc ;
5 close ;
6 n=0:1:20;
49
7 x=ones(1,length(n));
8 b=[0 1];
9 a=[1 -1 -1];
10 yanaly =0.447*(1.618) .^n -0.447*( -0.618) .^n;//
Ana l y t i c a l So l ut i o n
11 [ymat ,zf]= filter(b, a, x);
12 subplot (3,1,1);
13 plot2d3(n,x);
14 xlabel( n );
16 title( INPUT SEQUENCE (STEP FUNCTION) );
17 subplot (3,1,2);
19 xlabel( n );
22 subplot (3,1,3);
23 plot2d3(n,ymat ,zf);
24 xlabel( n );
Scilab code Exa 2.38 Determine Output Response
1 // Example 2. 38
us i ng s c i l a b
3 clear all;
4 clc ;
5 close ;
6 n=0:1:20;
50
Figure 2.6: Determine Output Response
51
7 x=n;
8 b=[0 1 1];
9 a=[1 -0.7 0.12];
10 yanaly =38.89*(0.4) .^n -26.53*(0.3) .^n -12.36+4.76*n;//
Ana l y t i c a l So l ut i o n
12 subplot (3,1,1);
13 plot2d3(n,x);
14 xlabel( n );
16 title( INPUT SEQUENCE (RAMP FUNCTION) );
17 subplot (3,1,2);
19 xlabel( n );
22 subplot (3,1,3);
23 plot2d3(n,ymat);
24 xlabel( n );
Scilab code Exa 2.40 Input Sequence Computation
1 // Example 2. 40
3 //h( n)=1 2 3 2 , y ( n) =[ 1 3 7 10 10 7 2 ]
4 clear all;
5 clc ;
6 close ;
7 z=%z;
8 a=z^6+3*(z^(5))+7*(z^(4))+10*(z^(3))+10*(z^(2))+7*(z
^(1))+2;
52
9 b=z^6+2*z^(5) +3*z^(4) +2*z^(3);
10 x =ldiv(a,b,4);
11 disp (x,x ( n)=);
Scilab code Exa 2.41.a z Transform of the Signal
1 // Example 2. 41 ( a )
3 //Z t r ans f or m o f n. ( 1) n u( n)
4 clear all;
5 clc ;
6 close ;
7 syms a n z;
8 x =(-1) ^n;
9 X= symsum (x*(z^(-n)),n ,0, %inf )
10 Y = diff (X,z);
12 disp (Y,Zt r ans f or m o f n. ( 1) n u( n) i s : );
Scilab code Exa 2.41.b z Transform of the Signal
1 // Example 2. 41 ( b)
3 //Z t r ans f or m o f n2 u( n)
4 clear all;
5 clc ;
6 close ;
7 syms n z;
8 x =1;
9 X= symsum (x*(z^(-n)),n ,0, %inf )
10 Y = diff(diff (X,z),z);
12 disp (Y,Zt r ans f or m o f n2 u( n) i s : );
53
Figure 2.7: Pole Zero Pattern of the System
Scilab code Exa 2.41.c z Transform of the Signal
1 // Example 2. 41 ( c )
3 //Z t r ans f or m o f ( 1) n . c os ( %pi /3n)
4 clc;
5 syms n z;
6 Wo=%pi/3;
8 X1=(-1)^n*symsum(x1*(z^-n),n,0,%inf);
10 X2=(-1)^n*symsum(x2*(z^-n),n,0,%inf);
11 X=(X1+X2)/2;
12 disp(X, X( z )= );
54
Scilab code Exa 2.45 Pole Zero Pattern of the System
1 // Example 2. 45
3 clear all;
4 clc ;
5 close ;
6 z=%z
7 H1Z =((z)*(z+1))/(z^2-z+0.5);
8 xset( window ,1);
9 plzr(H1Z);
Scilab code Exa 2.53.a z Transform of the Sequence
1 // Example 2. 53 ( a )
2 //Z t r ans f or m o f [ 3 1 2 5 7 0 1 ]
3 clear all;
4 clc ;
5 close ;
7 z=poly(0, z , r )
9 endfunction
10 x1=[3 1 2 5 7 0 1];
11 n=-3:3;
Scilab code Exa 2.53.b z Transform of the Signal
1 // Example 2. 53 ( b)
55
3 //Z t r ans f or m o f de l t a ( n)
4 clc;
5 syms n z;
6 x=1;
7 X=symsum(x*z^(-n),n,0,0);
8 disp(X, X( z )= );
Scilab code Exa 2.53.c z Transform of the Signal
1 // Example 2. 53 ( c )
4 clc;
5 syms n z k;
6 x=1;
7 X=symsum(x*z^(-n),n,k,k);
8 disp(X, X( z )= );
Scilab code Exa 2.53.d z Transform of the Signal
1 // Example 2. 53 ( d)
4 clc;
5 syms n z kc;
6 x=1;
7 X=symsum(x*z^(-n),n,-k,-k);
8 disp(X, X( z )= );
Scilab code Exa 2.54 z Transform of Cosine Signal
56
Figure 2.8: Impulse Response of the System
1 // Example 2. 54
3 //Z t r ans f or m o f c os (Won)
4 clc;
5 syms Wo n z;
10 X=(X1+X2)/2;
11 disp(X, X( z )= );
57
1 // Example 2. 58
us i ng s c i l a b
3 clear all;
4 clc ;
5 close ;
6 n=0:1:20;
7 x=[1 zeros (1,20)];
8 b=[1 -0.5];
9 a=[1 -1 3/16];
10 yanaly =0.5*(0.75) .^n+0.5*(0.25) .^n;// Ana l y t i c a l
So l ut i o n
12 subplot (3,1,1);
13 plot2d3(n,x);
14 xlabel( n );
16 title( INPUT SEQUENCE (IMPULSE FUNCTION) );
17 subplot (3,1,2);
19 xlabel( n );
22 subplot (3,1,3);
23 plot2d3(n,ymat);
24 xlabel( n );
58
Chapter 3
THE DISCRETE FOURIER
TRANSFORM
Scilab code Exa 3.1 DFT and IDFT
1 // Example 3 . 1
2 // Program t o Compute t he DFT o f a Sequence x [ n
] =[ 1 , 1 , 0 , 0 ]
3 // and IDFT o f a Sequence Y[ k ] =[ 1 , 0 , 1 , 0 ]
4 clear all;
5 clc ;
6 close ;
7 x = [1,1,0,0];
8 //DFT Computati on
9 X = fft (x , -1);
10 Y = [1,0,1,0];
11 //IDFT Computati on
12 y = fft (Y , 1);
13 // Di s pl ay s e que nc e X[ k ] and y [ n ] i n command window
14 disp(X,X[ k]=);
15 disp(y,y [ n]=);
59
Figure 3.1: DFT of the Sequence
Scilab code Exa 3.2 DFT of the Sequence
1 // Example 3 . 2
2 // Program t o Compute t he DFT o f a Sequence x [ n] =1 ,
0<=n<=2; and 0 o t he r wi s e
3 // f o r N=4 and N=8. Pl ot Magni tude and phase p l o t s o f
each .
4 clear all;
5 clc ;
6 close ;
7 //N=4
8 x1 = [1,1,1,0];
10 X1 = fft (x1 , -1);
11 //N=8
12 x2 = [1,1,1,0,0,0,0,0];
14 X2 = fft (x2 , -1);
15 // Di s pl ay s e que nc e X1 [ k ] and X2 [ k ] i n command window
16 disp(X1,X1 [ k]=);
60
17 disp(X2,X2 [ k]=);
18 // Pl o t s f o r N=4
19 n1 =0:1:3;
20 subplot (2,2,1);
21 a = gca ();
22 a.y_location = o r i g i n ;
23 a.x_location = o r i g i n ;
24 plot2d3(n1,abs(X1) ,2);
25 poly1=a.children (1).children (1);
26 poly1.thickness =2;
27 xtitle( N=4 , k , | X1( k ) | );
28 subplot (2,2,2);
29 a = gca ();
32 plot2d3(n1,atan(imag(X1),real(X1)) ,5);
35 xtitle( N=4 , k , <X1( k ) );
36 // Pl o t s f o r N=8
37 n2 =0:1:7;
38 subplot (2,2,3);
39 a = gca ();
42 plot2d3(n2,abs(X2) ,2);
45 xtitle( N=8 , k , | X2( k ) | );
46 subplot (2,2,4);
47 a = gca ();
50 plot2d3(n2,atan(imag(X2),real(X2)) ,5);
53 xtitle( N=8 , k , <X2( k ) );
61
Scilab code Exa 3.3 8 Point DFT
1 // Example 3 . 3
2 // Program t o Compute t he 8poi nt DFT o f t he Sequence
x [ n ] =[ 1 , 1 , 1 , 1 , 1 , 1 , 0 , 0 ]
3 clear all;
4 clc ;
5 close ;
6 x = [1,1,1,1,1,1,0,0];
8 X = fft (x , -1);
9 // Di s pl ay s e que nc e X[ k ] i n command window
10 disp(X,X[ k]=);
Scilab code Exa 3.4 IDFT of the given Sequence
1 // Example 3 . 4
2 // Program t o Compute t he IDFT o f t he Sequence X[ k
] =[ 5 , 0 , 1 j , 0 , 1 , 0 , 1+ j , 0 ]
3 clear all;
4 clc ;
5 close ;
6 j=sqrt(-1);
7 X = [5,0,1-j,0,1,0,1+j,0]
9 x = fft (X , 1);
10 // Di s pl ay s e que nc e s x [ n ] i n command window
11 disp(x,x [ n]=);
62
Figure 3.2: Plot the Sequence
Scilab code Exa 3.7 Plot the Sequence
1 // Example 3 . 7
2 // Program t o Compute c i r c u l a r c o nvo l ut i o n o f
f o l l o wi n g s e que nc e s
3 //x [ n ] =[ 1 , 2 , 2 , 1 , 0 ]
4 //Y[ k]=exp(j 4 pi k /5) . X[ k ]
5 clear all;
6 clc ;
7 close ;
8 x=[1,2,2,1,0];
9 X=fft(x,-1);
10 k=0:1:4;
11 j=sqrt(-1);
12 pi =22/7;
63
13 H=exp(-j*4*pi*k/5);
14 Y=H.*X;
16 y=fft(Y,1);
17 // Di s pl ay s e que nc e y [ n ] i n command window
18 disp(round(y),y [ n]=);
19 // Pl o t s
20 n=0:1:4;
21 a = gca ();
24 plot2d3(n,round(y) ,5);
27 xtitle( Pl ot o f s e que nc e y [ n ] , n , y [ n ] );
Scilab code Exa 3.9 Remaining Samples
1 // Example 3 . 9
2 // Program t o r e mai ni ng s ampl es o f t he s e que nc e
3 //X( 0 ) =12 ,X( 1 )=1+j 3 , X( 2 )=3+j 4 , X( 3 ) =1j 5 , X( 4 )=2+j 2 ,
X( 5 )=6+j 3 , X( 6 )=2j 3 , X( 7 ) =10
4 clear all;
5 clc ;
6 close ;
7 j=sqrt(-1);
8 z=1;
9 X(0+z)=12,X(1+z)=-1+j*3,X(2+z)=3+j*4,X(3+z)=1-j*5,X
(4+z)=-2+j*2,X(5+z)=6+j*3,X(6+z)=-2-j*3,X(7+z)
=10;
10 for a=9:1:14 do X(a)=conj(X(16-a)), end;
11 // Di s pl ay t he compl et e s e que nc e X[ k ] i n command
window
12 disp(X,X[ k]=);
64
Scilab code Exa 3.11 DFT Computation
1 // Example 3. 11
2 // Program t o Compute t he 8poi nt DFT o f t he
3 // x1 [ n ] =[ 1 , 0 , 0 , 0 , 0 , 1 , 1 , 1 ]
4 // x2 [ n ] =[ 0 , 0 , 1 , 1 , 1 , 1 , 0 , 0 ]
5 clear all;
6 clc ;
7 close ;
8 x1=[1,0,0,0,0,1,1,1];
9 x2=[0,0,1,1,1,1,0,0];
11 X1 = fft (x1 , -1);
12 X2 = fft (x2 , -1);
13 // Di s pl ay s e que nc e s X1 [ k ] and X2 [ k ] i n command
window
14 disp(X1,X1 [ k]=);
15 disp(X2,X2 [ k]=);
Scilab code Exa 3.13 Circular Convolution
1 // Example 3. 13
3 // x1 [ n] =[ 1 , 1 , 2 , 3 , 1]
4 // x2 [ n ] =[ 1 , 2 , 3 ]
5 clear all;
6 clc ;
7 close ;
8 x1=[1,-1,-2,3,-1];
9 x2=[1,2 ,3];
65
10 // Loop f o r z e r o paddi ng t he s ma l l e r s e que nc e out o f
t he two
11 n1=length(x1);
12 n2=length(x2);
13 n3=n2-n1;
14 if (n3 >=0) then
15 x1=[x1,zeros(1,n3)];
16 else
17 x2=[x2,zeros(1,-n3)];
18 end
20 X1=fft(x1 ,-1);
21 X2=fft(x2 ,-1);
22 Y=X1.*X2;
24 y=fft(Y,1);
26 disp(y,y [ n]=);
1 // Example 3. 14
3 // x1 [ n ] =[ 1 , 2 , 2 , 1 ]
4 // x2 [ n ] =[ 1 , 2 , 3 , 1 ]
5 clear all;
6 clc ;
7 close ;
8 x1=[1,2,2,1];
9 x2=[1,2,3,1];
11 X1=fft(x1 ,-1);
12 X2=fft(x2 ,-1);
13 Y=X1.*X2;
66
15 y=fft(Y,1);
17 disp(y,y [ n]=);
Scilab code Exa 3.15 Determine Sequence x3
1 // Example 3. 15
2 // Program t o Compute x3 [ n ] where X3 [ k]=X1 [ k ] . X2 [ k ]
3 // x1 [ n ] =[ 1 , 2 , 3 , 4 ]
4 // x2 [ n ] =[ 1 , 1 , 2 , 2 ]
5 clear all;
6 clc ;
7 close ;
8 x1=[1,2,3,4];
9 x2=[1,1,2,2];
11 X1=fft(x1 ,-1);
12 X2=fft(x2 ,-1);
13 X3=X1.*X2;
15 x3=fft(X3 ,1);
16 // Di s pl ay s e que nc e x3 [ n ] i n command window
17 disp(x3,x3 [ n]=);
1 // Example 3. 16
3 // x1 [ n ] =[ 1 , 1 , 2 , 1 ]
4 // x2 [ n ] =[ 1 , 2 , 3 , 4 ]
5 clear all;
67
6 clc ;
7 close ;
8 x1=[1,1,2,1];
9 x2=[1,2,3,4];
11 X1=fft(x1 ,-1);
12 X2=fft(x2 ,-1);
13 X3=X1.*X2;
15 x3=fft(X3 ,1);
16 // Di s pl ay s e que nc e x3 [ n ] i n command window
17 disp(x3,x3 [ n]=);
1 // Example 3. 17
2 // Program t o Compute y [ n ] where Y[ k]=X1 [ k ] . X2 [ k ]
3 // x1 [ n ] =[ 0 , 1 , 2 , 3 , 4 ]
4 // x2 [ n ] =[ 0 , 1 , 0 , 0 , 0 ]
5 clear all;
6 clc ;
7 close ;
8 x1=[0,1,2,3,4];
9 x2=[0,1,0,0,0];
11 X1=fft(x1 ,-1);
12 X2=fft(x2 ,-1);
13 Y=X1.*X2;
15 y=round(fft(Y,1));
17 disp(y,y [ n]=);
68
Scilab code Exa 3.18 Output Response
1 // Example 3. 18
2 // Program t o Compute out put r e s ponc e o f f o l l o wi n g
s e que nc e s
3 //x [ n ] =[ 1 , 2 , 3 , 1 ]
4 //h [ n ] =[ 1 , 1 , 1 ]
5 // ( 1 ) Li ne ar Convol ut i on
6 // ( 2 ) Ci r c ul a r Convol ut i on
7 // ( 3 ) Ci r c ul a r Convol ut i on wi th z e r o paddi ng
8 clear all;
9 clc ;
10 close ;
11 x=[1,2,3,1];
12 h=[1,1 ,1];
13 // ( 1 ) Li ne ar Convol ut i on Computati on
14 ylinear=convol (x,h);
15 // Di s pl ay Li ne ar Convol uted Sequence y [ n ] i n command
window
16 disp(ylinear , y l i n e a r [ n]=);
17 // ( 2 ) Ci r c ul a r Convol ut i on Computati on
18 //Now z e r o paddi ng i n h [ n ] s e que nc e t o make l e ng t h
o f x [ n ] and h [ n ] e qual
19 h1=[h,zeros (1,1)];
20 //Now Per f or mi ng Ci r c ul a r Convol ut i on by DFT method
21 X=fft(x,-1);
22 H=fft(h1 ,-1);
23 Y=X.*H;
24 ycircular=fft(Y,1);
25 // Di s pl ay Ci r c ul a r Convol uted Sequence y [ n ] i n
command window
26 disp(ycircular , y c i r c u l a r [ n]=);
27 // ( 3 ) Ci r c ul a r Convol ut i on Computati on wi th z e r o
Paddi ng
28 x2=[x,zeros (1,2)];
29 h2=[h,zeros (1,3)];
30 //Now Per f or mi ng Ci r c ul a r Convol ut i on by DFT method
31 X2=fft(x2 ,-1);
69
32 H2=fft(h2 ,-1);
33 Y2=X2.*H2;
34 ycircularp=fft(Y2 ,1);
35 // Di s pl ay Ci r c ul a r Convol uted Sequence wi th z e r o
Paddi ng y [ n ] i n command window
36 disp(ycircularp , y c i r c u l a r p [ n]=);
Scilab code Exa 3.20 Output Response
1 // Example 3. 20
2 // Program t o Compute Li ne ar Convol ut i on o f f o l l o wi n g
s e que nc e s
3 //x [ n ] =[ 3 , 1 , 0 , 1 , 3 , 2 , 0 , 1 , 2 , 1]
4 //h [ n ] =[ 1 , 1 , 1 ]
5 clear all;
6 clc ;
7 close ;
8 x=[3,-1,0,1,3,2,0,1,2,1];
9 h=[1,1 ,1];
10 // Li ne ar Convol ut i on Computati on
11 y=convol (x,h);
12 // Di s pl ay Sequence y [ n ] i n command window
13 disp(y,y [ n]=);
Scilab code Exa 3.21 Linear Convolution
1 // Example 3. 21
2 // Program t o Compute Li ne ar Convol ut i on o f f o l l o wi n g
s e que nc e s
3 //x [ n] =[ 1 , 2 , 1 , 2 , 3 , 2 , 3 , 1 , 1 , 1 , 2 , 1]
4 //h [ n ] =[ 1 , 2 ]
5 clear all;
6 clc ;
70
7 close ;
8 x=[1,2,-1,2,3,-2,-3,-1,1,1,2,-1];
9 h=[1 ,2];
10 // Li ne ar Convol ut i on Computati on
11 y=convol (x,h);
13 disp(y,y [ n]=);
Scilab code Exa 3.23.a N Point DFT Computation
1 // Example 3. 23 ( a )
3 //N poi nt DFT o f de l t a ( n)
4 clc;
5 syms n k N;
6 x=1;
7 X=symsum(x*exp(-%i*2* %pi*n*k/N),n,0,0);
8 disp(X, X( k )= );
Scilab code Exa 3.23.b N Point DFT Computation
1 // Example 3. 23 ( b)
3 //N poi nt DFT o f de l t a ( nno )
4 clc;
5 syms n k N no;
6 x=1;
7 X=symsum(x*exp(-%i*2* %pi*n*k/N),n,-no ,-no);
8 disp(X, X( k )= );
71
Scilab code Exa 3.23.c N Point DFT Computation
1 // Example 3. 23 ( c )
3 //N poi nt DFT o f an
4 clc;
5 syms a n k N;
6 x=a^n;
7 X=symsum(x*exp(-%i*2* %pi*n*k/N),n,0,N-1);
8 disp(X, X( k )= );
Scilab code Exa 3.23.d N Point DFT Computation
1 // Example 3. 23 ( d)
3 //N poi nt DFT o f x ( n) =1, 0<=n<=N/21
4 clc;
5 syms n k N;
6 x=1;
7 X=symsum(x*exp(-%i*2* %pi*n*k/N),n,0,(N/2) -1);
8 disp(X, X( k )= );
Scilab code Exa 3.23.e N Point DFT Computation
1 // Example 3. 23 ( e )
3 //N poi nt DFT o f x ( n)=exp ( %i 2%pi kon/N) ;
4 clc;
5 syms n k N ko;
6 x=exp(%i*2* %pi*ko*n/N);
7 X=symsum(x*exp(-%i*2* %pi*n*k/N),n,0,(N/2) -1);
8 disp(X, X( k )= );
72
Scilab code Exa 3.23.f N Point DFT Computation
1 // Example 3. 23 ( f )
3 //N poi nt DFT o f x ( n) =1, f o r n=even and 0 , f o r n=odd
4 clc;
5 syms n k N ;
6 x=1; //x ( 2 n) =1, f o r a l l n
7 X=symsum(x*exp(-%i*4* %pi*n*k/N),n,0,N/2-1);
8 disp(X, X( k )= );
Scilab code Exa 3.24 DFT of the Sequence
1 // Example 3. 24
2 // Program t o Compute t he DFT o f t he Sequence x [ n
]=( 1) n , f o r N=4
3 clear all;
4 clc ;
5 close ;
6 N=4;
7 n=0:1:N-1;
8 x=(-1)^n;
10 X = fft (x,-1);
11 // Di s pl ay Sequence X[ k ] i n command window
12 disp(X,X[ k]=);
Scilab code Exa 3.25 8 Point Circular Convolution
73
1 // Example 3. 25
2 // Program t o Compute t he 8poi nt Ci r c ul a r
Convol ut i on o f t he Sequences
3 // x1 [ n ] =[ 1 , 1 , 1 , 1 , 0 , 0 , 0 , 0 ]
4 // x2 [ n]= s i n ( 3 pi n/8)
5 clear all;
6 clc ;
7 close ;
8 x1=[1,1,1,1,0,0,0,0];
9 n=0:1:7;
10 pi =22/7;
11 x2=sin(3*pi*n/8);
13 X1=fft (x1 ,-1);
14 X2=fft (x2 ,-1);
15 // Ci r c ul a r Convol ut i on us i ng DFT
16 Y=X1.*X2;
18 y= fft (Y,1);
20 disp(y,y [ n]=);
Scilab code Exa 3.26 Linear Convolution using DFT
1 // Example 3. 26
2 // Program t o Compute t he Li ne ar Convol ut i on o f t he
f o l l o wi n g Sequences
3 //x [ n ] =[ 1 , 1 , 1]
4 //h [ n ] =[ 2 , 2 , 1 ]
5 clear all;
6 clc ;
7 close ;
8 x=[1,-1,1];
9 h=[2,2,1];
10 // Convol ut i on Computati on
74
11 y= convol(x,h);
13 disp(y,y [ n]=);
Scilab code Exa 3.27.a Circular Convolution Computation
1 // Example 3. 27 ( a )
2 // Program t o Compute t he Convol ut i on o f t he
3 // x1 [ n ] =[ 1 , 1 , 1 ]
4 // x2 [ n ] =[ 2 , 1 , 2]
5 clear all;
6 clc ;
7 close ;
8 x1=[1,1 ,1];
9 x2=[2,-1,2];
11 X1=fft (x1 ,-1);
12 X2=fft (x2 ,-1);
13 Y=X1.*X2;
14 y=fft (Y,1);
16 disp(y,y [ n]=);
Scilab code Exa 3.27.b Circular Convolution Computation
1 // Example 3. 27 ( b)
2 // Program t o Compute t he Convol ut i on o f t he
3 // x1 [ n] =[ 1 , 1 , 1 , 1 , 0]
4 // x2 [ n ] =[ 1 , 0 , 1 , 0 , 1]
5 clear all;
6 clc ;
75
7 close ;
8 x1=[1,1,-1,-1,0];
9 x2=[1,0,-1,0,1];
11 X1=fft (x1 ,-1);
12 X2=fft (x2 ,-1);
13 Y=X1.*X2;
14 y= fft (Y,1);
16 disp(y,y [ n]=);
Scilab code Exa 3.30 Calculate value of N
1 // Example 3. 30
2 // Program t o Ca l c ul a t e N f rom gi ve n data
3 //fm=5000Hz
4 // df =50Hz
5 // t =0. 5 s e c
6 clear all;
7 clc ;
8 close ;
9 fm=5000 //Hz
10 df=50 //Hz
11 t=0.5 // s e c
12 N1=2*fm/df;
13 N=2;
14 while N<=N1, N=N*2,end
15 // Di s pl a yi ng t he val ue o f N i n command window
16 disp(N,N=);
Scilab code Exa 3.32 Sketch Sequence
76
Figure 3.3: Sketch Sequence
1 // Example 3. 32
2 // Program t o pl o t t he r e s u l t o f t he gi ve n s e que nc e
3 //X[ k ] =[ 1 , 2 , 2 , 1 , 0 , 2 , 1 , 2 ]
4 //y [ n]=x [ n / 2 ] f o r n=even , 0 f o r n=odd
5 clear all;
6 clc ;
7 close ;
8 X=[1,2,2,1,0,2,1,2];
9 x = fft (X , 1);
10 y=[x(1) ,0,x(2) ,0,x(3) ,0,x(4) ,0,x(5) ,0,x(6) ,0,x(7) ,0,
x(8) ,0];
11 Y = fft (y , -1);
12 // Di s pl ay s e que nc e Y[ k ] and i n command window
13 disp(Y,Y[ k]=);
14 // Pl o t t i ng t he s e que nc e Y[ k ]
15 k=0:1:15;
16 a = gca ();
19 plot2d3(k,Y,2);
22 xtitle( Pl ot o f Y( k ) , k , Y( k ) );
77
Scilab code Exa 3.36 Determine IDFT
1 // Example 3. 36
2 // Program t o Compute t he IDFT o f t he f o l l o wi n g
Sequence
3 //X[ k] =[ 12 , 1. 5+ j 2 . 598 , 1. 5+ j 0 . 866 , 0 , 1. 5 j 0
. 866 , 1. 5 j 2 . 5 9 8 ]
4 clear all;
5 clc ;
6 close ;
7 j=sqrt(-1);
8 X=[12 , -1.5+j*2.598 , -1.5+j*0.866,0,-1.5-j*0.866, -1.5-
j*2.598];
10 x = fft (X , 1);
11 // Di s pl ay Sequence x [ n ] i n command window
12 disp(round(x),x [ n]=);
78
Chapter 4
THE FAST FOURIER
TRANSFORM
Scilab code Exa 4.3 Shortest Sequence N Computation
1 // Example 4 . 3
2 // Program t o c a l c u l a t e s h o r t e s t s e que nc e N such t hat
al gor i t hm B r uns // f a s t e r than A
3 clear all;
4 clc ;
5 close ;
6 i=0;
7 N=32; // Gi ven
8 // Ca l c ul a t i o n o f Twi ddl e f a c t o r exponent s f o r each
s t a g e
9 while 1==1
10 i=i+1;
11 N=2^i;
12 A=N^2;
13 B=5*N*log2(N);
14 if A>B then break;
15 end;
16 end
17 disp(N, SHORTEST SEQUENCE N = );
79
Scilab code Exa 4.4 Twiddle Factor Exponents Calculation
1 // Example 4 . 4
2 // Program t o c a l c u l a t e Twi ddl e f a c t o r exponent s f o r
each s t a g e
3 clear all;
4 clc ;
5 close ;
6 N=32; // Gi ven
7 // Ca l c ul a t i o n o f Twi ddl e f a c t o r exponent s f o r each
s t a g e
8 for m=1:5
9 disp(m, St age : m = );
10 disp( k = );
11 for t=0:(2^(m-1) -1)
12 k=N*t/2^m;
13 disp(k);
14 end
15 end
Scilab code Exa 4.6 DFT using DIT Algorithm
1 // Example 4 . 6
2 // Program t o f i nd t he DFT o f a Sequence x [ n
] =[ 1 , 2 , 3 , 4 , 4 , 3 , 2 , 1 ]
3 // us i ng DIT Al gor i t hm .
4 clear all;
5 clc ;
6 close ;
7 x = [1,2,3,4,4,3,2,1];
8 //FFT Computati on
80
9 X = fft (x , -1);
10 disp(X, X( z ) = );
Scilab code Exa 4.8 DFT using DIF Algorithm
1 // Example 4 . 8
2 // Program t o f i nd t he DFT o f a Sequence x [ n
] =[ 1 , 2 , 3 , 4 , 4 , 3 , 2 , 1 ]
3 // us i ng DIF Al gor i t hm .
4 clear all;
5 clc ;
6 close ;
7 x = [1,2,3,4,4,3,2,1];
9 X = fft (x , -1);
10 disp(X, X( z ) = );
Scilab code Exa 4.9 8 Point DFT of the Sequence
1 // Example 4 . 9
2 // Program t o f i nd t he 8poi nt DFT o f a Sequence x [ n
] =1 , 0<=n<=7
3 // us i ng DIT, DIF Al gor i t hm .
4 clear all;
5 clc ;
6 close ;
7 x = [1,1,1,1,1,1,1,1];
9 X = fft (x , -1);
10 disp(X, X( z ) = );
81
1 // Example 4. 10
2 // Program t o Compute t he 4poi nt DFT o f a Sequence x
[ n ] =[ 0 , 1 , 2 , 3 ]
3 // us i ng DITDIF Al gor i t hm .
4 clear all;
5 clc ;
6 close ;
7 x = [0,1,2,3];
9 X = fft (x , -1);
10 disp(X, X( z ) = );
Scilab code Exa 4.11 IDFT of the Sequence using DIT Algorithm
1 // Example 4. 11
2 // Program t o Compute t he IDFT o f a Sequence us i ng
DIT Al gor i t hm .
3 //X[ k ] = [ 7 , 0. 707 j 0 . 707 , j , 0. 707 j 0 . 707 , 1 , 0. 707+ j 0
. 7 0 7 , j , 0. 707+ j 0 . 7 0 7 ]
4 clear all;
5 clc ;
6 close ;
7 j=sqrt(-1);
8 X = [7,-0.707-j*0.707,-j,0.707 -j*0.707 ,1 ,0.707+j
*0.707 ,j, -0.707+j*0.707];
9 // I nv e r s e FFT Computati on
10 x = fft (X , 1);
11 disp(x, x ( n) = );
82
1 // Example 4. 12
2 // Program t o Compute t he 8poi nt DFT o f a Sequence
3 //x [ n ] =[ 0 . 5 , 0 . 5 , 0 . 5 , 0 . 5 , 0 , 0 , 0 , 0 ] us i ng r adi x 2 DIT
Al gor i t hm .
4 clear all;
5 clc ;
6 close ;
7 x=[0.5 ,0.5 ,0.5 ,0.5 ,0 ,0 ,0 ,0];
9 X = fft (x , -1);
10 disp(X, X( z ) = );
1 // Example 4. 13
2 // Program t o Compute t he 8poi nt DFT o f a Sequence
3 //x [ n ] =[ 0 . 5 , 0 . 5 , 0 . 5 , 0 . 5 , 0 , 0 , 0 , 0 ] us i ng r adi x 2 DIF
Al gor i t hm .
4 clear all;
5 clc ;
6 close ;
7 x=[0.5 ,0.5 ,0.5 ,0.5 ,0 ,0 ,0 ,0];
9 X = fft (x , -1);
10 disp(X, X( z ) = );
Scilab code Exa 4.14 DFT using DIT Algorithm
1 // Example 4. 14
[ n] =[ 1 , 1 , 1 , 1]
4 clear all;
83
5 clc ;
6 close ;
7 x=[1,-1,1,-1];
9 X =fft (x , -1);
10 disp(X, X( z ) = );
Scilab code Exa 4.15 DFT using DIF Algorithm
1 // Example 4. 15
[ n ] =[ 1 , 0 , 0 , 1 ]
3 // us i ng DIF Al gor i t hm .
4 clear all;
5 clc ;
6 close ;
7 x=[1,0,0,1];
9 X =fft (x , -1);
10 disp(X, X( z ) = );
Scilab code Exa 4.16.a 8 Point DFT using DIT FFT
1 // Example 4. 16 ( a )
2 // Program t o Eval uat e and Compare t he 8poi nt DFT o f
t he gi ve n Sequence
3 // x1 [ n] =1 , 3<=n<=3 us i ng DITFFT Al gor i t hm .
4 clear all;
5 clc ;
6 close ;
7 x1=[1,1,1,1,0,1,1,1];
9 X1 = fft (x1 , -1);
84
10 disp(X1, X1( k ) = );
Scilab code Exa 4.16.b 8 Point DFT using DIT FFT
1 // Example 4. 16 ( b)
2 // Program t o Eval uat e and Compare t he 8poi nt DFT o f
t he gi ve n Sequence
3 // x2 [ n] =1 , 0<=n<=6 us i ng DITFFT Al gor i t hm .
4 clear all;
5 clc ;
6 close ;
7 x2=[1,1,1,1,1,1,1,0];
9 X2 = fft (x2 , -1);
10 disp(X2, X2( k ) = );
Scilab code Exa 4.17 IDFT using DIF Algorithm
1 // Example 4. 17
2 // Program t o f i nd t he IDFT o f t he Sequence us i ng DIF
Al gor i t hm .
3 //X[ k]= [ 4 , 1 j 2 . 414 , 0 , 1 j 0 . 414 , 0 , 1+ j 0 . 414 , 0 , 1+ j 2
. 4 1 4 ]
4 clear all;
5 clc ;
6 close ;
7 j=sqrt(-1);
8 X= [4,1-j*2.414,0,1-j*0.414 ,0 ,1+j*0.414 ,0 ,1+j
*2.414];
10 x = fft (X , 1);
11 disp(x, x ( n) = );
85
Scilab code Exa 4.18 IDFT using DIT Algorithm
1 // Example 4. 18
2 // Program t o f i nd t he IDFT o f t he Sequence X[ k]=
[10 , 2+j 2 ,2,2j 2 ]
4 clear all;
5 clc ;
6 close ;
7 j=sqrt(-1);
8 X = [10,-2+j*2,-2,-2-j*2];
10 x = fft (X , 1);
11 disp(x, x ( n) = );
Scilab code Exa 4.19 FFT Computation of the Sequence
1 // Example 4. 19
2 // Program t o Compute t he FFT o f gi ve n Sequence x [ n
] =[ 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ] .
3 clear all;
4 clc ;
5 close ;
6 x = [1,0,0,0,0,0,0,0];
8 X = fft (x , -1);
9 disp(X, X( z ) = );
Scilab code Exa 4.20 8 Point DFT by Radix 2 DIT FFT
86
1 // Example 4. 20
2 // Program t o Compute t he 8poi nt DFT o f gi ve n
Sequence
3 //x [ n ] =[ 2 , 2 , 2 , 2 , 1 , 1 , 1 , 1 ] us i ng DIT, r adi x 2,FFT
Al gor i t hm .
4 clear all;
5 clc ;
6 close ;
7 x = [2,2,2,2,1,1,1,1];
9 X = fft (x , -1);
10 disp(X, X( z ) = );
Scilab code Exa 4.21 DFT using DIT FFT Algorithm
1 // Example 4. 21
2 // Program t o Compute t he DFT o f gi ve n Sequence
3 //x [ n] =[ 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1] us i ng DITFFT Al gor i t hm
.
4 clear all;
5 clc ;
6 close ;
7 x = [1,-1,-1,-1,1,1,1,-1];
9 X = fft (x , -1);
10 disp(X, X( z ) = );
Scilab code Exa 4.22 Compute X using DIT FFT
1 // Example 4. 22
3 //x [ n] =2n and N=8 us i ng DITFFT Al gor i t hm .
4 clear all;
87
5 clc ;
6 close ;
7 N=8;
8 n=0:1:N-1;
9 x =2^n;
11 X = fft (x , -1);
12 disp(X, X( z ) = );
Scilab code Exa 4.23 DFT using DIF FFT Algorithm
1 // Example 4. 23
3 //x [ n]=c os ( n pi /2) , and N=4 us i ng DIFFFT Al gor i t hm .
4 clear all;
5 clc ;
6 close ;
7 N=4;
8 pi =22/7;
9 n=0:1:N-1;
10 x =cos(n*pi/2);
12 X = fft (x , -1);
13 disp(X, X( z ) = );
1 // Example 4. 24
2 // Program t o Compute t he 8poi nt DFT o f gi ve n
Sequence
3 //x [ n ] =[ 0 , 1 , 2 , 3 , 4 , 5 , 6 , 7 ] us i ng DIF , r adi x 2,FFT
Al gor i t hm .
4 clear all;
88
5 clc ;
6 close ;
7 x = [0,1,2,3,4,5,6,7];
9 X = fft (x , -1);
10 disp(X, X( z ) = );
89
Chapter 5
INFINITE IMPULSE
RESPONSE FILTERS
Scilab code Exa 5.1 Order of the Filter Determination
1 // Example 5 . 1
2 //To Fi nd out t he or de r o f t he f i l t e r
3 clear all;
4 clc ;
5 close ;
6 ap=1; //db
7 as=30; //db
8 op=200; // rad / s e c
9 os=600; // rad / s e c
10 N=log(sqrt ((10^(0.1* as) -1) /(10^(0.1* ap) -1)))/log(os/
op);
11 disp(ceil(N), Order o f t he f i l t e r , N = );
Scilab code Exa 5.2 Order of Low Pass Butterworth Filter
1 // Example 5 . 2
90
2 //To Fi nd out t he or de r o f a Low Pass But t er wor t h
Fi l t e r
3 clear all;
4 clc ;
5 close ;
6 ap=3; //db
7 as=40; //db
8 fp=500; //Hz
9 fs =1000; //Hz
10 op=2*%pi*fp;
11 os=2*%pi*fs;
op);
Scilab code Exa 5.4 Analog Butterworth Filter Design
1 // Example 5 . 4
2 //To Des i gn an Anal og But t er wor t h Fi l t e r
3 clear all;
4 clc ;
5 close ;
6 ap=2; //db
7 as=10; //db
8 op=20; // rad / s e c
9 os=30; // rad / s e c
op);
12 s=%s;
13 HS=1/((s^2+0.76537*s+1)*(s^2+1.8477*s+1));// Tr ans f e r
Funct i on f o r N=4
14 oc=op /(10^(0.1* ap) -1)^(1/(2* ceil(N)));
15 HS1=horner(HS,s/oc);
16 disp(HS1 , Nor mal i zed Tr ans f e r Functi on , H( s ) = );
91
Scilab code Exa 5.5 Analog Butterworth Filter Design
1 // Example 5 . 5
2 //To Des i gn an Anal og But t er wor t h Fi l t e r
3 clear all;
4 clc ;
5 close ;
6 op=0.2* %pi;
7 os=0.4* %pi;
8 e1=0.9;
9 l1=0.2;
10 epsilon=sqrt (1/(e1^2) -1);
11 lambda=sqrt (1/(l1^2) -1);
12 N=log(lambda/epsilon)/log(os/op);
14 s=%s;
15 HS=1/((s^2+0.76537*s+1)*(s^2+1.8477*s+1));// Tr ans f e r
Funct i on f o r N=4
16 oc=op/epsilon ^(1/ ceil(N));
17 HS1=horner(HS,s/oc);
Scilab code Exa 5.6 Order of Chebyshev Filter
1 // Example 5 . 6
2 //To Fi nd out t he or de r o f t he Fi l t e r us i ng
Chebyshev Approxi mati on
3 clear all;
4 clc ;
5 close ;
6 ap=3; //db
92
7 as=16; //db
8 fp =1000; //Hz
9 fs =2000; //Hz
10 op=2*%pi*fp;
11 os=2*%pi*fs;
12 N=acosh(sqrt ((10^(0.1* as) -1) /(10^(0.1* ap) -1)))/acosh
(os/op);
Scilab code Exa 5.7 Chebyshev Filter Design
1 // Example 5 . 7
2 //To Des i gn an anal og Chebyshev Fi l t e r wi th Gi ven
S p e c i f i c a t i o n s
3 clear all;
4 clc ;
5 close ;
6 os=2;
7 op=1;
8 ap=3; //db
9 as=16; //db
10 e1=1/ sqrt (2);
11 l1=0.1;
12 epsilon=sqrt (1/(e1^2) -1);
13 lambda=sqrt (1/(l1^2) -1);
14 N=acosh(lambda/epsilon)/acosh(os/op);
Scilab code Exa 5.8 Order of Type 1 Low Pass Chebyshev Filter
1 // Example 5 . 8
2 //To Fi nd out t he or de r o f t he po l e s o f t he Type 1
Lowpass Chebyshev Fi l t e r
93
3 clear all;
4 clc ;
5 close ;
6 ap=1; //dB
7 as=40; //dB
8 op =1000* %pi;
9 os =2000* %pi;
(os/op);
Scilab code Exa 5.9 Chebyshev Filter Design
1 // Example 5 . 9
2 //To Des i gn a Chebyshev Fi l t e r wi th Gi ven
S p e c i f i c a t i o n s
3 clear all;
4 clc ;
5 close ;
6 ap=2.5; //db
7 as=30; //db
8 op=20; // rad / s e c
9 os=50; // rad / s e c
(os/op);
Scilab code Exa 5.10 HPF Filter Design with given Specications
1 // Example 5. 10
2 //To Des i gn a H. P. F. wi th gi ve n s p e c i f i c a t i o n s
3 clear all;
4 clc ;
94
5 close ;
6 ap=3; //db
7 as=15; //db
8 op=500; // rad / s e c
9 os =1000; // rad / s e c
op);
12 s=%s;
13 HS=1/((s+1)*(s^2+s+1));// Tr ans f e r Funct i on f o r N=3
14 oc=1000 // rad / s e c
15 HS1=horner(HS,oc/s);
Scilab code Exa 5.11 Impulse Invariant Method Filter Design
1 // Example 5. 11
2 //To Des i gn t he Fi l t e r us i ng I mpul s e I n v a r i e n t
Method
3 clear all;
4 clc ;
5 close ;
6 s=%s;
7 T=1;
8 HS=(2)/(s^2+3*s+2);
9 elts=pfss(HS);
10 disp(elts , Fa c t o r i z e d HS = );
11 //The po l e s comes out t o be at 2 and 1
12 p1=-2;
13 p2=-1;
14 z=%z;
15 HZ=(2/(1 -%e^(p2*T)*z^(-1))) -(2/(1-%e^(p1*T)*z^(-1)))
16 disp(HZ, HZ = );
95
1 // Example 5. 12
Method
4 clear all;
5 clc ;
6 close ;
7 s=%s;
8 HS=1/(s^2+ sqrt (2)*s+1);
9 pp=ilaplace(HS);
10 syms n z;
11 t=1;
12 X= symsum (pp*(z^(-n)),n ,0, %inf );
13 disp(X, Fa c t o r i z e d HS = );
1 // Example 5. 13
3 //To Des i gn t he 3 rd Order But t er wor t h Fi l t e r us i ng
I mpul s e I n v a r i e n t Method
4 clear all;
5 clc ;
6 close ;
7 s=%s;
8 HS=1/((s+1)*(s^2+s+1));
9 pp=ilaplace(HS);// I nv e r s e Lapl ac e
10 syms n z;
11 t=1;
12 X= symsum (pp*(z^(-n)),n ,0, %inf );//Z Transf orm
96
13 disp(X, H( z )= );
1 // Example 5. 15
Method
3 clear all;
4 clc ;
5 close ;
6 s=%s;
7 T=0.2;
8 HS=10/(s^2+7*s+10);
9 elts=pfss(HS);
10 disp(elts , Fa c t o r i z e d HS = );
11 //The po l e s comes out t o be at 5 and 2
12 p1=-5;
13 p2=-2;
14 z=%z;
15 HZ=T*(( -3.33/(1 -%e^(p1*T)*z^(-1)))+(3.33/(1 -%e^(p2*T
)*z^(-1))))
16 disp(HZ, HZ = );
Scilab code Exa 5.16 Bilinear Transformation Method Filter Design
1 // Example 5. 16
2 //To Fi nd out Bi l i n e a r Tr ans f or mat i on o f HS=2/(( s +1)
( s +2) )
3 clear all;
4 clc ;
5 close ;
6 s=%s;
7 z=%z;
97
8 HS=2/((s+1)*(s+2));
9 T=1;
10 HZ=horner(HS ,(2/T)*(z-1)/(z+1));
11 disp(HZ, H( z ) = );
Scilab code Exa 5.17 HPF Design using Bilinear Transform
1 // Example 5. 17
2 //To Des i gn an H. P. F. monotoni c i n passband us i ng
Bi l i n e a r Transf orm
3 clear all;
4 clc ;
5 close ;
6 ap=3; //db
7 as=10; //db
8 fp =1000; //Hz
9 fs=350; //Hz
10 f=5000;
11 T=1/f;
12 wp=2*%pi*fp;
13 ws=2*%pi*fs;
14 op=2/T*tan(wp*T/2);
15 os=2/T*tan(ws*T/2);
16 N=log(sqrt ((10^(0.1* as) -1) /(10^(0.1* ap) -1)))/log(op/
os);
18 s=%s;
19 HS=1/(s+1) // Tr ans f e r Funct i on f o r N=1
20 oc=op// rad / s e c
21 HS1=horner(HS,oc/s);
23 z=%z;
24 HZ=horner(HS ,(2/T)*(z-1)/(z+1));
25 disp(HZ, H( z ) = );
98
Scilab code Exa 5.18 Bilinear Transformation Method Filter Design
1 // Example 5. 18
2 //To Fi nd out Bi l i n e a r Tr ans f or mat i on o f H( s ) =( s
2+4. 525) /( s 2+0. 692 s +0. 504)
3 clear all;
4 clc ;
5 close ;
6 s=%s;
7 z=%z;
8 HS=(s^2+4.525) /(s^2+0.692*s+0.504);
9 T=1;
10 HZ=horner(HS ,(2/T)*(z-1)/(z+1));
11 disp(HZ, H( z ) = );
Scilab code Exa 5.19 Single Pole LPF into BPF Conversion
1 // Example 5. 19
2 //To Convert a s i n g l e Pol e LPF i nt o BPF
3 clear all;
4 clc ;
5 close ;
6 s=%s;
7 z=%z;
8 HZ =(0.5*(1+z^(-1)))/(1 -0.302*z^(-2));
9 T=1;
10 wu=3*%pi /4;
11 wl=%pi/4;
12 wp=%pi/6;
13 k=tan(wp/2)/tan((wu-wl)/2);
14 a=cos((wu+wl)/2)/cos((wu -wl)/2);
99
15 transf = -((((k-1)/(k+1))*(z^(-2))) -((2*a*k/(k+1))*(z
^(-1)))+1)/(z^(-2) -(2*a*k/(1+k)*z^(-1))+((k-1)/(k
+1)));
16 HZ1=horner(HZ,transf);
17 disp(HZ1 , H( z ) o f B. P. F = );
Scilab code Exa 5.29 Pole Zero IIR Filter into Lattice Ladder Structure
1 // Example 5. 29
2 // Program t o c onve r t gi ve n I I R pol e z e r o Fi l t e r i nt o
La t t i c e Ladder St r uc t ur e .
3 clear all;
4 clc ;
5 close ;
6 U=1; // Zero Adj us t
7 a(3+U,0+U)=1;
8 a(3+U,1+U)=13/24;
9 a(3+U,2+U)=5/8;
10 a(3+U,3+U)=1/3;
11 a(2+U,0+U)=1; // a (m, 0 ) =1
12 a(2+U,3+U)=1/3;
13 m=3,k=1;
14 a(m-1+U,k+U)=(a(m+U,k+U)-a(m+U,m+U)*a(m+U,m-k+U))
/(1-a(m+U,m+U)*a(m+U,m+U));
15 m=3,k=2;
/(1-a(m+U,m+U)*a(m+U,m+U));
17 m=2,k=1;
/(1-a(m+U,m+U)*a(m+U,m+U));
19 disp( LATTICE COEFFICIENTS );
20 disp(a(1+U,1+U), k1 );
21 disp(a(2+U,2+U), k2 );
22 disp(a(3+U,3+U), k3 );
23 b0=1;
100
24 b1=2;
25 b2=2;
26 b3=1;
27 c3=b3;
28 c2=b2-c3*a(3+U,1+U);
29 c1=b1 -(c2*a(2+U,1+U)+c3*a(3+U,2+U));
30 c0=b0 -(c1*a(1+U,1+U)+c2*a(2+U,2+U)+c3*a(3+U,3+U));
31 disp( LADDER COEFFICIENTS );
32 disp(c0, c0 = );
33 disp(c1, c1 = );
34 disp(c2, c2 = );
35 disp(c3, c3 = );
101
Chapter 6
FINITE IMPULSE
RESPONSE FILTERS
Scilab code Exa 6.1 Group Delay and Phase Delay
1 // Example 6 . 1
3 // Program t o Ca l c ul a t e Group Del ay and Phase Del ay
4 //y ( n) =0. 25 x ( n)+x ( n1) +0. 25 x ( n2)
5 clear all;
6 clc ;
7 close ;
8 //w=pol y ( 0 , w) ;
9 syms w;
10 theeta=-w;
11 gd= -diff (theeta ,w); //Group Del ay
12 pd=-theeta/w; // Phase Del ay
13 disp(gd, GROUP DELAY = );
14 disp(pd, PHASE DELAY = );
102
Figure 6.1: LPF Magnitude Response
103
Scilab code Exa 6.5 LPF Magnitude Response
1 // Example 6 . 5
2 // Program t o Pl ot Magni tude Responce o f i d e a l L. P. F.
wi th wc=0. 5 pi
3 //N=11
4 clear all;
5 clc ;
6 close ;
7 N=11;
8 U=6;
9 for n=-5+U:1:5+U
10 if n==6
11 hd(n)=0.5;
12 else
13 hd(n)=(sin(%pi*(n-U)/2))/(%pi*(n-U));
14 end
15 end
16 [hzm ,fr ]= frmag (hd ,256) ;
17 hzm_dB = 20* log10 (hzm)./ max ( hzm );
18 figure;
19 plot (2*fr , hzm_dB );
20 a= gca ();
21 xlabel ( Frequency w pi );
22 ylabel ( Magni tude i n dB );
23 title ( Frequency Response o f FIR LPF wi th N=11 );
24 xgrid (2);
Scilab code Exa 6.6 HPF Magnitude Response
1 // Example 6 . 6
2 // Program t o Pl ot Magni tude Responce o f i d e a l H. P. F.
wi th wc =0. 25 pi
104
Figure 6.2: HPF Magnitude Response
105
3 //N=11
4 clear all;
5 clc ;
6 close ;
7 N=11;
8 U=6;
9 for n=-5+U:1:5+U
10 if n==6
11 hd(n)=0.5;
12 else
13 hd(n)=(sin(%pi*(n-U))-sin(%pi*(n-U)/4))/(%pi*(n-U));
14 end
15 end
16 [hzm ,fr ]= frmag (hd ,256) ;
18 figure
19 plot (2*fr , hzm_dB )
20 a= gca ();
23 title ( Frequency Response o f FIR HPF wi th N=11 );
24 xgrid (2);
Scilab code Exa 6.7 BPF Magnitude Response
1 // Example 6 . 7
2 // Program t o Pl ot Magni tude Responce o f i d e a l B. P. F.
wi th
3 //wc1 =0. 25 pi and wc2 =0. 75 pi
4 //N=11
5 clear all;
6 clc ;
7 close ;
106
Figure 6.3: BPF Magnitude Response
107
8 N=11;
9 U=6;
10 for n=-5+U:1:5+U
11 if n==6
12 hd(n)=0.5;
13 else
14 hd(n)=(sin(%pi *3*(n-U)/4)-sin(%pi*(n-U)/4))/(%pi*(n-
U));
15 end
16 end
17 [hzm ,fr ]= frmag (hd ,256) ;
19 figure
21 a= gca ();
24 title ( Frequency Response o f FIR BPF wi th N=11 );
25 xgrid (2);
Scilab code Exa 6.8 BRF Magnitude Response
1 // Example 6 . 8
2 // Program t o Pl ot Magni tude Responce o f i d e a l B. R. F.
wi th
3 //wc1 =0. 33 pi and wc2 =0. 67 pi
4 //N=11
5 clear all;
6 clc ;
7 close ;
8 N=11;
9 U=6;
10 for n=-5+U:1:5+U
108
Figure 6.4: BRF Magnitude Response
109
11 if n==6
12 hd(n)=0.5;
13 else
14 hd(n)=(sin(%pi*(n-U))+sin(%pi*(n-U)/3)-sin(%pi *2*(n-
U)/3))/(%pi*(n-U));
15 end
16 end
17 [hzm ,fr ]= frmag (hd ,256) ;
19 figure
21 a= gca ();
24 title ( Frequency Response o f FIR BRF wi th N=11 );
25 xgrid (2);
Scilab code Exa 6.9.a HPF Magnitude Response using Hanning Window
1 // Example 6 . 9 a
3 // us i ng Hanni ng Window
4 //wc1 =0. 25 pi
5 //N=11
6 clear all;
7 clc ;
8 close ;
9 N=11;
10 U=6;
11 h_hann=window( hn ,N);
12 for n=-5+U:1:5+U
13 if n==6
14 hd(n)=0.75;
110
Figure 6.5: HPF Magnitude Response using Hanning Window
111
15 else
17 end
18 h(n)=h_hann(n)*hd(n);
19 end
20 [hzm ,fr ]= frmag (h ,256) ;
22 figure
24 a= gca ();
27 title ( Frequency Response o f FIR HRF wi th N=11
us i ng Hanni ng Window );
28 xgrid (2);
Scilab code Exa 6.9.b HPF Magnitude Response using Hamming Win-
dow
1 // Example 6 . 9 b
3 // us i ng Hamming Window
4 //wc1 =0. 25 pi
5 //N=11
6 clear all;
7 clc ;
8 close ;
9 N=11;
10 U=6;
11 h_hamm=window( hm ,N);
12 for n=-5+U:1:5+U
13 if n==6
14 hd(n)=0.75;
112
Figure 6.6: HPF Magnitude Response using Hamming Window
113
15 else
17 end
18 h(n)=h_hamm(n)*hd(n);
19 end
20 [hzm ,fr ]= frmag (h ,256) ;
22 figure
24 a= gca ();
27 title ( Frequency Response o f FIR HRF wi th N=11
us i ng Hamming Window );
28 xgrid (2);
Scilab code Exa 6.10 Hanning Window Filter Design
1 // Example 6. 10
2 // Program t o Pl ot Magni tude Responce o f gi ve n L. P. F.
wi th s p e c i f i c a t i o n s :
3 //N=7,w=pi /4
4 // Usi ng Hanni ng Window
5 clear all;
6 clc ;
7 close ;
8 N=7;
9 alpha =3;
10 U=1;
11 h_hann=window( hn ,N);
12 for n=0+U:1:6+U
13 if n==4
14 hd(n)=0.25;
114
Figure 6.7: Hanning Window Filter Design
115
15 else
16 hd(n)=(sin(%pi*(n-U-alpha)/4))/(%pi*(n-U-alpha));
17 end
18 h(n)=hd(n)*h_hann(n);
19 end
20 [hzm ,fr ]= frmag (h ,256) ;
22 figure
24 a= gca ();
27 title ( Frequency Response o f gi ve n LPF wi th N=7 );
28 xgrid (2);
Scilab code Exa 6.11 LPF Filter Design using Kaiser Window
1 // Example 6. 11
3 //wp=20rad / s ec , ws=30rad / s ec , wsf =100 rad / s e c
4 // as =44. 0dB, ap=0. 1dB
5 // Usi ng Kai s e r Window
6 clear all;
7 clc ;
8 close ;
9 wsf =100 // rad / s e c
10 ws=30; // rad / s e c
11 wp=20; // rad / s e c
12 as=44.0 //dB
13 ap=0.1 //dB
14 B=ws -wp;
15 wc =0.5*( ws+wp);
116
Figure 6.8: LPF Filter Design using Kaiser Window
117
16 wc1=wc*2*%pi/wsf;
17 delta1 =10^( -0.05* as);
18 delta2 =(10^(0.05* as) -1) /(10^(0.05* as)+1);
19 delta=min(delta1 ,delta2);
20 alphas =-20* log10(delta);
21 alpha =0.5842*( alphas -21) ^0.4+0.07886*( alphas -21)
22 D=(alphas -7.95) /14.36;
23 N1=wsf*D/B+1;
24 N=ceil(N1);
25 U=ceil(N/2);
26 win_l=window( kr ,N,alpha);
27 for n=-floor(N/2)+U:1: floor(N/2)+U
28 if n==ceil(N/2);
29 hd(n)=0.5;
30 else
31 hd(n)=(sin(%pi*(n-U)/2))/(%pi*(n-U));
32 end
33 h(n)=hd(n)*win_l(n);
34 end
35 [hzm ,fr ]= frmag (h ,256) ;
37 figure
39 a= gca ();
42 title ( Frequency Response o f gi ve n LPF us i ng Kai s e r
Window );
43 xgrid (2);
44 disp(h, Fi l t e r Co e f f i c i e n t s , h( n)=);
Scilab code Exa 6.12 BPF Filter Design using Kaiser Window
118
Figure 6.9: BPF Filter Design using Kaiser Window
119
1 // Example 6. 12
2 // Program t o Pl ot Magni tude Responce o f gi ve n B. P. F.
3 //wp1=40 pi rad / s ec , wp2=60 pi rad / s e c
4 // ws1=20 pi rad / s ec , ws2=80 pi rad / s e c
5 // as =30dB, ap=0. 5dB
6 //F=100 Hz
7 // Usi ng Kai s e r Window
8 clear all;
9 clc ;
10 close ;
11 wsf =200* %pi;// rad / s e c
12 ws1 =20* %pi;// rad / s e c
13 ws2 =80* %pi;// rad / s e c
14 wp1 =40* %pi;// rad / s e c
15 wp2 =60* %pi;// rad / s e c
16 as=30 //dB
17 ap=0.5 //dB
18 B=min(wp1 -ws1 ,ws2 -wp2);
19 wc1=wp1 -B/2;
20 wc2=wp2+B/2;
21 wc1=wc1*2*%pi/wsf;
22 wc2=wc2*2*%pi/wsf;
23 delta1 =10^( -0.05* as);
24 delta2 =(10^(0.05* as) -1) /(10^(0.05* as)+1);
25 delta=min(delta1 ,delta2);
26 alphas =-20* log10(delta);
27 alpha =0.5842*( alphas -21) ^0.4+0.07886*( alphas -21)
28 D=(alphas -7.95) /14.36;
29 N1=wsf*D/B+1;
30 N=ceil(N1);
31 U=ceil(N/2);
32 win_l=window( kr ,N,alpha);
33 for n=-floor(N/2)+U:1: floor(N/2)+U
34 if n==ceil(N/2);
35 hd(n)=0.4;
36 else
37 hd(n)=(sin (0.7* %pi*(n-U))-sin (0.3* %pi*(n-U)))/(%pi*(
120
Figure 6.10: Digital Dierentiator using Rectangular Window
n-U));
38 end
39 h(n)=hd(n)*win_l(n);
40 end
41 [hzm ,fr ]= frmag (h ,256) ;
43 figure
45 a= gca ();
48 title ( Frequency Response o f gi ve n LPF us i ng Kai s e r
Window );
49 xgrid (2);
121
Scilab code Exa 6.13.a Digital Dierentiator using Rectangular Window
1 // Example 6. 13 a
2 // Program t o Pl ot Magni tude Responce o f i d e a l
d i f f e r e n t i a t o r wi th s p e c i f i c a t i o n s :
3 //N=8,w=pi
4 // us i ng Re c t angul ar window
5 clear all;
6 clc ;
7 close ;
8 N=8;
9 alpha =7/2;
10 U=1;
11 h_Rect=window( r e ,N);
12 for n=0+U:1:7+U
13 hd(n)=-(sin(%pi*(n-U-alpha)))/(%pi*(n-U-alpha)*(n-U-
alpha));
14 h(n)=hd(n)*h_Rect(n);
15 end
16 [hzm ,fr ]= frmag (h ,256) ;
18 figure
20 a= gca ();
23 title ( Frequency Response o f gi ve n i d e a l
d i f f e r e n t i a t o r us i ng Re c t angul ar Window , N=8 );
24 xgrid (2)
122
Figure 6.11: Digital Dierentiator using Hamming Window
123
Scilab code Exa 6.13.b Digital Dierentiator using Hamming Window
1 // Example 6. 13 b
d i f f e r e n t i a t o r wi th s p e c i f i c a t i o n s :
3 //N=8,w=pi
4 // us i ng Hamming window
5 clear all;
6 clc ;
7 close ;
8 N=8;
9 alpha =7/2;
12 for n=0+U:1:7+U
13 hd(n)=-(sin(%pi*(n-U-alpha)))/(%pi*(n-U-alpha)*(n-U-
alpha));
14 h(n)=hd(n)*h_hamm(n);
15 end
16 [hzm ,fr ]= frmag (h ,256) ;
18 figure
20 a= gca ();
23 title ( Frequency Response o f gi ve n i d e a l
d i f f e r e n t i a t o r us i ng Hamming Window , N=8 );
24 xgrid (2)
Scilab code Exa 6.14.a Hilbert Transformer using Rectangular Window
124
Figure 6.12: Hilbert Transformer using Rectangular Window
Hi l be r t Tr ans f or mer
3 // us i ng Re c t angul ar Window
4 //N=11
5 clear all;
6 clc ;
7 close ;
8 N=11;
9 U=6;
11 for n=-5+U:1:5+U
12 if n==6
13 hd(n)=0;
14 else
15 hd(n)=(1-cos(%pi*(n-U)))/(%pi*(n-U));
16 end
18 end
19 [hzm ,fr]= frmag (h,256) ;
125
Figure 6.13: Hilbert Transformer using Blackman Window
20 figure
21 plot (2*fr ,-hzm);
22 a = gca ();
24 ylabel ( H( exp ( j w) ) );
25 title ( Frequency Response o f Hi l be r t Tr ans f or mer
wi th N=11 us i ng Re c t angul ar Window );
26 xgrid (2);
Scilab code Exa 6.14.b Hilbert Transformer using Blackman Window
Hi l be r t Tr ans f or mer
126
3 // us i ng Blackman Window
4 //N=11
5 clear all;
6 clc ;
7 close ;
8 N=11;
9 U=6;
10 for n=-5+U:1:5+U
11 h_balckmann(n) = 0.42+0.5* cos(2*%pi*(n-U)/(N-1))
+0.08* cos(4* %pi*(n-U)/(N-1));
12 if n==6
13 hd(n)=0;
14 else
15 hd(n)=(1-cos(%pi*(n-U)))/(%pi*(n-U));
16 end
17 h(n)=hd(n)*h_balckmann(n);
18 end
19 [hzm ,fr]= frmag (h,256) ;
20 figure
21 plot (2*fr ,-hzm);
22 a = gca ();
24 ylabel ( H( exp ( j w) ) );
25 title ( Frequency Response o f Hi l be r t Tr ans f or mer
wi th N=11 us i ng Blackman Window );
26 xgrid (2);
Scilab code Exa 6.15 Filter Coecients obtained by Sampling
1 // Example 6. 15
2 // Program t o de t e r mi ne f i l t e r c o e f f i c i e n t s obt ai ne d
by s ampl i ng :
3 //N=7,w=pi /2
4 clear all;
5 clc ;
127
6 close ;
7 N=7;
9 for n=0+U:1:N-1+U
10 h(n)=(1+2* cos(2*%pi*(n-U-3)/7))/N
11 end
12 disp(h, Fi l t e r Co e f f i c i e n t s , h( n)=)
Scilab code Exa 6.16 Coecients of Linear phase FIR Filter
1 // Example 6. 16
2 // Program t o de t e r mi ne f i l t e r c o e f f i c i e n t s obt ai ne d
by s ampl i ng :
3 //N=15
4 clear all;
5 clc ;
6 close ;
7 N=15;
9 for n=0:1:N-1
10 h(n+U)=(1+2* cos(2*%pi*(7-n)/N)+2* cos(4*%pi*(7-n)/N)
+2* cos(6*%pi*(7-n)/N))/N;
11 end
Scilab code Exa 6.17 BPF Filter Design using Sampling Method
1 // Example 6. 17
2 // Program t o de s i gn bandpass f i l t e r wi th f o l l o wi n g
s p e c i f i c a t i o n s :
3 //N=7, f c 1 =1000Hz , f c 2 =3000Hz , F=8000Hz
4 clear all;
5 clc ;
128
Figure 6.14: Frequency Sampling Method FIR LPF Filter
6 close ;
7 N=7;
9 for n=0:1:N-1
10 h(n+U)=2*( cos(2*%pi*(3-n)/N)+cos (4*%pi*(3-n)/N))/N;
11 end
Scilab code Exa 6.18.a Frequency Sampling Method FIR LPF Filter
129
2 // Program t o de s i gn L. P. F. f i l t e r wi th f o l l o wi n g
3 //N=15 , wc=pi /4
4 clear all;
5 clc ;
6 close ;
7 N=15;
8 U=1;
9 for n=0+U:1:N-1+U
10 h(n)=(1+ cos (2*%pi*(7-n)/N))/N;
11 end
12 [hzm ,fr ]= frmag (h ,256) ;
14 figure;
16 a= gca ();
20 xgrid (2)
Scilab code Exa 6.18.b Frequency Sampling Method FIR LPF Filter
2 // Program t o de s i gn L. P. F. f i l t e r wi th f o l l o wi n g
3 //N=15 , wc=pi /4
4 clear all;
5 clc ;
6 close ;
7 N=15;
8 U=1;
9 for n=0+U:1:N-1+U
130
Figure 6.15: Frequency Sampling Method FIR LPF Filter
10 h(n)=(1+ cos (2*%pi*(7-n)/N)+cos (4* %pi*(7-n)/N))/N;
11 end
12 [hzm ,fr ]= frmag (h ,256) ;
14 figure;
16 a= gca ();
20 xgrid (2)
Scilab code Exa 6.19 Filter Coecients Determination
131
Figure 6.16: Filter Coecients Determination
132
1 // Example 6. 19
3 //N=13 ,w=pi /6
4 clear all;
5 clc ;
6 close ;
7 alpha =6;
8 U=1;
9 for n=0+U:1:12+U
10 if n==7
11 hd(n)=0.167;
12 else
14 end
15 end
16 [hzm ,fr ]= frmag (hd ,256) ;
18 figure
20 a= gca ();
24 xgrid (2)
25 disp(hd, Fi l t e r Co e f f i c i e n t s , h( n)=);
Scilab code Exa 6.20 Filter Coecients using Hamming Window
1 // Example 6. 20
3 //N=13 ,w=pi /6
133
Figure 6.17: Filter Coecients using Hamming Window
134
4 // Usi ng Hamming Window
5 clear all;
6 clc ;
7 close ;
8 N=13;
9 alpha =6;
10 U=1;
12 for n=0+U:1:12+U
13 if n==7
14 hd(n)=0.167;
15 else
17 end
18 h(n)=hd(n)*h_hamm(n);
19 end
20 [hzm ,fr ]= frmag (h ,256) ;
22 figure
24 a= gca ();
28 xgrid (2)
Scilab code Exa 6.21 LPF Filter using Rectangular Window
1 // Example 6. 21
135
Figure 6.18: LPF Filter using Rectangular Window
3 //N=7, f c =1000Hz , F=5000Hz
4 clear all;
5 clc ;
6 close ;
7 N=7;
8 U=4;
10 for n=-3+U:1:3+U
11 if n==4
12 hd(n)=0.4;
13 else
14 hd(n)=(sin (2*%pi*(n-U)/5))/(%pi*(n-U));
15 end
17 end
18 [hzm ,fr ]= frmag (h ,256) ;
20 figure
136
22 a= gca ();
26 xgrid (2)
Scilab code Exa 6.28 Filter Coecients for Direct Form Structure
1 // Example 6. 28
2 // Program t o c a l c u l a t e FIR Fi l t e r c o e f f i c i e n t s f o r
t he d i r e c t f orm s t r u c t u r e
3 // k1=1/2 , k2=1/3 , k3=1/4
4 clear all;
5 clc ;
6 close ;
7 U=1;
8 k1=1/2;
9 k2=1/3;
10 k3=1/4;
11 a(3+U,0+U)=1;
12 a(1+U,1+U)=k1;
13 a(2+U,2+U)=k2;
14 a(3+U,3+U)=k3;
15 m=2,k=1;
16 a(m+U,k+U)=a(m-1+U,k+U)+a(m+U,m+U)*a(m-1+U,m-k+U);
17 m=3,k=1;
19 m=3,k=2;
21 disp(a(3+U,0+U), a ( 3 , 0 ) );
22 disp(a(3+U,1+U), a ( 3 , 1 ) );
23 disp(a(3+U,2+U), a ( 3 , 2 ) );
24 disp(a(3+U,3+U), a ( 3 , 3 ) );
137
Scilab code Exa 6.29 Lattice Filter Coecients Determination
1 // Example 6. 29
2 // Program t o c a l c u l a t e gi ve n FIR Fi l t e r s La t t i c e
f orm c o e f f i c i e n t s .
3 clear all;
4 clc ;
5 close ;
7 a(3+U,0+U)=1;
8 a(3+U,1+U)=2/5;
9 a(3+U,2+U)=3/4;
10 a(3+U,3+U)=1/3;
11 a(2+U,0+U)=1; // a (m, 0 ) =1
12 a(2+U,3+U)=1/3;
13 m=3,k=1;
/(1-a(m+U,m+U)*a(m+U,m+U));
15 m=3,k=2;
/(1-a(m+U,m+U)*a(m+U,m+U));
17 m=2,k=1;
/(1-a(m+U,m+U)*a(m+U,m+U));
19 disp(a(1+U,1+U), k1 );
20 disp(a(2+U,2+U), k2 );
21 disp(a(3+U,3+U), k3 );
138
Chapter 7
FINITE WORD LENGTH
EFFECTS IN DIGITAL
FILTERS
Scilab code Exa 7.2 Subtraction Computation
1 // Example 7 . 2
2 //To Compute Subt r a c t i o n
3 // ( a ) 0. 25 f rom 0 . 5
4 clear all;
5 clc ;
6 close ;
7 a=0.5;
8 b=0.25;
9 c=a-b;
10 disp(c, = ,b, ,a, PART 1 );
11 // ( a ) 0 . 5 f rom 0. 25
12 d=b-a;
13 disp(d, = ,a, ,b, PART 2 );
139
Scilab code Exa 7.14 Variance of Output due to AD Conversion Process
1 // Example 7. 14
2 //To Compare t he Var i e nc e o f Output due t o A/D
Conver s i on pr o c e s s
3 //y ( n) =0. 8y ( n1)+x ( n)
4 clear all;
5 clc ;
6 close ;
7 n=8; // Bi t s
8 r=100; // Range
9 Q=2*r/(2^n); // Quant i z at i on Step Si z e
10 Ve=(Q^2) /12;
11 Vo=Ve *(1/(1 -0.8^2));
12 disp(Q, QUANTIZATION STEP SIZE = );
13 disp(Ve, VARIANCE OF ERROR SIGNAL = );
14 disp(Vo, VARIANCE OF OUTPUT = );
140
Chapter 8
MULTIRATE SIGNAL
PROCESSING
Scilab code Exa 8.9 Two Component Decomposition
1 // Example 8 . 9
3 // Devel op a two component de c ompos i t i on f o r t he
t r a n s f e r f unc t i o n
4 // and de t e r mi ne P0( z ) and P1( z )
5 clear all;
6 clc ;
7 close ;
8 syms z a n;
9 HZ=(z)/(z-a);
10 hn=a^n;// I nv e r s e Z Transf orm o f HZ
11 h2n=a^(2*n);
12 P0=symsum(h2n*z^(-n),n,0,%inf);
13 h2n1=a^(2*n+1);
14 P1=symsum(h2n1*z^(-n),n,0,%inf);
15 disp(P0, P0( Z) = );
16 disp(P1, P1( Z) = );
141
Scilab code Exa 8.10 Two Band Polyphase Decomposition
1 // Example 8. 10
2 // Devel op a two band pol yphas e de c ompos i t i on f o r t he
t r a n s f e r f unc t i o n
3 //H( z )=z2+z+2/z 2+0. 8 z +0. 6
4 clear all;
5 clc ;
6 close ;
7 z=%z;
8 HZ=(z^2+z+2)/(z^2+0.8*z+0.6);
9 HZa=horner(HZ,-z);
10 P0 =0.5*( HZ+HZa);
11 P1 =0.5*(HZ -HZa);
12 disp(P1/z, + ,P0 , H( z ) = )
142
Chapter 9
STATISTICAL DIGITAL
SIGNAL PROCESSING
Scilab code Exa 9.7.a Frequency Resolution Determination
1 // Example 9 . 7 ( a )
2 // Program To Determi ne Frequency Re s o l ut i o n o f
Ba r t l e t t ,
3 // Welch ( 50% Over l ap ) and BlackmannTukey Methods
4 clear all;
5 clc;
6 close;
7 // Data
8 Q=10; // Qual i t y Fact or
9 N=1000; // Sampl es
10 //FREQUENCY RESOLUTION CALCULATION
11 K=Q;
12 rb =0.89*(2* %pi*K/N);
13 rw =1.28*(2* %pi*9*Q)/(16*N);
14 rbt =0.64*(2* %pi *2*Q)/(3*N);
16 disp(rb, Re s o l ut i o n o f Ba r t l e t t Method);
17 disp(rw, Re s o l ut i o n o f Welch ( 50% o ve r l a p ) Method);
18 disp(rbt , Re s o l ut i o n o f BlackmannTukey Method);
143
Scilab code Exa 9.7.b Record Length Determination
1 // Example 9 . 7 ( b)
2 // Program To Determi ne Record Length o f Ba r t l e t t ,
3 // Welch ( 50% Over l ap ) and BlackmannTukey Methods
4 clear all;
5 clc;
6 close;
7 // Data
8 Q=10; // Qual i t y Fact or
9 N=1000; // Sampl es
10 //RECORD LENGTH CALCULATION
11 lb=N/Q;
12 lw=16*N/(9*Q);
13 lbt =3*N/(2*Q);
15 disp(lb, Record Length o f Ba r t l e t t Method);
16 disp(lw, Record Length o f Welch ( 50% o ve r l a p ) Method
);
17 disp(lbt , Record Length o f BlackmannTukey Method);
Scilab code Exa 9.8.a Smallest Record Length Computation
1 // Example 9 . 8 ( a )
2 // Program To Determi ne Sma l l e s t Record Length o f
Ba r t l e t t Method
3 clear all;
4 clc;
5 close;
6 // Data
7 fr =0.01; // Frequency Re s o l ut i o n
144
8 N=2400; // Sampl es
9 //RECORD LENGTH CALCULATION
10 lb =0.89/ fr;
12 disp(lb, Record Length o f Ba r t l e t t Method);
Scilab code Exa 9.8.b Quality Factor Computation
1 // Example 9 . 8 ( b)
2 // Program To Determi ne Qual i t y Fact or o f Ba r t l e t t
Method
3 clear all;
4 clc;
5 close;
6 // Data
7 fr =0.01; // Frequency Re s o l ut i o n
8 N=2400; // Sampl es
9 lb =0.89/ fr;
10 //QUALITY FACTOR CALCULATION
11 Q=N/lb;
13 disp(Q, Qual i t y Fact or o f Ba r t l e t t Method);
145
Chapter 11
DIGITAL SIGNAL
PROCESSORS
Scilab code Exa 11.3 Program for Integer Multiplication
1 // Program 11. 3
2 // Program To Ca l c ul a t e t he val ue o f t he f unc t i o n
3 //Y=AB
4 clear all;
5 clc;
6 close;
7 // I nput Data
8 A=input( Enter I nt e g e r Number A = );
9 B=input( Enter I nt e g e r Number B = );
10 // Mu l t i p l i c a t i o n Computati on
11 Y=A*B;
13 disp(Y,Y = AB = );
Scilab code Exa 11.5 Function Value Calculation
146
1 // Program 11. 5
2 // Program To Ca l c ul a t e t he val ue o f t he f unc t i o n
3 //Y=AX1+BX2+CX3
4 clear all;
5 clc;
6 close;
7 // Data
8 A=1;
9 B=2;
10 C=3;
11 X1=4;
12 X2=5;
13 X3=6;
14 //Compute t he f unc t i o n
15 Y=A*X1+B*X2+C*X3;
17 disp(Y,Y = AX1+BX2+CX3 = );
147

Digital Signal Processing - R. Babu

Uploaded by

Document Informationclick to expand document information

Copyright:

Available Formats

Digital Signal Processing - R. Babu

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Digital Signal Processing - R. Babu

Uploaded by

Copyright:

Available Formats

Scilab Textbook Companion for

You might also like