ES4C8 Signal and Image Processing
Lab 1 – Digital Signal Processing Laboratory
Objectives
The activities in this laboratory session enable you to:
1. Use MATLAB commands specific to Digital Signal Processing
2. Use Difference Equations to obtain the Impulse Response of a filter and calculate Inverse Z-Transforms and output sequences
3. Determine the Transfer Function of a filter from its Pole-Zero plot
4. Estimate the influence of poles and zeros on the Frequency Response magnitude/phase.
Resources
MATLAB v.2017b and MATLAB Signal Processing Toolbox R2017b
Notes
1. The work you do, the notes you take and the m-files and plots you save during the laboratory serve as course material. The coursework and the exam will contain theoretical and mathematical questions (not M-code) from the material covered in the laboratory (marked Q on the labsheet).
2. MATLAB code is optimised for the use of matrices, therefore some code uses matrices to replace the standard approach of execution loops; as a result, the default indexing for a matrix starts at 1(one), whereas the mathematics in the lectures assume 0(zero). ⧫
3. The MATLAB symbol for the imaginary number sqrt(-1) is i (not j as used in the lectures, although when you type it recognises it). For clarity, do not use i or j as the name for a variable or an index. ⧫
4. Pay attention to the dimensions of matrices and vectors and to the precision of variables you use, as there may be accumulated errors due to lack of precision. Use commands such as size(my_variable), isempty(my_variable) and whos to check.
5. To speed up the completion of exercises you need to copy and paste the code from the lab sheet.
6. Pay attention to the left quotation mark ‘ and the straight apostrophe ‘ characters when copying and pasting text. MATLAB recognises only the latter. The same applies for the modifier letter circumflex accent ˆ and the circumflex accent ^ , or the n-dash –, the m-dash—, the hyphen – and the minus sign -. ⧫
7. It is important to initialize all variables at the start of any iterative algorithms, since memory may contain prior random values when allocated by the processor to a variable.
8. It is important to write non executable comments in your code (in MATLAB preceded by the percent sign % character) for explaining the purpose of the code.
9. For reasons of brevity, the code relating to labelling the plots has been omitted from Exercise 2 onwards. You are expected, for the assignment, to include graph axes with units, titles and legends, as appropriate for your plots, as well as more explicative comments in the code.
10. Save the code for each exercise in a separate m-file (e.g. lab1e21.m) and each figure as an image file in Portable Network Graphics *.png format (e.g. lab1ex21.png).
11. To get a mark for the completion of the lab, please show a demonstrator your results in the form of saved figures (*.png files) with your ID number and written answers for the questions Q 2.1, Q 2.3b, Q2. 4c and Q 3b. This mark represents 2% of your assignment mark.
ES4C8 Lab1 2018-19.docx 1 of 6 Copyright © 2015 School of Engineering
Exercise 1. Digital Signals
1.1. A question of precision
Consider the following hypothetical MATLAB code, in which we take the square root of a number x for N times iteratively and then square the result N times iteratively [1]. With reference to Notes 4 and 5, debug the code and explain the result. Use the command whos to check for precision ⧫.
x = 11; %choose the number
N = 50; %choose the number of iterations %calculate the square root iteratively for k = 1:N
x = sqrt (x);
fprintf (1, ‘ %d %.20f \n ‘, k, x); end
%retrieve the original number by calculating the square iteratively
for k = 1:N
x = x ˆ 2;
fprintf (1, ‘ %d %.20f \n ‘, k, x); end
1.2. Sampled signals
Type in the following code and observe the commands necessary to plot different types of signals (‘continuous’, sampled and zero order hold – ZOH approximation), as well as the numbering of the x-axis. Make a note of the axes labelling, the plot title commands and the positioning of the ID number in the figure, which you need to include in further answers. Use the File>>Save As… from the figure menu button to save your figure as a png image.
omega = pi; %real frequency[rad/s]=>f=0.5[Hz], T=2[s]
t = [0:1/20:5]; %time length of signal[s], stepped in small increments s = cos (omega*t); %‘continuous’ signal
subplot(2,2,1)
plot (t,s);
xlabel(‘time[s]’); ylabel(‘Amplitude’); title(‘Continuous signal’);
grid on;
Ts = 0.2; %sampling period[s/sample]=>sampling frequency fs=5[Hz],T/Ts=10 n=[0:20*Ts:length(s)]; %sampling positions
x = s(n+1); %
subplot(2,2,2)
stem (x); %use ‘stem’to plot individual samples
xlabel(‘Sample no.(Matlab index)’); ylabel(‘Amplitude’); title(‘Sampled signal’);
grid on;
n_s =[0: length(x)-1]; %true sample number subplot(2,2,3)
stem (n_s,x); %see note 2 ⧫
xlabel(‘Sample no. ‘); ylabel(‘Amplitude’); title(‘Sampled signal’); grid on;
subplot (2,2,4) stairs (n_s,x);
xlabel(‘Sample no.’); ylabel(‘Amplitude’); title(‘ZOH signal’);
grid on;
gtext (‘ID uXXXXXXX’);%change the ID uXXXXXXX by your own ID number
Q 1.2: What does ‘continuous’ mean in this instance?
ES4C8 Lab1 2018-19.docx 2 of 6 Copyright © 2015 School of Engineering
Exercise 2. Linear Shift Invariant (LSI) Systems
2.1. Difference Equation
Consider the Difference Equation 𝑦(𝑛) = 0.9𝑥(𝑛) + 0.8𝑥(𝑛 − 1) − 0.5𝑦(𝑛 − 1) describing an LSI system. A simple realisation of the Impulse Response is given by the output sequence 𝑦(𝑛) this system when the input sequence 𝑥(𝑛) is the unit impulse {1, 0, 0, … , 0, … }.
NT = 10;
x = zeros (NT, 1); y = zeros (NT, 1); x(1) = 1; forn=1:NT
%choose the number of samples
%initialise the memory for the input sequence %initialise the memory for the output sequence %input signal is a unit impulse; see note 2 ⧫ %seenote2 ⧫
y(n) = 0.9 * x(n); if (n > 1)
y(n) = y(n) + 0.8 * x(n-1)- 0.5 * y(n-1);
end end
figure; stem (y); %output signal is the Impulse Response h(n)
The more generic code below is for the generalised Difference Equation
𝑦(𝑛)=[𝑎0 𝑥(𝑛)+ 𝑎1𝑥(𝑛 −1)+⋯+ 𝑎𝑁𝑥(𝑛−𝑁)]−[𝑏1𝑦(𝑛 −1)+⋯+ 𝑏𝑀𝑦(𝑛−𝑀)]. Complete the code to implement the filter above, by giving appropriate values to the multipliers 𝑎 and 𝑏.
NT = 10;
% insert here the values for multipliers ak
a=
N = length (a)- 1;
% insert here the values for multipliers bk;
% pay attention to the sign of the feedback path filter coefficients b=
M = length (b)-1; x =
y =
x(1) =
for n = 1:NT y(n) = 0;
%initialise the memory for the input sequence
%initialise the memory for the output sequence
%initialise input signal as a unit impulse
for k = 0:N
if ((n-k) > 0)
y(n) = y(n) + x(n-k)* a(k+1);
end end
for k = 1:M
if ((n-k) > 0)
end end
end
… …
y(n) = y(n) – y(n-k) * b(k+1);
% insert here a plotting command, axes and title
% insert here your ID number on the graph
Q 2.1: What are the first 5 values for the filter IR, namely {𝒉(𝟎), 𝒉(𝟏), 𝒉(𝟐), 𝒉(𝟑), 𝒉(𝟒)} ? ES4C8 Lab1 2018-19.docx 3 of 6 Copyright © 2015 School of Engineering
An alternative MATLAB implementation of the filter described by the Difference Equation above is done using the function filter().
NT = 10;
a = [0.9 0.8];
b = [1 0.5]; % b(1) is always 1; see note2 ⧫
x = zeros (NT, 1);
x(1) = 1;
y = filter (a, b, x); % y has the same number of samples as x help filter % check the implementation of the filter command
2.2. Linearity, Superposition and Time Invariance
Now consider the sequence 𝑥(𝑛) = {0.6, −0.2, 0.8} as the input to the filter in Section 2.1.
Q 2.2: Using the properties of an LSI system, prove mathematically and computationally that the output of the system is, in fact, the summation of the individual, scaled and delayed, impulse responses of each individual element of 𝐱(𝐧).
2.3. Z-Transform and Transfer Function
𝑍𝑇 𝑍𝑇
Using the shift property of the Z-Transform (if 𝑥(𝑛) ⇔ 𝑋(𝑧), then 𝑥(𝑛 − 𝑚) ⇔ 𝑧−𝑚𝑋(𝑧)), the Z-
Transform of the generalised Difference Equation becomes
𝑌(𝑧) = [𝑎0 𝑋(𝑧) + 𝑎1𝑧−1𝑋(𝑧) + ⋯ + 𝑎𝑁𝑧−𝑁𝑋(𝑧)] − [𝑏1𝑧−1𝑌(𝑧) + ⋯ + 𝑏𝑀𝑧−𝑀𝑌(𝑧)] or 𝑌(𝑧) =
𝑋(𝑧)
𝑎0+ 𝑎1𝑧−1+⋯+ 𝑎𝑁𝑧−𝑁 = 𝐻(𝑧), where 𝐻(𝑧) = ∑∞ h(𝑛)𝑧−𝑛 is the transfer function. 1+𝑏1𝑧−1+⋯+ 𝑏𝑀𝑧−𝑀 𝑛=0
Q 2.3a: Compute the first 6 values for the filter IR coefficients 𝒉(𝒏) and state if the filter is FIR or IIR for the filter given by the TF 𝐻(𝑧)=0.9+0.8𝑧−1−0.4𝑧−2+0.1𝑧−3 . Calculate the
1−𝑧−1
coefficients by polynomial division and prove they are the same.
We can consider the previous expression as 𝐻(𝑧) = 𝐻1(𝑧)𝑋1(𝑧), where
𝐻 (𝑧) = 0.9+ 0.8𝑧−1 − 0.4𝑧−2 + 0.1𝑧−3 and 𝑋 (𝑧) = 1 , in other words the output of the
1 1 1−𝑧−1
filter given by 𝐻1(𝑧) when the input is the unit step sequence 𝑥1(𝑛), ie {1, 1, 1, . . . , 1, . . . } ⇔ 1−𝑧−1.
Complete the code below to include the new filter and the new input sequence.
Q 2.3b: Prove that the output sequence is the same as in Q2.3a.
𝑍𝑇 1
NT = 10;
a = % insert here the values for multipliers ak b = % insert here the values for multipliers bk x = ones (NT, 1); %input signal is a unit step;
y=
ES4C8 Lab1 2018-19.docx 4 of 6 Copyright © 2015 School of Engineering
2.4. Transfer Function, Poles and Zeros
Consider a system with a single, real, pole given by the TF (𝑧) = 𝑧 . For this recursive system 𝑧−𝑝
parameter p controls how much of the previous output is fed back and determines the stability of the system. The code below calculates the filter IR coefficients and shows the Plot-Zero plot for 𝑝 = 0.9 .
NT = 40;
z = 0;
p = 0.9;
a = [1 0]; b = [1 -p];
x = zeros (NT, 1); x(1) = 1; y = filter (a, b, x); figure;
subplot(4,2,1); stem(y)
subplot(4,2,2);zplane(z,p)
%- plot a unit circle, hold on, then plot the pole coordinates)
%only in the DSP System Toolbox (alternative
Q 2.4a: Complete/modify the code above to display on the same graph results for 𝑝 = {0.6, 1.1, −0.8, −1.2}. Comment on the influence of the pole on the filter IR and its stability.
Q 2.4b: Complete/modify the code above to display the graph results for the TF (𝑧) = 1 , 𝑧−𝑝
for 𝑝 = {0.6, 1.1, −0.8, −1.2}. Comment on the influence in terms of the zero (in this case the missing 𝑧 of the numerator) on the filter IR and its stability.
Now consider a system with complex poles. Because the coefficients of the time domain terms in the difference equation must be real, the coefficients of 𝑧 , when written out as a polynomial, must in turn be real. This in turn means that the roots of the numerator and denominator polynomials in 𝑧 (the zeros and poles respectively) must be complex conjugates.
The code below calculates the TF coefficients for 𝑝 = 0.9𝑒𝑗𝜋/10. Defining 𝑝 as the value of the pole location, the poly()function may be used as shown below to expand the factored equation to obtain the denominator coefficients.
10 10 4 20
Q 2.4c: Complete/modify the code above to display results for 𝑝 = {𝑒𝑗𝜋, 1.2𝑒𝑗𝜋, 𝑒𝑗𝜋, 𝑒𝑗𝜋} and
its complex conjugate 𝑝∗ for the TF 𝐻(𝑧) = 𝑧2 . Comment on the influence of the pole (𝑧−𝑝)(𝑧−𝑝∗)
on the filter IR and its stability. Calculate, then verify on the plot, the frequency of the
𝑜𝑠𝑐 𝒐𝒔𝒄 oscillator relative to the sampling frequency 𝑓 or 𝝎 .
p = 0.9*exp(j*pi/10);
a = [1 0 0];
b = poly([p conj(p)]); %expand (z-p)(z-p*)
roots(b) %the pole location can be checked
ES4C8 Lab1 2018-19.docx 5 of 6 Copyright © 2015 School of Engineering
𝑓𝝎 𝑠𝒔
Exercise 3. Frequency Response
The relationship between the discrete frequency Ω and the real frequency 𝜔 of a signal is Ω = 𝜔𝑇
𝑠 The frequency response is equal to the TF on the unit circle, 𝐻(Ω) = 𝐻(𝑒𝑗Ω) = 𝐻(𝑧)|𝑧=𝑒𝑗Ω . The
[rad/sample]=[rad/s][s/sample], where 𝑇 is the sampling interval. 𝑠
code below calculates the magnitude and phase of the FR of the filter given by the TF 𝐻(z) = 1 , when 𝑟 =0.9 and 𝜃 = 𝜋/10 .
1−2𝑟 cos(𝜃)𝑧−1+𝑟2𝑧−2
r = 0.9; % pole radius theta = pi/10; % pole angle %coefficients of z in the TF
a = [1 0 0];
b = [1 -2*r*cos(theta) r*r];
%frequency response
Omega=0:pi/100:2*pi; %discrete frequency range L = length(Omega);
H = zeros(L,1); for k =1:L
expvec = exp(j*[2:-1:0]*Omega(k)); %current value for z=exp(j*Omega)
H(k) = sum(a.*expvec)./ sum(b.*expvec); end
subplot(1,2,1); plot (Omega, abs(H)); grid on; %magnitude response (MR) set( gca , ‘xlim’ , [0 2*pi]); set( gca , ‘xtick’ , [0:4]*pi/2); subplot(1,2,2); plot (Omega, angle(H)); grid on; %phase response (PR) set( gca, ‘xlim’, [0 2*pi]);set( gca,’xtick’,[0:4]*pi/2);
Q 3a: Modify the code above to display results for 𝑟 = {0.5, 0.1} and 𝜃 = {𝜋/4, 𝜋/20 } . Comment on the influence of the magnitude and phase of the FR.
An alternative MATLAB implementation of the FR is done using the function freqz(). This function is implemented using Fast Fourier Transform (FFT algorithms).
r = 0.9; % pole radius
theta = pi/10;
a = [1 0 0];
b = [1 -2*r*cos(theta) r*r];
[H,Omega] = freqz(a,b,100); %only in the DSP System Toolbox fs = 1000; %sampling frequency [Hz]
[H1,omega] = freqz(a,b,100,fs);
… % insert here plotting commands for MR and PR % do not forget axes, title and ID number
Q 3b: Modify the code above to represent the FR over the entire frequency range [0:2𝜋]
taking into account the symmetry properties of the FFT. Calculate the true frequency of the
poles when the sampling frequency is 𝑓 = 1000 Hz. 𝑠
References: 1. Leis J. W., Digital Signal Processing Using MATLAB for Students and Researchers, 2011 (ebook)
ES4C8 Lab1 2018-19.docx 6 of 6 Copyright © 2015 School of Engineering