K. S. SCHOOL OF ENGINEERING & MANAGEMENT # 15, Mallasandra Mallasandra,, Off Kanakapur Kanakapura a Road, Bangal Bangalore ore-56 -56006 0062, 2, Karnat Karnataka aka,, India. India.
DEPARTMENT OF ELECTRONICS & COMMUNICATION ENGINEERING
Digit igital al Sig Signa nall Pro Proces cessing ing Lab Manu Manua al Sub Code: 10ECL57 Sem:V
Prepa repare red d By Mr. Ravikiran B. A., Asst. Professor Mrs. Vidhya R., Asst. Professor
Table of Contents PART – A
1
Verification of Sampling theorem.
1
2
Impulse response of a given system
4
3
Linear convolution of two given sequences.
7
4
Circular co convolution of two given sequences
11
5
Autocorrelation of a given sequence and verification of its properties.
15
6
Cross-correlation of given sequences and verification of its properties.
18
7
Solving a given difference equation.
21
8
Computatio Computation n of N point DFT of a given sequence sequence and to plot magnitude and phase spectrum.
23
9
Linear convolution of two sequences using DFT and IDFT.
26
10
Circular convolution of two given sequences using DFT and IDFT
29
11
Design Design and implementat implementation ion of FIR filter to meet given given specifications.
32
12
Design and implementation of IIR filter to meet given specifications.
36
PART – B
About the DSP Trainer Kit
44
Using Code Composer Studio
47
1
Linear convolution of two given sequences.
56
2
Circular co convolution of two given sequences
58
3
Computation of N point DFT of a given sequence.
60
4
Impulse Response of the First Order and Second Order Systems
62
Viva Questions
64
th
DSP Laboratory (10ECL57)
5 Sem
KSSEM, Bangalore
PROGRAM 1 VERIFICATION VERIFICATION OF SAMPLING THEOREM Aim: To write the MATLAB code for verifying Sampling Theorem. Generate Generate a sinuso sinusoidal idal wave of 1kHz. 1kHz. Calcu Calculate late the Nyquist Nyquist frequency frequency,, and and verify verify Sampling Theorem, showing output waveforms for undersampled, oversampled and right sampled cases.
Theory: Sampling is the process process of converting an continuous time signal signal into a discrete time signal. In sampling, the values values of the continuous time signal signal is recorded at at discrete intervals of of time (usually equidistant). equidistant). The number of samples samples taken during one second is called called the sampling rate. Sampling is described by the relation:
( ) =
(
)
−∞ <
<∞
Where ( ) is the discrete-time signal obtained by sampling the analog signal every T = 1/ is known as the Sampling Frequency. seconds. The Sampling Theorem states that : “A bandlimited signal can be reconstructed exactly if it is sampled at a rate atleast twice the maximum frequency component in it."
( ) = sin( ) = sin(2 ) with maximum Assume a band-limited signal frequency component ′ ′. The theorem says that, for a good reconstruction of the original continuous time signal, the sampling sampling frequency must be at least 2 . This frequency is known
as the “Nyquist Rate”.
Samplin Sampling g this this sign signal al at at
gives gives us the the discre discrete te time time signa signal: l:
(
)=
sin
2
Now, assuming the sampling frequency is more than the Nyquist Frequency, the continuous time signal can be reconstructed accurately using the interpolation function:
( ) =
sin sin 2 2
Then, approximated approximated recovered recovered signal can be be written as:
( ) =
Dept of ECE
−
Page | 1
th
DSP Laboratory (10ECL57)
5 Sem
KSSEM, Bangalore
Whenever Whenever the Sampling Sampling frequency frequency is greate greaterr than than or equal equal to the Nyqu Nyquist ist Frequ Frequency ency,, the signal can be reconstructed faithfully, capturing all the essential properties of the original < 2 , we encounter a problem called “Aliasing”, continuous-time signal. However, when where distortion is caused by high frequencies frequencies overlapping overlapping low frequencies. frequencies. A lot of data is lost in this t his process and the signal cannot be recovered.
MATLAB CODE: % Experiment 1 : Sampling Theorem Verification clear all all; ; close all all; ; clc; % Signal Parameters f1 = 1000; f2 = 1900; fmax fmax = max(f max(f1, 1,f2 f2); ); T = 1/min(f1,f2); 1/min(f1,f2); t = 0:0.01*T:2*T; 0:0.01*T:2*T;
% % % % %
Signal 1 Frequency = 1kHz Signal 2 Frequency = 1.4 kHz Maximum frequency component of signal Signal Period should cover entire length Time index
% Generate the original signal and plot it: x = cos(2*pi*t*f1)+ cos(2*pi*t*f2); subplot(2,2,1); plot(t,x); grid on on; ; title('Continuous title('Continuous signal'); signal'); xlabel('t' xlabel('t'); ); ylabel('x(t)' ylabel('x(t)'); ); % Oversampling Condition: fs1 = 10*fmax;
% Composite Signal
% Oversampling (fs > 2f)
n1 = 0:1/fs1:2*T; % Time scale x1 = cos(2*pi*f1*n1)+cos(2*pi*f2*n1); % Generating sampled signal subplot(2,2,2); stem(n1,x1); hold on on; ; plot(n1,x1,'r' plot(n1,x1,'r'); ); grid on on; ; hold off off; ; title('Oversampl title('Oversampling ing Condition : Fs = 10F'); xlabel('n' xlabel('n'); ); ylabel('x(n)' ylabel('x(n)'); ); % Right Sampling Condition: fs2 = 2*fmax;
% Nyquist Sampling (fs = 2f)
n2 = 0:1/fs2:2*T; x2 = cos(2*pi*f1*n2)+cos(2*pi*f2*n2); subplot(2,2,3); stem(n2,x2); hold on on; ; plot(n2,x2,'r' plot(n2,x2,'r'); ); grid on on; ; hold off off; ; title('Sampling title('Sampling at Nyquist Frequency : Fs = 2F'); 2F'); xlabel('n' xlabel('n'); ); ylabel('x(n)' ylabel('x(n)'); );
Dept of ECE
Page | 2
DSP Laboratory (10ECL57) % Under Sampling Condition: fs3 = 1.2*fmax;
th
5 Sem
KSSEM, Bangalore
% Undersampling (fs < 2f)
n3 = 0:1/fs3:2*T; 0:1/fs3:2*T; x3 = cos(2*pi*f1*n3)+cos(2*pi*f2*n3); subplot(2,2,4); stem(n3,x3); hold on on; ; plot(n3,x3,'r' plot(n3,x3,'r'); ); grid on on; ; hold off off; ; title('Undersampling title('Undersampling Condition : Fs = 1.2 f'); f'); xlabel('n' xlabel('n'); ); ylabel('x(n)' ylabel('x(n)'); );
OUTPUT:
Dept of ECE
Page | 3
th
DSP Laboratory (10ECL57)
5 Sem
KSSEM, Bangalore
PROGRAM 2 IMPULSE RESPONSE OF A GIVEN SYSTEM Aim: To write the MATLAB code to find the impulse response of a given second-order system whose whose difference equation equation representation is given. given. Assume a second-order system represented by the following difference equation:
( ) =
( )+
( − 1) +
( − 2) +
( − 1) +
( − 2)
Theory: Impulse response of a system is defined as the output of a given system, when the input applied to the system, system, is in the form of of a unit unit impulse, or a Dirac delta function. The impulse response completely completely characterizes characterizes the behaviour of any LTI LTI system. The impulse response response is often determined determined from knowledge knowledge of of the system system configuration configuration and dynamics, or can be measured by applying and approximate impulse to the system input. Discrete-time LTI systems can can also be described using using Difference Difference Equations. A linear constant-coefficient constant-coefficient difference equation can be of the form:
[ − ] =
[ − ]
Where the integer N is termed the order of the difference equation, and corresponds to the maximum memory involving the system output. The order generally represents the number of energy storage devices in a physical system. We can calculate calculate the impulse response of the system using using Z-transforms as shown shown in the following example: Consider a difference equation:
( ) =
( ) + 0.2 ( − 1) 1) − 1.5 ( − 2) 2 ) − 3 ( − 1) 1) + 0.12 ( − 2)
This can be rewritten in the standard form as:
( ) + 3 ( − 1) − 0.12 ( − 2) =
( ) + 0.2 ( − 1 ) − 1.5 ( − 2)
Finding the Z-transform of the equation:
( ) + 3
( ) − 0.12
( ) = ( ) + 0.2
( ) − 1.5
( )
Or:
( )[1 + 3
Dept of ECE
− 0.12
] = ( )[1 + 0.2
− 1.5
]
Page | 4
th
DSP Laboratory (10ECL57)
5 Sem
KSSEM, Bangalore
Transfer Function of the system can be obtained as :
( ) =
( ) [1 + 3 = ( ) [1 + 0.2
− 0.12 − 1. 1.5
] ]
By long divisio division, n, we get:
( ) = 1 − 2.8
+ 7.02
− 21.4
+ 65.03
By taking Inverse-Z transform, we can obtain the Impulse Response as:
ℎ[ ] = [1
− 2.8
7.02
21.4
65.03]
MATLAB CODE: % Experiment Experiment 2 : Impulse Impulse Response Response of a Given Second-Order System clear all all; ; close all all; ; clc; % b a N
Accept Input and Output signal Co-efficients: = input('Enter input('Enter the coefficients of x(n) in 1-D Matrix Form: '); = input('Enter input('Enter the coefficients coefficients of y(n) in 1-D Matrix Form: '); = input('Enter input('Enter the number of samples of impulse response desired: ');
% Calculate Impulse Response using IMPZ function: % [H,T] = IMPZ(B,A,N) computes N samples of the impulse response, using % coefficients B and A from difference equation representation. [h,t] = impz(b,a,N); %Plot and Display impulse response co-efficients: stem(t,h); title('Impulse title('Impulse Response Plot'); Plot'); ylabel('h(n)' ylabel('h(n)'); ); xlabel('n' xlabel('n'); ); disp('Impulse disp('Impulse Response Coefficients:'); Coefficients:'); disp(h);
OUTPUT: Enter Enter the coefficie coefficients nts of x(n) in 1-D 1-D Matrix Form: Form: [1 0.2 0.2 -1.5] Enter the coefficients coefficients of y(n) y(n) in 1-D Matrix Form: Form: [1 3 -0.12] Enter the number of samples of impulse response desired: 5 Impulse Response Coefficients: 1.0000 -2.8000 7.0200 -21.3960 65.0304 Dept of ECE
Page | 5
DSP Laboratory (10ECL57)
Dept of ECE
th
5 Sem
KSSEM, Bangalore
Page | 6
th
DSP Laboratory (10ECL57)
5 Sem
KSSEM, Bangalore
PROGRAM 3 LINEAR CONVOLUTION OF TWO GIVEN SEQUENCES MATLAB code to perform Linear Linear Convolution upon upon two given discrete Aim: To write the MATLAB time signal signals. s. Theory: Convolution is the process process used used to find the response response of a Linear Time Invariant system to a given given inpu input, t, assum assuming ing we alre already ady know know the impuls impulsee respon response se of that that syste system. m. In case case of continuous-time signals, we can find the system response using the Convolution Integral, while in case of discrete-time systems, the response can be calculated using the Convolution Sum.
( ) be two discrete-time signals. The convolution sum of the two signals Let ( ) and can be calculated using the formula: ( ) = 1( ) ∗ 2( ) =
1( ) 2( 2( − )
If 1( ) is a MM- point point sequenc sequencee and and 2( ) is an N – point sequence, then the convolved sequence, ( ) is a (M+N-1) – point sequence. We can perform perform the convolution convolution by different methods: “ CONV” ” function : 1. Using MATLAB’s “CONV MATLAB has a built- in function called “conv” function, which basically performs a linear convolution of any two given sequences. sequences. 2. Using Using the the Line Linear ar Conv Convolu oluti tion on Sum Sum form formula ula : Here Here,, we we use use the the con convo volu luti tion on sum sum for formu mula la and and sub subst stit itut utee val value uess of of and and in the the expression, and calculate the values of the convolved signal. Alternatively, we can perform the signal inversion-time shift-superposition method, by which we can calculate the resultant signal values.
Assume two discrete-time sequences by:
1 and 2 in a Linear Time Invariant System, given
1( ) = { = {1,2,−1,3 1,2,−1,3}} and 2( ) = { = {2,3 2,3,, −2} −2} We see that length of sequence 1 is (M = 4) and that of sequence 2 is (N = 3). Therefore, the length of the convolved sequence will be (M+N-1 = 6). Using any of the above given methods, methods, we see that the resultant convolved convolved sequence can can be given by:
( ) = 1( ) ∗ 2( ) = { 2
Dept of ECE
7
2
− 1 11
− 6}
Page | 7
DSP Laboratory (10ECL57)
th
5 Sem
KSSEM, Bangalore
MATLAB CODE: 1. Using “conv” function: %% Linear Convolution using CONV command clear all all; ; close all all; ; clc; % Accept Accept inpu input t sign signal al seque sequences nces x1 = input('Enter input('Enter Input Sequence for Signal x1(n): '); '); x2 = input('Enter input('Enter Input Sequence for Signal x2(n): ' '); ); %Perform Linear Convolution using CONV command y=conv(x1,x2); %Plot Input and Convolved Signals subplot(3,1,1); stem(x1); title('Input title('Input Signal x1(n)'); x1(n)'); xlabel('n' xlabel('n'); ); ylabel('x1(n)' ylabel('x1(n)'); ); subplot(3,1,2); stem(x2); title('Input title('Input Signal x2(n)'); x2(n)'); xlabel('n' xlabel('n'); ); ylabel('x2(n)' ylabel('x2(n)'); ); subplot(3,1,3); stem(y); title('Convolved title('Convolved Signal y(n) = x1(n)*x2(n)'); x1(n)*x2(n)'); xlabel('n' xlabel('n'); ); ylabel('y(n)' ylabel('y(n)'); ); % Display the convolved Sequence in Command Window disp('Convolved disp('Convolved sequence:'); sequence:'); disp(y);
2. Using Using Conv Convolu olutio tion n Sum Sum formu formula la:: %% Linear Convolution without using CONV command clear all all; ; close all all; ; clc; x1 = input('Ent input('Enter er Input Input Sequenc Sequence e for Signal Signal x1(n x1(n): ): '); n1 = length(x1); x2 = input('Enter input('Enter Input Sequence for Signal x2(n): '); '); n2=length(x2); N = n1+n2-1; T = 1:N;
%Length of Convolved Sequence % Create Time Index
%Zero padding to make sequences of length N x1=[x1 x1=[x1 zeros zeros(1, (1,N-n N-n1)] 1)]; ; x2=[x2 zeros(1,N-n2)]; %Initializing Output sequence of zeros. y = zeros(1,N); %Performing Linear Convolution:
Dept of ECE
Page | 8
DSP Laboratory (10ECL57)
th
5 Sem
KSSEM, Bangalore
for n = 1:N % y(n) = 0R; for k = 1:n y(n)=y(n)+x1(k)*x2(n-k+1); end end % Plot Plot Input Input and and Output Output Sequ Sequence ences: s: subplot(3,1,1); stem(T,x1); title('Input title('Input Signal x1(n)'); x1(n)'); xlabel('n' xlabel('n'); ); ylabel('x1(n)' ylabel('x1(n)'); ); subplot(3,1,2); stem(T,x2); title('Input title('Input Signal x2(n)'); x2(n)'); xlabel('n' xlabel('n'); ); ylabel('x2(n)' ylabel('x2(n)'); ); subplot(3,1,3); stem(T,y); title('Convolved title('Convolved Signal y(n) = x1(n)*x2(n)'); xlabel('n' xlabel('n'); ); ylabel('y(n)' ylabel('y(n)'); ); % Display the convolved Sequence in Command Window disp('Convolved disp('Convolved sequence:'); sequence:'); disp(y);
OUTPUT: Enter Enter Input Input Sequenc Sequence e for Signa Signal l x1(n): x1(n): [1 2 -1 3] Enter Enter Input Input Sequenc Sequence e for Signa Signal l x2(n): x2(n): [2 3 -2] Conv Co nvol olve ved d se sequ quen ence ce: : 2 7 2 -1 11 -6
Dept of ECE
Page | 9
DSP Laboratory (10ECL57)
Dept of ECE
th
5 Sem
KSSEM, Bangalore
Page | 10
th
DSP Laboratory (10ECL57)
5 Sem
KSSEM, Bangalore
PROGRAM 4 CIRCULAR CONVOLUTION OF TWO GIVEN SEQUENCES SEQUENCES MATLAB code to perform Circular Convolution Convolution upon two given discrete discrete Aim: To write the MATLAB time signals. Theory: The Circular convolution, also known known as cyclic convolution, convolution, of two aperiodic functions occurs occurs when one of them is convolved in the normal way with a periodic summation of the other function. Circular convolution is only only defined for finite length functions (usually (usually equal in length), continuous or discrete in time. In circular convolution, it is as if the finite length functions repeat in time, periodically. Because the input functions are now periodic, the convolved output is also periodic. Circular convolution sum can be calculated using the formula:
( ) = 1( ) ∗ 2( ) = = 0,1, … . ,
For
1( ) 2 (
− )
−1
Circular convolution can be performed in different ways : 1. Using the the expression expression for linear linear convolutio convolution n sum, but assuming assuming the the signal signal repeats repeats periodically. This can be done by changing changing the negative negative indices of (n-k) to repetitions repetitions of the latter portions of the original aperiodic signal. 2. Convolutio Convolution n in time domain domain correspon corresponds ds to multiplicati multiplication on in frequency frequency domain. domain. To make use of this property, we can can calculate the DTFT of each of the aperiodic signals, signals, multiply these in the frequency domain, and find the IDFT of the product, to get the periodic convolved signal in time domain. Let us take the case of two discrete-time aperiodic signals given by:
1( ) = {2,1,2,1}
2( ) = {1,2,3,4}
and
Using Using the the form formula ula with with N = 4. For For m = 0:
(0) =
1( ) 2 (− )
= 14
For m = 1:
(1) =
1( ) 2 (1 − )
= 16
For m = 2:
Dept of ECE
Page | 11
th
DSP Laboratory (10ECL57)
5 Sem
KSSEM, Bangalore
(2) =
1( ) 2 (2 − )
= 14
(3) =
1( ) 2 (3 − )
= 16
For m = 3:
( ) = {14,16,14,16}
So, we get the circular convolution sum as: MATLAB CODE:
1. Using Using Conv Convolu olutio tion n Sum Sum Form Formula ula:: %% Circular Convolution using Formula clear all all; ; close all all; ; clc; x1 = input('Enter input('Enter Input Sequence for Signal x1(n): '); '); n1 = length(x1); x2 = input('Ent input('Enter er Input Input Sequenc Sequence e for Sign Signal al x2(n): x2(n): '); n2=length(x2); N = max(n1,n2); T = 1:N;
% Length of Convolved Sequence % Create Time Index
%Zero padding to make sequences of length N x1=[x1 zeros(1,N-n1)]; x2=[x2 zeros(1,N-n2)]; %Initializing Output sequence of zeros. y = zeros(1,N); %Performing Linear Convolution: for m=1:N for n=1:N i=m-n+1; %(m-n+1) since we're taking index from 1 if(i<=0) if (i<=0) i=N+i; end y(m)=y(m)+x1(n)*x2(i); %Convolution Sum Formula end end % Plot Input and Output Sequences: subplot(3,1,1); stem(T,x1); title('Input title('Input Signal x1(n)'); x1(n)'); xlabel('n' xlabel('n'); ); ylabel('x1(n)' ylabel('x1(n)'); ); subplot(3,1,2); stem(T,x2);
Dept of ECE
Page | 12
DSP Laboratory (10ECL57)
th
5 Sem
KSSEM, Bangalore
title('Input title('Input Signal x2(n)'); x2(n)'); xlabel('n' xlabel('n'); ); ylabel('x2(n)' ylabel('x2(n)'); ); subplot(3,1,3); stem(T,y); title('Convolved title('Convolved Signal y(n) = x1(n)*x2(n)'); x1(n)*x2(n)'); xlabel('n' xlabel('n'); ); ylabel('y(n)' ylabel('y(n)'); ); % Display the convolved Sequence in Command Window disp('Convolved disp('Convolved sequence:'); sequence:'); disp(y);
function. 2. Using “cconv” function. %% Circular Circular Conv Convolut olution ion using using CCONV CCONV command command clear all all; ; close all all; ; clc; % Accept input signal sequences x1 = input('Enter input('Enter Input Sequence for Signal x1(n): '); '); x2 = input('Enter input('Enter Input Sequence for Signal x2(n): ' '); ); n=max(length(x1),length(x2)); %Perform Line %Perform Linear ar Convoluti Convolution on using CONV CONV command command y=cconv(x1,x2,n); %Plot Input and Convolved Signals subplot(3,1,1); stem(x1); title('Input title('Input Signal x1(n)'); x1(n)'); xlabel('n' xlabel('n'); ); ylabel('x1(n)' ylabel('x1(n)'); ); subplot(3,1,2); stem(x2); title('Input title('Input Signal x2(n)'); x2(n)'); xlabel('n' xlabel('n'); ); ylabel('x2(n)' ylabel('x2(n)'); ); subplot(3,1,3); stem(y); title('Convolved title('Convolved Signal y(n) = x1(n)*x2(n)'); x1(n)*x2(n)'); xlabel('n' xlabel('n'); ); ylabel('y(n)' ylabel('y(n)'); ); % Display the convolved Sequence in Command Window disp('Convolved disp('Convolved sequence:'); sequence:'); disp(y);
OUTPUT: Enter Input Sequence for Signal x1(n): [2 1 2 1] Enter Input Sequence for Signal x2(n): [1 2 3 4] Convolved sequence: 14 16 14 16
Dept of ECE
Page | 13
DSP Laboratory (10ECL57)
Dept of ECE
th
5 Sem
KSSEM, Bangalore
Page | 14
th
DSP Laboratory (10ECL57)
5 Sem
KSSEM, Bangalore
PROGRAM 5 AUTOCORRELATION AUTOCORRELATIO N OF A GIVEN SEQUENCE write the MATLAB MATLAB code code to perform perform Autocorre Autocorrelation lation on a given given signal signal and to to verify Aim: To write its properties. Theory: In signal signal processi processing, ng, Correlation Correlation is a measure measure of similarity similarity of of two waveforms, waveforms, as a function function of a time-lag time-lag applie applied d to one of them. them. This This is also also known known as a sliding sliding dot produc productt or sliding sliding inner-product. It is commonly used for searching searching a long-signal for a shorter, known feature. feature. Auto-correlation is a special special form of Correlation, in which a signal signal is cross-convolved with itself. itself. It is a mathe mathematic matical al tool for finding finding repeati repeating ng patterns patterns,, such as the presence presence of a periodic signal which has been buried under noise, or identifying the missing fundamental frequency in a signal implied by its harmonic frequencies. It is often used in signal processing for analyzing functions or series series of values, such as time domain domain signals. Auto-correlation of a signal ( ) is given by the formula:
( ) =
( ) ( − )
= 0,±1, 0,±1, ±2, ±2, …
In case of finite-duration signals, we can write: | |
( ) = Where = ,
=0
≥0
( )+
( )+
( − )] =
=
This can be written as:
|
= 0,
and
Assuming our signal is of the form Energy in the signal is given by: [
( ) ( − ) =
( − )
( )+
(0)+ ( )| ≤ |
< 0.
( − )+ 2
(0) + 2 (0)| =
( ) ( − )
() ≥ 0
.
That is, autocorrelation sequence of a signal attains its maximum value at zero lag. This is consistent with the notion that a signal matches perfectly with itself at zero shift. Assume a signal
( ) = {1,2,3,4}. Its autocorrelation sequence can be calculated as: ( ) = {4,11,20,30, {4,11,20,30, 20,11,4} 20,11,4}
Dept of ECE
Page | 15
th
DSP Laboratory (10ECL57)
5 Sem
KSSEM, Bangalore
Energy of the signal is given by:
=
( ) = 1 + 2 + 3 + 4 = 30 =
(0)
Hence, it is an energy signal. We see that the Total Total energy of the signal signal is equal to the amplitude of the autocorrelation signal at the origin. MATLAB CODE: % Experiment 5 : Autocorrelation of a Signal. clear all all; ; close all all; ; clc; % Accept user-defined input sequence and define time index x = input('Ent input('Enter er a finite-l finite-lengt ength h signal signal sequ sequence ence : '); n = 0:length(x)-1; % Perform Autocorrelation using xcorr function. rxx = xcorr(x,x); % Generate Time Index for Autocorrelation sequence, about origin n2 = -length( -length(x)+1 x)+1:len :length( gth(x)-1 x)-1; ; disp('Autocorrelation disp('Autocorrelation Sequence : '); '); disp(int8(rxx)); % Plot the Original Signal and Autocorrelation Sequence subplot(2,1,1); stem(n,x); title('Input title('Input Signal'); Signal'); xlabel('n' xlabel('n'); ); ylabel('x(n)' ylabel('x(n)'); ); subplot(2,1,2); stem(n2,rxx); title('Autocorrelation title('Autocorrelation Sequence'); Sequence'); xlabel('n' xlabel('n'); ); ylabel('rxx(l)' ylabel('rxx(l)'); ); grid on on; ; % Veri Verifyin fying g Autocorrela Autocorrelation tion propert properties: ies: E = sum(x.^2); % Energy of signal. mid = ceil(length(rxx)/2); % Find index of centre of sequence E0 = rxx(mid); % Detect Amplitude of Sequence midpoint fprintf('Ene fprintf('Energy rgy of Input Input Sign Signal al : %d\n %d\n' ',E); fprintf('Amplitude fprintf('Amplitude of Midpoint of Autocorrelation Sequence : %d\n',E0); %d\n' ,E0); % Verify Autocorrelation Property by comparing Energy values if int8(E0) == int8(E) %Type conversion for approximation disp('Autocorrelation disp('Autocorrelation Energy Property is verified'); else disp('Autocorrelation disp('Autocorrelation Energy Property is not verified'); end % Verify that the Signal is even. rxx_r = rxx(mid:length(rxx)); %Right Side of AC Sequence
Dept of ECE
Page | 16
th
DSP Laboratory (10ECL57)
5 Sem
rxx_l = rxx(mid:-1:1);
KSSEM, Bangalore %Left Side of AC Sequence
if rxx_r == rxx_l disp('Autocorrelation disp('Autocorrelation Sequence is Even. Hence, verified.'); else disp('Autocorrelation disp('Autocorrelation Sequence is not Even. Hence, not verified.'); verified.' ); end
OUTPUT: Enter a finite-length signal sequence : [1 2 3 4] Autocorrelation Sequence : 4 11 20 30 20 11
4
Energy of Input Signal : 30 Amplitude of Midpoint of Autocorrelation Sequence : 30 Autocorrelation Energy Property is verified Autocorrelation Sequence is Even. Hence, verified.
Dept of ECE
Page | 17
th
DSP Laboratory (10ECL57)
5 Sem
KSSEM, Bangalore
PROG PROGRA RAM M6 CROSS - CORRELATI CORRELATION ON OF A GIVEN GIVEN SEQUEN SEQUENCE CE code to perform cross-correlation on on a given signal signal and to verify Aim: To write the MATLAB code its properties. Theory:
Cross-cor Cross-correlati relation on is is a process, process, in which which a signal signal is convolve convolved d with another another signal. signal. It is commonly used for searching a long-signal for a shorter, known feature. It also has applications in pattern recognition, single particle analysis, electron tomographic averaging, crypta cryptana nalys lysis, is, and and neuro neurophy physio siology logy.. CrossCross-co corre rrelati lation on of two signa signals ls 1( ) and 2( ) is given by the formula:
( ) =
( ) ( − )
= 0, ±1, ±1, ±2, ±2, …
In case of finite-duration signals, we can write: | |
( ) = Where = ,
=0
≥0
( )+
( )+
( − )] =
=
This can be written as:
= 0,
and
Assuming our signal is of the form Energy in the signal is given by: [
( ) ( − ) =
( − )
( )+
(0)+
() ≤
( − )+ 2
(0) + 2
( 0)
< 0.
( 0) =
( ) ( − )
() ≥ 0
.
Note that the shape of the autocorrelation sequence does not change with amplitude scaling of input signals. Only the t he amplitude of the autocorrelation sequence changes accordingly. Cross-correlation satisfies the property:
( ) =
Dept of ECE
( − )
Page | 18
th
DSP Laboratory (10ECL57)
5 Sem
KSSEM, Bangalore
MATLAB CODE: % Experiment 6 : Cross-correlation of two Signals. clear all all; ; close all all; ; clc; % Accept user-defined input sequences and define time index for it x1 = input('Enter input('Enter a finite-length signal sequence X1(n): '); n1 = 0:length(x1)-1; 0:length(x1)-1; x2 = input('Enter input('Enter a finite-length signal sequence X2(n): '); n2 = 0:length(x2)-1; 0:length(x2)-1; x= max(x1,x2); % Perfor Perform m Cross Cross - Cor Corre relat latio ion n using using xcorr xcorr funct functio ion. n. rxy = xcorr(x1,x2); % rxy(l) ryx = xcorr(x2,x1); % ryx(l) % Generat Generate e Time Time Index Index for Cros Cross s - Corr Correlat elation ion seque sequence, nce, abou about t origin origin n3 = -length( -length(x)+1 x)+1:len :length( gth(x)-1 x)-1; ; disp('Cr disp('Cros oss s - Cor Correl relati ation on Se Seque quenc nce e rxy( rxy(l): l): '); disp(int8(rxy)); disp('Cr disp('Cros oss s - Cor Correl relati ation on Se Seque quenc nce e ryx( ryx(l): l): '); disp(int8(ryx)); % Plot Plot the the Origi Original nal Sign Signal al and and Cross Cross - Corr Correlat elation ion Sequ Sequence ence subplot(3,1,1); stem(n1,x1); title('Input title('Input Signal 1'); 1'); xlabel('n' xlabel('n'); ); ylabel('x1(n)' ylabel('x1(n)'); ); subplot(3,1,2); stem(n2,x2); title('Input title('Input Signal 2'); 2'); xlabel('n' xlabel('n'); ); ylabel('x2(n)' ylabel('x2(n)'); ); subplot(3,1,3); stem(n3,rxy); title(' title('Cr Cro oss - Cor orre rel lat atio ion n Se Sequ que enc nce' e'); xlabel('n' xlabel('n'); ); ylabel('rxy(l)' ylabel('rxy(l)'); ); grid on on; ;
% Verifying
Cross-correlation Cross-correla tion properties:
E1 = sum(x1.^2); E2 = sum(x2.^2); mid = ceil(length(rxy)/2); E0 = abs(max(rxy)); abs(max(rxy)); fprintf('Energy fprintf('Energy of Input Signal fprintf('Energy fprintf('Energy of Input Signal fprintf('M fprintf('Max ax Am Ampl plit itud ude e of Cr Cros oss s
% Energy of signal 1. % Energy of signal 2. % Find Find index index of of cent centre re o of f sequenc sequence e % Detect Max Amplitude of Sequence X1 : %d\n',E1); %d\n',E1); X2 : %d\n',E2); %d\n',E2); - Co Corr rrel elat atio ion n Se Sequ quen ence ce : %d %d\n \n' ',E0);
% Verify Verify Cross Cross - Corr Correlat elation ion Prope Property rty by by compari comparing ng Energy Energy valu values es % Max amplitude of Sequence should be less than sqrt(E1*E2). if int8(E0) <= int8(sqrt(E1*E2)) %Type conversion to 8-bit int disp('Cr disp('Cross oss - Cor Corre relat latio ion n Energy Energy Prope Property rty is veri verifie fied' d'); else disp('Cr disp('Cross oss - Cor Corre relat latio ion n Energy Energy Prope Property rty is not not verifie verified' d'); end
Dept of ECE
Page | 19
DSP Laboratory (10ECL57)
th
5 Sem
KSSEM, Bangalore
% Verify Signal property : rxy(l)=ryx(-l). if rxy == fliplr(ryx) disp('Sin disp('Since ce rxy( rxy(l) l) = ryx(ryx(-l), l), Cros Cross s - Corr Correlat elation ion prop property erty is verified.'); verified.' ); else disp('C disp('Cro ross ss - Co Corr rrel elat atio ion n pr prop oper erty ty is is not not veri verifi fied ed.' .'); end
OUTPUT: Enter a finite-length signal sequence X1(n): [4 3 2 1] Enter a finite-length signal sequence X2(n): [1 2 3 4] Cross Cross - Correlatio Correlation n Sequen Sequence ce rxy(l): rxy(l): 16 24 25 20 10 4 1 Cross Cross - Correlatio Correlation n Sequen Sequence ce ryx(l): 1 4 10 20 25 24 16 Energy of Input Signal X1 : 30 Energy of Input Signal X2 : 30 Max Max Amplitu Amplitude de of Cross Cross - Correl Correlati ation on Sequ Sequenc encee : 25 Cross Cross - Correlatio Correlation n Energy Energy Property Property is verifie verified d Since Since rxy(l) = ryx(-l), ryx(-l), Cross - Correlatio Correlation n property property is verified. verified.
Dept of ECE
Page | 20
DSP Laboratory (10ECL57)
th
5 Sem
KSSEM, Bangalore
PROGRAM 7 SOLVING A GIVEN DIFFERENCE EQUATION Aim: To write the MATLAB code to solve a given difference equation, given the co-efficients and initial values. Let us consider the difference equation as y (n) – 3/2 y (n-1) + ½ y (n-2) = x (n). Given n x(n) = (1/4) *u(n). Assume initial conditions as y(-1) = 4, y(-2) = 10. Theory: n
Consider x(n) = (1/4) *u(n). Let n take values from 0 to 5, n=0:5 n=0, x(0)=1 n=1, x(1)=0.25 n=2, x(2)=0.0625 n=3, x(3)=0.0156 n=4, x(4)=0.0039 n=5, x(5)=0.0010 For n=0; y(0) - 3/2 y(0-1) y(0-1) + 1/2 y(0-2) y(0-2) = x(0) Substituting the initial conditions and the value of x(0) in the above equation we get, y(0) (0) = 1 + 6 - 5 = 2 Similarly, For n=1; n=1; y(1) y(1) = 0.25 + 3 - 2 = 1.2500 1.2500 For n=2; y(2) y(2) = 0.0625 0.0625 + 1.875 -1 = 0.9375 0.9375 For n=3; n=3; y(3) = 0.0156 0.0156 + 1.4062 1.40625 5 - 0.625 0.625 = 0.7969 0.7969 For n=4; n=4; y(4) = 0.0039 0.0039 + 1.1953 1.19535 5 - 0.46625 0.46625 = 0.7365 0.7365 For n=5; n=5; y(5) = 0.0010 0.0010 + 1.0957 1.09575 5 - 0.3982 0.3982 = 0.6982 0.6982 MATLAB CODE: %Experiment 7 : Difference Equation Solving clear all all; ; close all all; ; clc; %Accept Difference Equation Coefficients from Input b = input('Enter input('Enter the coefficients of input x(n) : '); a = input('Ent input('Enter er the c coeff oefficie icients nts of outpu output t y(n) : '); y = input('Enter input('Enter the initial conditions y(-1), y(-2),... : '); %Calculate Initial Conditions using filtic z = filtic(b,a,y); %Ignore if No Initial Conditions %Enter Input sequence samples. x = [1 1/4 1/16 1/64 1/256 1/1024]; n = 0:length(x)-1; %Time Base for plotting
Dept of ECE
Page | 21
DSP Laboratory (10ECL57)
th
5 Sem
KSSEM, Bangalore
%Calculate output using initial conditions Yout = filter(b,a,x,z); %Use filter(b,a,x) if no IC %Display output sequence disp('Difference disp('Difference Equation Solution : y(n) : '); '); disp(Yout); %Plot Input and Output subplot(2,1,1); stem(n,x); title('Input title('Input Sequence x(n)'); x(n)'); xlabel('n' xlabel('n'); ); ylabel('x(n)' ylabel('x(n)'); ); subplot(2,1,2); stem(n,Yout); grid on on; ; title('Output title('Output Sequence y(n)'); y(n)'); xlabel('n' xlabel('n'); ); ylabel('y(n)' ylabel('y(n)'); );
OUTPUT: Enter the coefficients of input x(n) : 1 Enter Enter the coeffic coefficient ients s of output output y(n) : [1 [1 -3/2 1/2] 1/2] Enter the initial conditions y(-1), y(-2),... : [4 10] Difference Equation Solution : y(n) : 2.0000 1.2500 0.9375 0.7969 0.6982
Dept of ECE
0.7305
Page | 22
th
DSP Laboratory (10ECL57)
5 Sem
KSSEM, Bangalore
PROGRAM 8 COMPUTATI COMPUTATION ON OF N- POINT POINT DFT point DFT of a given given sequence sequence and to plot magnitude magnitude and phase Aim: Computation of N point spectrum. Theory: DFT stands for Discrete Fourier Transform. It is used to find the amplitude and phase spectrum of a discrete time sequence. sequence. The N-point DFT of a finite-length sequence x(n) of length N, is given by X(k) as N-DFT
x(n)
X(k)
Basic equation to find the DFT of a sequence is given below.
( ) =
( )
For example: Find the DFT of the sequence x(n) = {0,1,2,3} In this case N=4. For k=0,
( )
(0) =
= ( 0) + ( 1) + ( 2) + ( 3) = 0 +1 +1 + 2 + 3 = 6
For k=1,
( )
(1) =
= ( 1)
.
+ ( 1)
+ ( 2)
+ ( 3)
= - j- 2 + 3j = -2 + j2 Similarly, For k=2,
(2) =
( )
= ( 1)
.
+ ( 1)
+ ( 2)
+ ( 3)
= 0-1+ 0-1+2-3 2-3 = -2 For k=3, X(3) X(3) = -2 – j2 Hen Hence, X(k X(k) = {6, {6, -2+ -2+j2, j2, -2, -2, -2-j -2-j2 2}
Dept of ECE
Page | 23
th
DSP Laboratory (10ECL57)
5 Sem
KSSEM, Bangalore
MATLAB CODE: %Experiment 8 : N-Point DFT clear all all; ; close all all; ; clc; %Accept Input sequence from user xn = input ('Enter ('Enter the sequence x(n) : '); '); xn=xn'; N = length(xn); length(xn); Xk = zeros(N, 1); %Initialize zero matrix for DFT sequence %Calculate DFT using formula n = 0:N-1; for k = 0:N-1 Xk(k+1) = exp(-j*2*pi*k*n/N)*xn; end %Display DFT Sequence disp('DSP disp('DSP Sequence : X(k) :'); :'); disp(int8(Xk)); %Plot Signals n = 0:N-1;
%Time base
% Input Sequence subplot (2,2,[1:2]); stem(n, xn); title('Input title('Input Sequence x(n)'); x(n)'); xlabel('n' xlabel('n');ylabel( );ylabel('x(n)' 'x(n)'); ); % Output Magnitude Plot subplot (2,2,3); stem(n, abs(Xk)); grid on on; ; title('Magnitude title('Magnitude Plot of DFT : |X(k)|'); |X(k)|'); xlabel('n' xlabel('n');ylabel( );ylabel('|X(k)|' '|X(k)|'); ); % Output Phase Plot subplot(2,2,4); stem(n, angle(Xk)'); grid on on; ; title('Pha title('Phase se Plot Plot of of DFT DFT : angl angle(X( e(X(k))' k))'); xlabel('n' xlabel('n');ylabel( );ylabel('Angle' 'Angle'); );
OUTPUT: Enter the sequence x(n) : [0 1 2 3] DSP Sequence : X(k) : 6.0000 + 0.0000i -2.0000 + 2.0000i -2.0 -2.000 000 0 - 0.00 0.0000 00ii -2.0 -2.000 000 0 - 2.00 2.0000 00ii
Dept of ECE
Page | 24
DSP Laboratory (10ECL57)
Dept of ECE
th
5 Sem
KSSEM, Bangalore
Page | 25
th
DSP Laboratory (10ECL57)
5 Sem
KSSEM, Bangalore
PROGRAM 9 LINEAR CONVOLUTION USING DFT DFT AND IDFT IDFT Aim: To calculate the Linear Convolution of two sequences using DFT and IDFT Theory: An interesting interesting property of the Discrete Fourier Fourier Transforms, Transforms, is the effect effect it has has on convolution. Convolution of two signals in the time domain translates to a multiplication of their Fourier transforms in the frequency domain. In this procedure, we find the discrete Fourier transforms of the individual signals, multiply them, and apply an Inverse Fourier Transform upon the product, to get the convolved signal in the time domain.
If x(n) and h(n) are the two sequences of length ‘l’ and ‘m’ respectively. respectivel y. then X(k) and H(k) their DFT’s of length N=L+M-1. Y(k)=x(k)h(k) Therefore the linear convolution of two sequence is the N point IDFT of Y(k). Ex: Find the linear convolution of x(n)={1,2} and h(n)={1,2,3} using DFT and IDFT method. Soln: Given,
x(n)={1,2} h(n)={1,2,3} here L = 2, M = 3
N=L+M-1
Therefore, N=4
x(n)={1,2,0,0} and h(n)={1,2,3,0} Finding X(k) using DIT FFT algorithm: X(k) = {3 {3 , 1-2j , -1 , 1+2j 1+2j } Finding H(k) using DIT FFT algorithm H(k) = {6 , -2-2j , 2 , -2+2j } Product Y(k) is calculated as : Y(k) = X(k)H(k)
( ) = {18 ,−6 ,−6 + 2j ,−2 ,−2 ,6 − 2j } Finding y(n) using DIT FFT algorithm:
y(n) (n) = {1 ,4,7 ,6 } Dept of ECE
Page | 26
DSP Laboratory (10ECL57)
th
5 Sem
KSSEM, Bangalore
MATLAB CODE % Expe Experime riment nt 9 : Linear Linear Convo Convoluti lution on using using Four Fourier ier Tran Transfor sforms ms clear all all; ; close all all; ; clc; %Accept input sequences x1 = input('Ent input('Enter er Input Input Sequ Sequence ence f for or Signal Signal x1(n x1(n): ): '); n1 = length(x1); x2 = input('Enter input('Enter Input Sequence for Signal x2(n): ' '); ); n2=length(x2); N = n1+n2-1; T = 1:N;
% Length of convolved sequence
%Calculate N-point DFT and IDFT. y1=fft(x1,N); % N-point DFT of x1 y2=fft(x2,N); % N-point DFT of x2 y3=y1.*y2; % Multiplication in time domain y=ifft(y3,N); % N-point IDFT of y to recover result % Plot Input and Output Sequences: subplot(3,1,1); stem(x1); title('Input title('Input Signal x1(n)'); x1(n)'); xlabel('n' xlabel('n'); ); ylabel('x1(n)' ylabel('x1(n)'); ); subplot(3,1,2); stem(x2); title('Input title('Input Signal x2(n)'); x2(n)'); xlabel('n' xlabel('n'); ); ylabel('x2(n)' ylabel('x2(n)'); ); subplot(3,1,3); stem(y); title('Convolved title('Convolved Signal y(n) = x1(n)*x2(n)'); x1(n)*x2(n)'); xlabel('n' xlabel('n'); ); ylabel('y(n)' ylabel('y(n)'); ); % Display the convolved Sequence in Command Window disp('Convolved disp('Convolved sequence:'); sequence:'); disp(y);
OUTPUT Enter Input Sequence for Signal x1(n): [1 2] Enter Input Sequence for Signal x2(n): [1 2 3] Convolved sequence: 1 4 7 6
Dept of ECE
Page | 27
DSP Laboratory (10ECL57)
Dept of ECE
th
5 Sem
KSSEM, Bangalore
Page | 28
DSP Laboratory (10ECL57)
th
5 Sem
KSSEM, Bangalore
PROGRAM 10 CIRCULAR CONVOLUTION USING DFT AND IDFT Aim: To calculate the Circular Convolution of two sequences using DFT and IDFT Theory: Convolution in time domain domain corresponds corresponds to multiplication in frequency frequency domain. To make use of this property, we can calculate the DTFT of each of the aperiodic signals, multiply these in the frequency domain, and find the IDFT of the product, to get the periodic convolved signal in time domain. Example: Find the circular convolution convolution of x(n)={1,2,3,4} and h(n)={4,3,2} h(n)={4,3,2} using DFT and and IDFT method. Solution: Solution: Given Given two signals, signals, x(n)={1,2 x(n)={1,2,3,4} ,3,4} and h(n)={4,3, h(n)={4,3,2} 2} Finding X(k) using DIT FFT algorithm X(k) X(k) = {10 {10 , -2-2 -2-2jj , -2 , -2-2j -2-2j } Finding H(k) using DIT FFT algorithm H(k) = { 9 , 2-3j , 3 , 2+3j } Y(k) = X(k)H(k) Y(k) ={ 90 90 , 2+10j 2+10j , -6 , 2-10j } Finding y(n) using DIT FFT algorithm
y(n) (n) = {2 {22 2 ,19 ,19 ,20 ,20 ,29 ,29 } MATLAB CODE % Experim Experiment ent 10 10 - Circ Circular ular Conv Convolut olution ion using using Four Fourier ier Trans Transform forms s clear all all; ; close all all; ; clc; %Accept input sequences x1 = inpu input t('Enter Input Sequence for Signal x1(n): '); '); n1 = length(x1); x2 = input('Enter input('Enter Input Sequence for Signal x2(n): '); '); n2=length(x2); N=max(n1,n2); T = 1:N;
% Length of convolved sequence
x1=[x1 zeros(1,N-n1)]; x2=[x2 zeros(1,N-n2)];
Dept of ECE
Page | 29
DSP Laboratory (10ECL57)
th
5 Sem
KSSEM, Bangalore
%Calculate N-point DFT and IDFT. y1=fft(x1,N); % N-point DFT of x1 y2=fft(x2,N); % N-point DFT of x2 y3=y1.*y2; % Multiplication in time domain y=ifft(y3,N); % N-point N-point IDFT of y to reco recover ver resul result t % Plot Input and Output Sequences: subplot(3,1,1); stem(x1); title('Input title('Input Signal x1(n)'); x1(n)'); xlabel('n' xlabel('n'); ); ylabel('x1(n)' ylabel('x1(n)'); ); subplot(3,1,2); stem(x2); title('Input title('Input Signal x2(n)'); x2(n)'); xlabel('n' xlabel('n'); ); ylabel('x2(n)' ylabel('x2(n)'); ); subplot(3,1,3); stem(y); title('Convolved title('Convolved Signal y(n) = x1(n)*x2(n)'); x1(n)*x2(n)'); xlabel('n' xlabel('n'); ); ylabel('y(n)' ylabel('y(n)'); ); grid on on; ; % Display the convolved Sequence in Command Window disp('Convolved disp('Convolved sequence:'); sequence:'); disp(y);
OUTPUT Enter Input Sequence for Signal x1(n): [1 2 3 4] Enter Enter Input Input Sequenc Sequence e for Signal Signal x2(n) x2(n): : [4 3 2] Convolved sequence: 22 19 20 29
Dept of ECE
Page | 30
DSP Laboratory (10ECL57)
Dept of ECE
th
5 Sem
KSSEM, Bangalore
Page | 31
th
DSP Laboratory (10ECL57)
5 Sem
KSSEM, Bangalore
PROGRAM 11 DESIGN AND IMPLEMENTATION OF FIR FILTER Aim: To design and implement a FIR Filter for the given specifications. Theory: A linear-phase is required throughout the passband of the filter to preserve the shape of the given signal in the passband. A causal IIR filter cannot give linear-phase characteristics and only special types of FIR filters that exhibit center symmetry in its impulse response give the linear-space. An FIR filter with wit h impulse response h(n) can be obtained as follows: h(n) = hd(n) 0≤n≤N-1
= 0 otherwise ……………….(a) The impulse response h d(n) is truncated at n = 0, since we are interested in causal FIR Filter. It is possible to write above equation alternatively as h(n) = hd(n)w(n) ……………….(b) where w(n) is said to be a rectangular window defined by w(n) = 1 0≤n≤N -1 = 0 otherwise Taking DTFT on both the sides of equation(b), we get
H(ω) = Hd(ω)*W(ω)
Hamming window: window: The impulse response of an N-term Hamming window is defined as follows:
( ) =
0. 0.54 – 0.46 0
(2п / ( − 1) 1))
0≤ ℎ
≤
−1
Problem: Using MATLAB design an IIR filter fil ter to meet the following specifications choosing Hamming window: • • • •
Dept of ECE
Window length, N = 27 Stop band attenuation = 50dB Cut-off frequency = 100 Hz Sampling frequency = 1000 Hz
Page | 32
DSP Laboratory (10ECL57)
th
5 Sem
KSSEM, Bangalore
MATLAB CODE % Experiment 11 : Designing a FIR Filter (Hamming Window) close all all; ; clear all all; ; clc; % Accept Filter Parameters from User N = input('Enter input('Enter the window length N : '); '); fc = inpu input t('Enter the cut-off frequency fc (Hz) : '); '); Fs = input('Enter input('Enter the sampling frequency Fs (Hz) : '); '); Wc = 2*fc/ Fs; Wh = hamming(N); window
%Nyquist Frequency %Create a N-point symmetric Hamming
% Generate Generate a FIR filter filter based on on Hamming Hamming Window Window created created b = fir1(N-1, Wc ,Wh); % Calculate Frequency Response of Filter designed. % h - Fre Freque quency ncy Res Respon ponse se val values ues w = Fr Frequ equen encie cies s [h,w] = freqz(b,1,256); mag = 20*log10(abs(h)); %Magnitude of Response % Di Disp spla lay y Val Value ues s disp('Hamming disp('Hamming Window Co-efficients : '); '); disp(Wh); disp('Unit disp('Unit Sample Response of FIR Filter h(n) : '); disp(b); % Plot Frequency response of Butterworth Filter. freqz(b); title('Hamming title('Hamming Filter Frequency Response'); Response');
OUTPUT Enter the window length N : 27 Enter the cut-off frequency fc (Hz) : 100 Enter the sampling frequency Fs (Hz) : 1000 Hamming Window Co-efficients : 0.0800 0.0934 0.1327 0.1957 0.2787 0.3769 0.4846 0.5954 0.7031 0.8013
Dept of ECE
Page | 33
DSP Laboratory (10ECL57)
th
5 Sem
KSSEM, Bangalore
0.8843 0.9473 0.9866 1.0000 0.9866 0.9473 0.8843 0.8013 0.7031 0.5954 0.4846 0.3769 0.2787 0.1957 0.1327 0.0934 0.0800 Unit Sample Response of FIR Filter h(n) : Columns 1 through 7 0.0019 0.0 0.0023 0.0 0.0022 -0.0000 -0.0058 -0.0142 -0.0209 Columns 8 through 14 -0.0185 0.0000 0.0374
0.0890
0.1429
Columns 15 through 21 0.18 0.1840 40 0.1 0.142 429 9 0.08 0.0890 90
0.0 0.037 374 4
0.00 0.0000 00 -0.0 -0.018 185 5 -0.0 -0.020 209 9
Columns 22 through 27 -0.0 -0.014 142 2 -0.0 -0.005 058 8 -0.0 -0.000 000 0
0.00 0.0022 22
Dept of ECE
0.00 0.0023 23
0.1840
0.1994
0.00 0.0019 19
Page | 34
DSP Laboratory (10ECL57)
Dept of ECE
th
5 Sem
KSSEM, Bangalore
Page | 35
DSP Laboratory (10ECL57)
th
5 Sem
KSSEM, Bangalore
PROGRAM 12 DESIGN AND IMPLEMENTATION OF IIR FILTER Aim: To design and implement a IIR Filter for the given specifications. Theory: A desired frequency response is approximated by a transfer function expressed as a ratio of polynomials. This type of transfer function yields an impulse response of infinite duration. Therefore, the analog filters are commonly referred to as infinite impulse response (IIR) filters. The main classes of analog filters are 1.Butterworth Filter. 2.Chebyshev 2.Chebyshev Filter. These filters differ in the nature of their magnitude responses as well as in their design and implementation. BUTTERWORTH FILTERS: FILTERS: Butterworth filters have very smooth passband, which we pay for with a relatively rel atively wide transition region. A Butterworth filter is characterized by its magnitude frequency response, where N is the order of the filter and Ωc √2 times the dc gain (Ω=0) magnitude is 1/ √2
Dept of ECE
is defined as the cutoff frequency where the filter
Page | 36
th
DSP Laboratory (10ECL57)
5 Sem
KSSEM, Bangalore
Butterworth filter tables N=1; N=1; (s + 1) 1) 2
N=2; (s +0.5 s +1) 2
N=3; (s +s+1)(s+1) 2
2
N=4; (s +0.76536s+1)(s +1.864776s+2) 2 2 N=5; N=5; (s+1)( (s+1)( s +0.6180s+1)( +0.6180s+1)( s +1.6180s+1) CHEBYSHEV FILTERS: Chebyshev Chebyshev filters are equiripple in either the passband or stopband. Hence the magnitude response oscillates between the permitted minimum and maximum values in the band a number of times depending upon the order of filters. There are two types of chebyshev filters. The chebyshev I filter is equiripple in passband and monotonic in the stopband, whereas Chebyshev Chebyshev II is just the opposite. The Chebyshev Chebyshev low-pass filter has a magnitude response response given by
Dept of ECE
Page | 37
DSP Laboratory (10ECL57)
th
5 Sem
KSSEM, Bangalore
PROBLEM 1: BUTTERWORTH FILTER: Using MATLAB design an IIR filter filt er with passband edge frequency 1500Hz and stop band edge at 2000Hz for a sampling frequency of 8000Hz, variation of gain within pass band 1db and stopband attenuation of 15 db. Use Butterworth prototype design and Bilinear Transformation. MATLAB CODE : % Experiment 12 : Design of IIR Filter (Butterworth Filter) clear all all; ; close all all; ; clc; % Accept Input Parameters from user Rp = input('Enter input('Enter Passband Attenuation in dB Rs = input('Enter input('Enter Stopband Attenuation in dB fp = input('Enter input('Enter Passband Frequency in Hz : fs = input('Ent input('Enter er Stopba Stopband nd Freq Frequenc uency y in Hz Hz : Fs = input('Enter input('Enter Sampling Frequency in Hz :
: '); '); : '); '); '); '); '); '); ');
% Calculate Sampled Frequency values Wp=2* fp / Fs ; Ws=2* fs / Fs ; % Calculate Butterworth filter order and cutoff frequency: % N = Minimum order of Filter Wn = Cutoff Frequencies [N,Wn] = buttord(Wp,Ws,Rp,Rs); % Butterworth Filter Design (z = zeros [z,p] = butter(N,Wn);
p = poles)
% Display Filter parameters : disp('Order disp('Order of Butterworth Filter : '); '); disp(N); disp('Butterworth disp('Butterworth Window Cutoff Frequency : '); disp(Wn); % Plot Frequency Response of Filter freqz(z,p); title('Butterworth title('Butterworth frequency response'); response');
OUTPUT: Enter Passband Attenuation Attenuation in dB : 1 Enter Stopband Attenuation in dB : 15 Enter Passband Frequency in Hz : 1500 Enter Stopband Stopband Frequency in Hz : 2000 2000 Enter Sampling Frequency in Hz : 8000 Order of Butterworth Filter : 6 Dept of ECE
Page | 38
DSP Laboratory (10ECL57)
th
5 Sem
KSSEM, Bangalore
Butterworth Window Co-efficient : 0.4104
PROBLEM PROBLEM 2: CHEBYSHEV CHEBYSHEV FILTER: FILTER: Using MATLAB design design an IIR filter with passband edge frequency 1500Hz 1500Hz and stop band edge at 2000Hz for a sampling frequency of 8000Hz, variation of gain within pass band 1 db and stopband attenuation of 15 db. Use Chebyshev prototype design and Bilinear Transformation.
DESIGN: W1 = (2*pi* F1 )/ Fs
= 2*pi*100)/4000 2*pi*100)/4000 = 0.05Π rad W2 = (2*pi* F2 )/ Fs = 2*pi*500)/4000 =0.25Π rad Prewarp: T=1sec
Ω1 = 2/T tan (w1/2) (w1/2) = 0.157 rad/sec Ω2 = 2/T tan (w2/2) = 0.828 rad/sec
Dept of ECE
Page | 39
th
DSP Laboratory (10ECL57)
5 Sem
KSSEM, Bangalore
Order:
έ = 10 − 1 έ = 0.765 -As/20
20/20
A= 10 , A = 10 , A=10 g= (A2 (A2 - 1) / έ , g = 13.01
Ωr= Ω2 / Ω1 Ωr=0.828/0.157 = 5.27 rad \sec n = log log10
+
(
− )
/ log10{Ωr + √(
– )}
n= 1.388 Therefore n= 2. Cut-off Frequency:
Ωc = Ωp = Ω1 = 0.157 rad \sec Normalized Transfer Function: 2
2
H(s)=[bo / 1+ έ ] / [ s +b1s+b0] 2 = 0.505/[ s +0.8s+0.036] Denormalized Transfer Function: H(s)= Hn(s) |s-s/Ωc H(s)= Hn(s) |s-s/0.157 2 H(s) = 0.0125 / [s +0.125s+0.057] Apply BLT: H(Z) = H(s)|s=(2/T)[(1-z-1)/(1+z-1)] -1
-2
H(Z) = 0.0125+0.025Z + 0.01 0.0125 25 Z -1 4.2769-7.96Z + 3.76Z-2 -1
H(Z) = 0.0029+0.0052Z + 0.0029 Z -1 -2 1-1.86Z + 0.88Z
Dept of ECE
-2
Page | 40
DSP Laboratory (10ECL57)
th
5 Sem
KSSEM, Bangalore
MATLAB CODE: % Experiment 12B : Design of IIR Filter (Chebyshev Filter) clear all all; ; close all all; ; clc; % Accept Input Parameters from user Rp = inpu input t('Enter Passband Attenuation in dB Rs = input('Enter input('Enter Stopband Attenuation in dB fp = input('Enter input('Enter Passband Frequency in Hz : fs = input('Enter input('Enter Stopband Frequency in Hz : Fs = input('Enter input('Enter Sampling Frequency in Hz :
: '); '); : '); '); '); '); '); '); '); ');
% Calc Calculat ulate e Sampled Sampled Freque Frequency ncy values values Wp=2* fp / Fs ; Ws=2* fs / Fs ; % Calculate Chebyshev filter order and cutoff Frequency: % N = Minimum order of Filter Wn = Cutoff Frequencies [N,Wn]=cheb1ord(Wp,Ws,Rp,Rs); % Ch Chebyshev Fi Filter De Design (z (z = zeros [z,p]=cheby1(N,Rp,Wn);
p = p ol es )
% Display Filter parameters : disp('Order disp('Order of Chebyshev Filter : '); '); disp(N); disp('Chebyshev disp('Chebyshev Window Cutoff Frequency : '); '); disp(Wn); % Plot Frequency Response of Filter : freqz(z,p); title('Che title('Chebysh byshev ev Freq Frequenc uency y respo response' nse');
OUTPUT: Enter Passband Attenuation Attenuation in dB : 1 Enter Stopband Attenuation in dB : 15 Enter Passband Frequency in Hz : 1500 Enter Stopband Frequency in Hz : 2000 Enter Sampling Frequency in Hz : 8000 Order of Chebyshev Filter : 4 Chebyshev Chebyshev Window Cutoff Frequency : 0.3750
Dept of ECE
Page | 41
DSP Laboratory (10ECL57)
Dept of ECE
th
5 Sem
KSSEM, Bangalore
Page | 42
DSP Laboratory (10ECL57)
th
5 Sem
KSSEM, Bangalore
PART B: Exercises using the DSP Kit
Dept of ECE
Page | 43
DSP Laboratory (10ECL57)
V Semester
KSSEM, Bangalore
TMS320C6748 DSP BOARD
Package content:
The C6748 DSP Experime Experimenter nter Kit packaged packaged for Universitie Universitiess (TMDSEXP (TMDSEXPL138L138-UNV) UNV) provid provides es everyt everything hing academics academics need need to get starte started d with teaching teaching and projec projects ts using using the TI TMS320C6000 DSP. The TI TMS320C6000 family of DSPs is popular for use in Real-time DSP coursesThe Experimenter Kit for Universities is a low-cost, flexible development platform for the OMAP-L138 which is a low-power dual core processor based on C6748 DSP plus an ARM926EJ-S 32-bit RISC MPU.
The C6748 DSP kit has a TMS320C6748 DSP onboard that allows full -speed verification of code with Code Composer Studio. The C6748 DSP kit provides: A USB Interface 128MB DDRAM and ROM An analog interface circuit for Data conversion (AIC) An I/O port Embedded JTAG emulation support
Connectors on the C6748 DSP kit provide DSP external memory interface (EMIF) and peripheral signals that enable its functionality to be expanded with custom or third party daughter boards.
Dept of ECE
Page | 44
DSP Laboratory (10ECL57)
V Semester
KSSEM, Bangalore
The C6748 DSP kit includes a stereo codec. This analog interface circuit (AIC) has t he following characteristics: High-Performance Stereo Codec • Interfaces directly to digital or analog microphones • Supports 8-96 ksps sampling rates • High High SNR (100-102dB DAC, 92dB ADC) • Integrated PLL supporting a wide range of audio clocks • Low-power headphone, speaker and playback modes for portable s ystems • Programmable digital audio effectsinclude 3D sound, bass, t reble, EQ and de-emphasis Softwa Software re Control Control Via TI McASPMcASP-Com Compat patible ible Multip Multiprot rotoco ocoll Serial Serial Port.G Port.Glue lueles lesss Interface to TI McASPs. Audio-Da Audio-Data ta Input/Outp Input/Output ut Via TI McASP-Com McASP-Compatible patible Programma Programmable ble Audio Audio Interface Interface 16/20/24/32-Bit Word Lengths. The C6748DSP kit has the following f ollowing features: The 6748 DSP KIT kit is a low-cost standalone development platform that enables customers to evaluate and develop applications applications for the TI C67XX DSP family. family. The DSP KIT also also serve servess as a hardware referen reference ce design for the TMS320C6 TMS320C6748 748 DSP. DSP. Schematics, Schematics, logic logic equations equations and application notes are available to ease hardware development and reduce time to market. An on-boa on-board rd AIC31 AIC3106 06 codec codec allows allows the DSP DSP to transm transmit it and receiv receivee analog analog signals signals.. McASP McASP is used for the codec control control interface and for data. data. Analog Analog audio I/O is done through through two 3.5mm 3.5mm audio audio jacks that that correspo correspond nd to line input, input, and line. The analog analog output output is driven driven to the line out .McASP1 can be be re-routed to the expansion connectors in software. The DSP KIT includes 2 LEDs and 8 DIP switches as a simple way to provide the user with interactive feedback. An includ included ed 5V externa externall power supply supply is is used to power power the board. board. On-boa On-board rd voltag voltagee regula regulator torss provid providee the 1.26V 1.26V DSP DSP core core voltag voltage, e, 3.3V 3.3V digital digital and 3.3V analog voltag voltages. es. A voltag voltagee superv superviso isorr monitor monitorss the intern internally ally generate generated d voltag voltage, e, and will will hold the board board in reset reset until the supplies are within operating specifications and the reset button is released. Code Composer communicates with the DSP KIT through an embedded JTAG emulator with a USB USB host host interf interfac ace. e. The The DSP DSP KIT KIT can also also be used used with with an exte extern rnal al emulat emulator or throug through h the external JTAG connector. TMS320C6748 DSP Features
Highest-Performance Floating-Point Digital Digital Signal Processor (DSP): Eight 32-Bit Instructions/Cycle
Dept of ECE
Page | 45
DSP Laboratory (10ECL57)
V Semester
KSSEM, Bangalore
32/64-Bit Data Word 375/456-MHz C674x Fixed/Floating-Point 3648/2746 C674x MIPS/MFLOPS MIPS/MFLOPS Up to 3648/2746 Peripheral Set, Optimized for Audio Rich Peripheral Highly Optimized C/C++ Compiler Extended Temperature Devices Available DSP Core Advanced Very Long Instruction Word (VLIW) TMS320C67x™ TMS320C67x™ DSP Core Independent Functional Functional Units: Eight Independent Two ALUs (Fixed-Point) Four ALUs (Floating- and Fixed-Point) Fixed-Point) Two Multipliers (Floating- and Fixed-Point) Architecture With 64 32-Bit General-Purpose General-Purpose Registers Registers Load-Store Architecture Packing Reduces Reduces Code Size Instruction Packing Instructions Conditional All Instructions Instruction Set Features Instructions for IEEE 754 Native Instructions Single- and Double-Precision 16-, 32-Bit 32-Bit Data) Byte-Addressable (8-, 16-, 8-Bit Overflow Protection Bit-Field Extract, Set, Clear; Bit-Counting; Bit-Counting; Normalization Saturation; Bit-Field 67x cache memory. L1P Program Cache (Direct-Mapped) (Direct-Mapped) 32K-Byte L1P L1D Data Data Cache (2-Way) 32K-Byte L1D 256K-Byte 256K-Byte L2 unified Memory RAM\Cache. Real-Time Clock With 32 KHz Oscillator and Separate Separate Power Rail. Rail. Three 64-Bit 64-Bit General-Purpose General-Purpose Timers Digital Audio Interface Interface Transmitter (DIT) Supports: Integrated Digital S/PDIF, IEC60958-1, AES-3, CP-430 Formats Up to 16 transmit pins Enhanced Channel Status/User Data Extensive Error Checking and Recovery 2 Two Inter-Integrated Inter-Integrated Circuit Bus (I C Bus™) Bus™) . 3 64-Bit General-Purpose General-Purpose Timers (each configurable as 2 32-bit timers) Flexible Phase-Locked-Loop Phase-Locked-Loop (PLL) Based Based Clock Generator Generator Module
IEEE-1149.1 IEEE-1149.1 (JTAG
3.3-V 3.3-V I/Os, I/Os, 1.2 -V Internal Internal (GDP & PYP) 3.3-V I/Os, I/Os, 1.4-V Internal Internal (GDP)(300 MHz only) LCD Controller Two Serial Peripheral Peripheral Interfaces (SPI) (SPI) Each With Multiple Chip-Selects Chip-Selects Two Multimedia Card (MMC)/Secure (MMC)/Secure Digital (SD) Card Card Interface with Secure Data Data I/O (SDIO) Interfaces One Multichannel Audio Serial Serial Port. Port. Two Multichannel Multichannel Buffered Serial Ports Ports
Dept of ECE
) BoundaryBoundary-ScanScan-Comp Compatible atible
Page | 46
DSP Laboratory (10ECL57)
V Semester
KSSEM, Bangalore
General Procedure to work on Code Composer Studio V4 for non-real time projects
1. Launch ccs Launc Launch h the the CCS v4 ico icon n from from the the Deskt Desktop op or or goto goto All Programs ->Texas Instruments >CCSv4
2. Choose Choose the locat location ion for for the workspac workspace, e, where where your your project project will will be saved saved..
3. Clic Click k the the CCS CCS icon icon from from the the welc welcom omee pag page to go the the work workbe benc nch, h, it it is ma mark rkeed in the the belo below w picture.
Dept of ECE
Page | 47
4. Configure ccs to your target A. From From the Target Target menu sele select ct New targe targett Configu Configurati ration on Target -> New target Configuration. It will will open a window window like given given below. below.
Specif Specify y any arbitr arbitrary ary targe targett name. name. For For Eg., Eg., 6748co 6748confi nfig.c g.ccxm cxmll ( Extension should be .ccxml).. Click Finish Finish then then you will will get get config configura uratio tion n wind window ow for the create created d targ target et .
B. 1. Sele Select ct the the Connection: Texas instruments XDS100v1 USB Emulator 2. Sele Select ct the the Device: TMS320C6748. To make make sear search ch easi easier er type type 6748 in devi device ce bloc block. k. *
Check the box TMS320C674 8 and finally Save. Click Save
Dept of ECE
Page | 48
C. Next
goto goto Adva Advanc nced ed tab give give the the suita suitable ble Gel file file path path as shown shown bel below. ow.
Follow the path. “C:\6748support\c6748.gel”
D. Go to to view view option option and and selec selectt the Targe Targett Config Configura uratio tion: n: View->Target Configuration. A wiza wizard rd will will open open in the the work worksp spac acee expa expand nd the the User User Defi Define ned d folde folderr and and you can can find find your our target target,, Right Right click click on 6748confi 6748config.c g.ccxml cxml and selec selectt the Launc Launch h Selecte Selected d Configurat Configuration. ion.
E. Conn Connec ectt CCS CCS to the the targ arget Goto Target -> Connect Target
Now our targe targett is succe successf ssfull ully y confi configu gure red d and conne connect cted. ed. In futu future re we no need need to repe repeat at thes thesee step steps. s. If we are are work workin ing g with with the the same same hard hardwa ware re we can just open the already already configur configured ed target target and launch launch the selecte selected d configur configuratio ation n and conne connect ct it.
Dept of ECE
Page | 49
5. Creating the New Project Step P1: Chang Changee the Perspe Perspecti ctive ve Debug to C/C++ from from the the rig right corn corner er of the the CCS CCS Step P2: Go to to File File New CCS Proje Project. ct.
Step P3: Spec Specify ify the the name name of the the proj projec ectt in the the space space provided provided .and .and Click Next
Eg., Eg., Projec Projectt Nam Name: e: Hello Hello World World
Dept of ECE
Page | 50
Select the project type
Projec Projectt Type: Type: C6000 C6000 Next Click Next
*Howev *However er our our targe targett is based based on on C6000 C6000 family family,, Based Based on the the fami family ly we need need to selec selectt the Project Type. Set Set the the project settings window window as shown shown below. below.
Click finish
Dept of ECE
Page | 51
6. To writ writee the the prog progra ram m sele select ct one one new new file file.. A.Go to to File File New Sourc Sourcee File. File.
B. Spec Specif ify y the the arbi arbitr trary ary sourc sourcee file file name name.. It should should be in the the sour source ce fold folder er (cur (curre rent nt proj projec ectt name.).
Note: Exte Extens nsio ion n of the the sour source ce file file must ust be the the lang langua uag ge what what we pref prefer erre red d to writ writee the the code code.. Eg: For c-> c-> .c C++ C+ + -> .cpp .cpp Assembly Assembly -> .asm
C.Type Type your C – Code
Dept of ECE
Page | 52
7. BUILD Build Build the the progr program am to check check your code code Have Have any any erro errors rs and and warning warnings. s. Go to Project
Build Build activ activee projec project. t.
If your our code ode doesn window that Build Comple Complete te “Build
t have have
any any erro errors rs and and warn warnin ings gs,, a me mess ssag agee will will be prin printe ted d in the the cons consol olee
for for the proje projecc t”
Problems window window displa display y err error orss or warnin warnings gs,, if any. any.
Dept of ECE
Page | 53
8. DEBUG Afte Afterr succ succes essf sful ul Buil Build, d, Debu Debug g your our code. code. Dur During ing this this step step ccs ccs will will conn connec ectt to targ target et and and it will load program program on to target. target. Debug Active Active proje project ct.. Goto Target Debug
During During debug debug ccs ccs will will displa display y follo followin wing g err error or messa message ge..
Now press press reset reset butto button n on the the 6748 hardwar hardwaree , then then click click on retry. retry.
Once Once you clic click k on on retry retry ccs ccs will will load load prog progra ram m on on to proc proces esso sorr and and the then n ccs ccs will will guid guidee us to debug debug mode, mode, if it is not not done automat automatical ically. ly. Chang Changee the Perspe Perspecti ctive ve C/C++ to Debug from from the the righ rightt cor corne nerr of the the CCS. CCS.
Dept of ECE
Page | 54
9.RUN Now Now you can can run run the the code code,, by sele select ctin ing g the the opt optio ion n run run fro from m the the dial dialog og box box else else you can can Go to Target Run
Once Once you run run the prog progra ram m the the outp output ut will will be prin printe ted d in the the Cons Consol olee
Dept of ECE
Page | 55
DSP Laboratory (10ECL57)
V Semester
KSSEM, Bangalore
PROGRAM 1 LINEAR CONVOLUTION OF TWO GIVEN SEQUENCES Aim: To write write the the C code code to perfo perform rm Linear Linear Convo Convoluti lution on upon upon two two given given disc discret retee time signals. C Code: /* PROGRAM PROGRAM TO TO IMPLEME IMPLEMENT NT LINEAR LINEAR CONVOLU CONVOLUTION TION */ #include int y[20]; main() {
in t m=6; int n= n=6;
/*Lengt h o f i/ p sample s sequence*/ /*Length of of im impulse re response Co Coefficients*/
int i=0,j; /*Input Signal Samples*/ int x[15]={1,2,3,4,5,6,0,0,0,0,0,0}; /*Impu /*Impulse lse Respon Response se Coeffi Coefficie cients nts*/ */ int h[15]={1,2,3,4,5,6,0,0,0,0,0,0}; /*Calculate Values*/ for(i=0;i
Dept of ECE
Page | 56
DSP Laboratory (10ECL57)
V Semester
KSSEM, Bangalore
OUTPUT: Sequence 1: 1 2 3 Sequence 2: 1 2 3 Convolved Sequence: 1 4 10
4
5
6
4
5
6
20
35
56
70
76
73
60
36
Procedure:
1. Open Open Code Code Compos Composer er Studio Studio v4. 2. Connec Connectt the the target target (step (step 4). 4). 3. Create Create new project project with with name as sine sine wave wave (step (step 5). 5). 4. Create Create new source source file file and type type the code code linear. linear.cc and save save it.(step it.(step 6) 5. Now perform perform steps steps 7, 8and 9.Once 9.Once you you run the program program you can can watch convolut convolution ion result in the console window. 6. To vie view w the the valu values es in in grap graph h a. Go to Tools Graph Single Time 7. Set the Graph Graph Prop Propert erties ies as show shown n Buffer size : 11 DSP Data Type : 16-Bit Unsigned Int Start Address: y Display Data Size : 11 8. Click : Ok 9. Note Note down down the output output wavefo waveform. rm.
NOTE : Follow the same procedure for the rest of the experiments.
Dept of ECE
Page | 57
DSP Laboratory (10ECL57)
V Semester
KSSEM, Bangalore
PROGRAM 2: CIRCULAR CONVOLUTION OF TWO GIVEN SEQUENCES Aim: To write write the C code to perform perform Circular Circular Convol Convolution ution upon two given discrete discrete time signals. C Code: /* PROGRAM TO IMPLEMENT CIRCULAR CONVOLUTION */
#include int m,n,x[30],h[30],y[30],i,j,temp[30],k,x2[30],a[30]; void main() { printf("Ente printf("Enter r the length of the first sequence\n") sequence\n"); ; scanf("%d",&m); printf("\nEnter the length of the second sequence\n"); scanf("%d",&n); printf("\nEnter the first sequence\n"); for(i=0;i
/* If If leng length th of of both both seq seque uenc nces es are are not not equ equal al */ */
{ if(m>n)
/* Pad the smaller sequence with zero */
{ for(i=n;i
/* folding h(n) to h(-n) */
a[j]=h[n-j];
Dept of ECE
Page | 58
DSP Laboratory (10ECL57)
V Semester
KSSEM, Bangalore
/*Circular convolution*/ for(i=0;i
OUTPUT: Enter the length of the first sequence 4 Ente En ter r the the le leng ngth th of th the e se seco cond nd se sequ quen ence ce
4 Enter the first sequence 2 1 2 1 Enter Enter the second second sequen sequence ce
1
2
3
4
Convolved Sequence:
14
16
Dept of ECE
14
16
Page | 59
DSP Laboratory (10ECL57)
V Semester
KSSEM, Bangalore
PROGRAM 3: COMPUTATION OF N- POINT DFT OF A GIVEN SEQUENCE Aim: Computation of N point DFT of a given given sequence. sequence. C Code: #include #include int N,k,n,i; float pi=3.1416,sumre=0, sumim=0,out_real[8]={0.0}, out_imag[8]={0.0}; int x[32]; void main(void) { printf("Ente printf("Enter r the length length of the the sequence\ sequence\n"); n"); scanf("%d",&N); printf("\nEnter the sequence\n"); for(i=0;i
Dept of ECE
Page | 60
DSP Laboratory (10ECL57)
V Semester
KSSEM, Bangalore
OUTPUT: Enter the length of the sequence 4 Enter the sequence 0 1 2 3 DFT of the Sequence 6.0 -2.0 + 2.0i -2.0 - 0.0i -2.0 - 2.0i
Dept of ECE
Page | 61
DSP Laboratory (10ECL57)
V Semester
KSSEM, Bangalore
PROGRAM 4: IMPULSE RESPONSE OF FIRST ORDER AND SECOND ORDER SYSTEM Aim: To Calculate the Impulse Response of the First Order and Second Order Systems. S ystems. System: 1. First Order System: C Code: #include #defin e Orde r 1 #defin e Le n 10
/*First Orde r Filter */ /*Lengt h o f Impuls e Response*/
float y[Len] = {0,0,0}, sum; main() { int j,k; floa fl oat t b[Or b[Orde der+ r+1] 1]={ ={1, 1, 3} 3}; ; /*In /* Inpu put t Coe Coeff ffic icie ient nts* s*/ / floa fl oat t a[ a[Or Orde der+ r+1] 1]={ ={1, 1, 0.2}; 0.2}; /* /*Ou Outp tput ut Co Coef effi fici cien ents ts*/ */ printf("Impulse Response Coefficients:\n"); for(j=0;j=0) sum=sum+(b[k]*y[j-k]); } if(j<=Order) y[j]=a[j]-sum; else y[j]=-sum; printf("Response [%d] = %f\n",j,y[j]); } }
OUTPUT: Impulse Response Coefficients: Response [0] = 1.000000 Respon Response se [1] [1] = -2.800 -2.800000 000 Response [2] = 8.400000 Respon Response se [3] [3] = -25.19 -25.19999 9999 9 Response [4] = 75.599998 Respon Response se [5] = -226.7 -226.7999 99988 88 Response [6] = 680.399963 Response Response [7] = -2041.19 -2041.199951 9951 Response Response [8] = 6123.599 6123.599609 609 Response Response [9] = -18370.7 -18370.79882 98828 8
Dept of ECE
Page | 62
DSP Laboratory (10ECL57)
V Semester
KSSEM, Bangalore
2. Seco Second nd Ord Order er Syst System em::
C Code: #include #defin e Orde r 2 #defin e Le n 10
/*Secon d Ord e r Filter*/ /*Lengt h o f Impuls e Response*/
float y[Len] = {0,0,0}, sum; main() { int j,k; floa fl oat t b[ b[Or Orde der+ r+1] 1]={ ={1, 1, 3, 3, -0 -0.1 .12} 2}; ; /* /*In Inpu put t Coef Coeffi fici cien ents ts*/ */ float float a[Order a[Order+1] +1]={1 ={1, , 0.2, -1.5}; -1.5}; /*Outp /*Output ut Coeffic Coefficien ients* ts*/ / printf("Impulse Response Coefficients:\n"); for(j=0;j=0) sum=sum+(b[k]*y[j-k]); } if(j<=Order) y[j]=a[j]-sum; else y[j]=-sum; printf("Response [%d] = %f\n",j,y[j]); } }
OUTPUT: Impulse Response Coefficients: Response [0] = 1.000000 Respo Response nse [1] [1] = -2.800 -2.800000 000 Response [2] = 7.020000 Respo Response nse [3] [3] = -21.39 -21.39600 6000 0 Response [4] = 65.030403 Respo Response nse [5] [5] = -197.6 -197.6587 58722 22 Response [6] = 600.779785 Response Response [7] = -1826.05 -1826.058350 8350 Response [8] = 5550.268555 Response Response [9] = -16869.9 -16869.93359 33594 4
Dept of ECE
Page | 63
DSP Laboratory (10ECL57)
V Semester
KSSEM, Bangalore
SAMPLE VIVA QUESTIONS:
1. What is sampling theorem? 2. What do you you mean by process process of reconstruction? 3. What are techniques of reconstructions? 4. What do you mean Aliasing? What is the condition to avoid aliasing for sampling? 5. Write the conditions of sampling. 6. How many types of sampling there? 7. Explain the statementstatement- t= 0:0.000005:0.05 0:0.000005:0.05 8. In the above example what does colon (:) and semicolon (;) denote? 9. What is a) Undersampling Undersampling b) Nyquist Nyquist Plot c) Oversampling. Oversampling. 10. Write Write the MATLAB MATLAB program program for Oversa Oversampling mpling.. 11. What is the use of command ‘legend’? 12. Write the difference between built in function, plot and stem describe the function. 13. What is the function of built in function and subplot? 14. What is linear l inear convolution? 15. Explain how convolution syntax built in function works. 16. How to calculate the beginning beginning and end of the sequence for the two sided sided controlled output? 17. What is the total output length of linear convolution sum. 18. What is an LTI system? 19. Describe impulse response of a function. 20. What is the t he difference between convolution and filter? 21. Where to use command filter or impz, and what is the difference between these two? 22. What is the use of function command ‘deconv’? 23. What is the t he difference between linear and circular convolution? 24. What do you mean by statement subplot (3,3,1). 25. What do you mean by command “mod” and where it is used? 26. What do you mean by Autocorrelation and Crosscorrelation sequences? sequences? 27. What is the t he difference between Autocorrelatio and Crsscorrelation. 28. List all the properties of autocorrelation and Crosscorrelaion sequence. 29. Where we use the inbuilt function ‘xcorr’ and what is the purpose of using this function? 30. How to calculate calculate output of DFT DFT using MATLAB? 31. What do you mean by filtic command, explain. 32. How to calculate output length l ength of the linear and circular convolution. 33. What do you mean by built in function ‘fliplr’ and where we need to use this. 34. What is steady state response? 35. Which built in function is used to solve a given difference equation? 36. Explain the concept of difference equation. 37. Where DFT is used? Dept of ECE
Page | 64
DSP Laboratory (10ECL57)
V Semester
KSSEM, Bangalore
38. What is the t he difference between DFT and IDFT? 39. What do you mean by built in function ‘abs’ and where it is used? 40. What do you mean by phase spectrum and magnitude spectrum/ give comparison. 41. How to compute maximum length N for a circular convolution convolution using DFT and IDFT.(what is command). 42. Explain Explain the statement- y=x1.*x2 43. What is FIR and IIR filter filt er define, and distinguish between these two. 44. What is filter? 45. What is window method? How you you will design an FIR filter using window method. 46. What are low-pass and band-pass filter and what is the difference between these two? 47. Explain Explain the the comma command nd : N = ceil(6.6 ceil(6.6 *pi/tb) *pi/tb) 48. Write down commonly used window function characteristics. 49. What is the MATLAB MATLAB command command for Hamming Hamming window? window? Explain. Explain. 50. What do you mea by cut-off frequency? 51. What do you mean by command butter, cheby1? 52. Explain the command command in detaildetail- [N,wc]=buttord(2*fp/fs,2 [N,wc]=buttord(2*fp/fs,2*fstp/fs,rp,As) *fstp/fs,rp,As) 53. What is CCS? CCS? Explain in detail detail to execute a program program using CCS. 54. Why do we need of CCS? 55. How to execute a program using ‘dsk’ and ‘simulator’? 56. Which IC is used in CCS? Explain the dsk, dsp kit. 57. What do you mean by noise? 58. Explain the program for linear convolution for your given sequence. 59. 59. Why we are using command ‘float’ in CCS programs. 60. Where we use ‘float’ and where we use ‘int’? 61. Explain the commandcommand- i= (n-k) % N 62. Explain the entire CCS program in step by step of execution. 63. What is MATLAB? What does it stand for?
Dept of ECE
Page | 65