INTRODUCTION TO MATLAB FOR DIGITAL SIGNAL PROCESSING
Objective: To learn the basics of the MATLAB software for DSP applications.
System Requirement: Any IBM PC compatible MATLAB software
Background information:
MATLAB is both a computer programming language, and an environment for simulating,
analysing and designing signals and systems. ‘MATLAB’ is a contraction of ‘MATrix
LABoratory’. The basic entity in MATLAB is a matrix, an m �n array of variables, written in
general as:
�a11 a12 K a1n �
�a a21 K a2 n �
A = �21 � or A=�
aij �
� � where 1 �i �m and 1 �j �n.
�M M O M�
� �
am1 am 2
� K amn �
That is, i is the row-index and j is the column-index.
MATLAB uses a computing environment for numerical computation and visualisation. It
integrates numerical analysis, matrix computation, signal processing, and graphical display in
an intuitive and easily-used environment. There are three main windows in the MATLAB
environment: the “command” window, the “program” window, and the “graphics” window.
MATLAB can be operated in two ways: in the interpreter mode and in the compiler mode.
In the interpreter mode, a single calculation or command, typed in the command window, can
be executed. This is sometimes useful, for example, when you want to use the command
window as a calculator, or when you want to investigate array values whilst debugging a
program. There is no way to store programs in this mode, although there is limited access to
previous commands via the up-arrow (↑) key. In the compiler mode, sequences of MATLAB
commands are assembled in an “M-file”, which is a MATLAB stored program. The program
can be executed in the usual way. Creating, editing and debugging M-files in the program
window in compiler mode is the customary mode of operation.
MATLAB also features sets of application-specific functions called toolboxes. Toolboxes
are comprehensive collections of MATLAB M-files that extend the MATLAB environment
to solve particular classes of problems, including Digital Signal Processing.
Procedure 1: Interpreter Mode
1.1. Single command
In the MATLAB window, simple calculations (or operations) can be executed by
typing the commands directly. For example, you can type:
>>15+35-40*30/50
>>10^2-9\45
What did you observe in the MATLAB window? You can see that MATLAB can be
used as a calculator!
_____________________________________________________________________
S. Sanei/K. Lever MATLAB - 1 October 2004
1.2. Creating variables
To clear all the variables in the workspace, type:
>>clear
To list out all the variables in the workspace, type:
>>who
You can see that there is now no variable in the workspace.
Now, type:
>>A=6
What response did you get?
_____________________________________________________________________
Now, again type:
>>who
What response did you get now?
_____________________________________________________________________
In MATLAB, all variables are matrix arrays (simple scalar variables, and row and
column vectors can be thought of as special cases of a matrix).
Type:
>>B=[1 2 3; 4 5 6]
1 2 3�
�
The MATLAB window will show the matrix B: B = �
4 5 6�
� �
Type:
>>C=[1 2 3; 4 5 6];
Now C has the same value as B. However, with “;” added at the end of the command
line, MATLAB no longer echoes what you typed in the command window.
At this point, can you list out all the variables that you have in your workspace?
Which MATLAB command should you use?
_____________________________________________________________________
To see the contents of the variables, type their names, and record what you observe.
>>A
>>B
>>C
_____________________________________________________________________
What would happen (and why) if you were to type the following command? Try it.
S. Sanei/K. Lever MATLAB - 2 October 2004
>>A;
_____________________________________________________________________
Type:
>>C=[3 4 5 6 7 8 9]
As you can see, C is assigned a new value.
Type:
>>D=[3:1:9]
Is D identical to C?
_____________________________________________________________________
In the expression, [3:1:9], the “3” represents the first element of a vector, the “1” is
the incremental step, and the “9” is the last element of the vector. Now type:
>>E=[9:-1:3]
What is the difference between D and E?
____________________________________________________________________
Type:
>>F=[9:-2:3]
What is the difference between E and F?
_____________________________________________________________________
You can also extract an element or a sub-matrix from a matrix. Type the following:
>>D(3)
>>D(2:4)
>>B(2,2)
>>B(1,2:3)
>>B(1,:)
>>B(2,:)
>>B(:,1)
>>B(:,3)
Explain how each of these values is obtained. What does “:” do in this context?
_________________________________________________________________
1-3. Basic calculations - Matrix multiplication
Type:
>>G=A*8;
Write down the value of G. Can you explain how this value comes about?
_____________________________________________________________________
Type:
S. Sanei/K. Lever MATLAB - 3 October 2004
>>D
>>D’
What is the difference between D and D'?
_____________________________________________________________________
Type:
>>D*D'
Write down the answer. Can you explain how this value comes about?
_____________________________________________________________________
Type:
>>TMP=[2 4 6; 3 5 7]'
What is the value of this variable, “TMP”?
_____________________________________________________________________
Type:
>>B*TMP
What is the answer? Verify the above result by hand.
_____________________________________________________________________
Note that in matrix multiplication, the dimensions of the two matrices should always
fit each other, i.e. mn followed by nm.
You can also find the inverse of a square matrix with the help of MATLAB
Type:
>>H=[2 4; 3 1]
>>INVH=inv(H)
What is the answer? Verify your result by hand.
_____________________________________________________________________
Type:
>>I=H*INVH
What is the answer? Explain the result.
_____________________________________________________________________
S. Sanei/K. Lever MATLAB - 4 October 2004
1-4. Basic calculations – Element-wise multiplication and division
In this mode, two matrices with the same dimension (unless one of them is scalar) are
multiplied together element by element. For example, [1 2].*[3 4]=[3 8].
To illustrate the effect of element-wise multiplication, Type:
>>B*TMP
>>B*TMP'
>>B.*TMP
>>B.*TMP'
Notice that two of the above expressions result in an error. Why?
_____________________________________________________________________
For the other two expressions, can you verify your result by hand?
_____________________________________________________________________
Why do you need to transpose matrix “TMP” when doing element-wise
multiplication?
_____________________________________________________________________
Type:
>>B./TMP'
What is the answer? Verify your result by hand.
_____________________________________________________________________
Note that addition and subtraction are always element-wise (unlike multiplication and
division).
1.5. Built-in MATLAB functions
There are many built-in functions in the MATLAB library. For instance; we can
simply call the "roots" function to calculate the roots of polynomials, such as
x 2 3x 2 = 0 .
Type:
>>A=[1 3 2];
>>roots(A)
What does the variable A represent?
_____________________________________________________________________
What are the roots of the above equation? Verify your result by hand.
_____________________________________________________________________
S. Sanei/K. Lever MATLAB - 5 October 2004
A short description is usually given for each built-in function. To access that
description, we can use the command “help”:
>>help roots
What happens?
_____________________________________________________________________
Some other builtin functions can help you extract or insert numbers in a vector.
Type:
>>J=[ones(1,4),[2:2:11],zeros(1,3)]
>>length(J)
>>J(2:2:length(J))
Explain the results you obtained at each stage.
_____________________________________________________________________
Many tools for drawing graphics are also available. As an example, see how a simple
threedimensional graph can be displayed.
Type and watch:
>>[X,Y,Z]=sphere(20);
>>mesh(X,Y,Z)
>>hold
>>mesh(X+10,Y+10,Z)
Explain what “hold” does.
Procedure 2: Compiler Mode
2-1. Simple M-files
An ".m" file (or M-file) is a collection of MATLAB commands and can be executed
in the MATLAB program window. To create an ".m" file, go to the "File" menu and,
select "New". A program editor window is opened. In this window type in the
following commands:
t=[1:1:100];
s=sin(4*pi*t/100);
plot(t,s);
Go to the "File" menu, select "Save as" and type in the name of the file as “testmat”
(or any other filename that you wish). This file will be saved as "testmat.m". In the
MATLAB command window, type:
>>testmat
S. Sanei/K. Lever MATLAB - 6 October 2004
What happened? Can you explain what the program does?
_____________________________________________________________________
2-2. Loops
In MATLAB programs (M-files), you can also construct loops. Try the following
program (you may name it testloop.m), and describe what happens:
A=[1 3 2; 1 –3 2; 2 5 2];
for index=1:3,
Aroots=roots(A(index,:))
end
_____________________________________________________________________
2-3. More Practice
Study the following program to work out what it does, and try it:
t=[1:1:100];
s=sin(4*pi*t/100);
subplot(2,1,1),
plot(t,s);
axis([0 100 -1.5 1.5]);
grid on;
rec_s=abs(s);
subplot(2,1,2),
plot(t,rec_s);
axis([0 100 -1.5 1.5]);
grid on;
What waveform is signal "rec_s?" What do the commands "subplot," "axis," "grid on"
do? (Hint: use online help).
_____________________________________________________________________
Try a second example:
t=[1:1:100];
s=sin(4*pi*t/100);
subplot(2,1,1),
plot(t,s);
axis([0 100 -1.5 1.5]);
grid on;
doub_s=s.*s;
subplot(2,1,2),
plot(t,doub_s);
axis([0 100 -1.5 1.5]);
grid on;
S. Sanei/K. Lever MATLAB - 7 October 2004
What waveform is signal "doub_s?" Is it a sine wave? Why is the frequency of
“doub_s” twice the frequency of s? Can you derive that?
_____________________________________________________________________
Try another example: write and run an Mfile as follows:
t=-2:0.05:3;
x=sin(2*pi*0.789*t)
plot(t,x), grid on
title(‘Test plot of Sinusoid’)
xlabel(‘Time(sec)’)
ylabel(‘Amplitude(volt)’)
Now amend the program to plot a new signal y = 0.75cos(20.789t) on top of the
previous one (type “help plot” to see the plot format to find out how to show multiple
signals in the same frame). Use the “save and run” ikon in the toolbar when amending
existing programs.
Try the following M-file
step = 0.1; u = [ ];v = [ ]; w = [ ]; s=0;
while( length(u)<=300)
u = [u;s];
v = [v;s^2];
w = [w;sin(s)];
s = s+step;
end
[u,v,w]
plot3(u,v,w,'blue')
grid
What does the program do?
Exercise 1
Plot a signal u(t) which is the sum of a sine signal x 1(t), with amplitude of 0.8 volt,
frequency of 1000 Hz and phase shift of zero and another sine wave x2(t) with
amplitude of 0.6 volt, frequency of 2500 Hz and phase shift of 45 {i.e. u(t) = x1(t) +
x2(t)}. Hint: Use time index, t=[0:1:39]/20000.
Exercise 2
Plot the magnitude and phase of the following expression for 0 2:
H = 1 0.3e j
S. Sanei/K. Lever MATLAB - 8 October 2004
Note that MATLAB® has two built-in constants, namely: j and pi where j = 1
and pi = 3.1415926 … . Some other MATLAB® functions that you may use here are:
exp, abs and angle. Use “subplot” to show both plots at the same time. Use the
“help” command to learn more about all these functions.
2-4. Functions
M-files can be either scripts or functions. Scripts are simply program files containing a
sequence of MATLAB statements. Functions (called from a script or another function) make
use of their own local variables, accept input arguments and deliver output values back to the
point where they were called. A line at the top of a function M-file contains the syntax
definition. The name of a function, as defined in the first executable line of the M-file, should
be the same as the name of the file without the “.m” extension.
Example:
The following script, “mainFIR.m”, is a “main program”, and calls the specially-written (not
built-in) functions “delta.m” and “lp2hp.m”. Functions can call other functions, (both built-in
and specially-written) and in this example “lp2hp.m” calls “delta.m”. Check that functions
“delta.m” and “lp2hp.m” are not built-in functions, by typing “help delta” and “help lp2hp”.
The purpose of the main program is to design an FIR linear phase filter, using the built-in
function “fir1.m” (type “help fir1” in the command window to understand the syntax of the
fir1.m function). The impulse response and the amplitude of the frequency response of this
lowpass filter are then evaluated and plotted. The linear phase lowpass filter is then converted
to a linear phase highpass filter by means of a simple transformation between the impulse
responses. In an obvious notation: hH [n] = d [n M ] hL [ n] , where M is the midpoint of the
impulse responses. (This only works when the common length of the impulse responses is
odd, and when the lowpass filter has 0 dB gain at DC: why?).
The main program and the two functions need to be in the same directory, where they can
“see one another”. Note the use of non-executable comments (% … ) to annotate the code:
these should be plentiful enough to make the code self-explanatory long after it was first
written.
% mainFIR.m
% This main program designs a LPF,
% converts it to a HPF, and
% plots impulse and frequency
% responses of both filters.
clear;clc;
for k = 1:2
figure(k)
clf
end
% Choose signal duration, and
% set up time and frequency axes
N = 512; % duration
t = [0:N-1]; % time axis
S. Sanei/K. Lever MATLAB - 9 October 2004
f = [-N/2:N/2-1]./N;% frequency axis
% Set up impulse at origin
% using "own" function
d = delta(1,N); % NB: no zero index in MATLAB
% Design a 21-tap (order 20)
% FIR LPF filter
% using built-in function fir1(.,.).
% The cutoff frequency is 0.2fs (NOT 0.4fs).
b = fir1(20,2.*0.2); % feedforward coeffs.
Lb = length(b);
a = [1]; % no feedback coeffs for FIR filter
% Find impulse and frequency response of LPF
hL = filter(b,a,d); % impulse response
AL = 20.*log10(fftshift(abs(fft(hL))));
% Plot LPF characteristics
figure(1)
subplot(2,1,1);
stem(t(1:Lb+1),hL(1:Lb+1),'k.')
xlabel('time, n')
ylabel('h_L[n]')
title('Impulse and frequency response of LPF')
subplot(2,1,2);
plot(f,AL,'k-')
axis([min(f) max(f) -80 10]);
xlabel('normalised frequency, f')
ylabel('A_L[n], dB')
grid on
% Find impulse and frequency response of
% corresponding HPF using "own" function
hH = lp2hp(Lb,hL);
AH = 20.*log10(fftshift(abs(fft(hH))));
% Plot HPF characteristics
figure(2)
subplot(2,1,1);
stem(t(1:Lb+1),hH(1:Lb+1),'k.')
xlabel('time, n')
ylabel('h_H[n]')
title('Impulse and frequency response of HPF')
subplot(2,1,2);
plot(f,AH,'k-')
axis([min(f) max(f) -80 10]);
xlabel('normalised frequency, f')
ylabel('A_H[n], dB')
grid on
S. Sanei/K. Lever MATLAB - 10 October 2004
% delta.m
% Set up a delta-function impulse
function d = delta(location,duration)
d = zeros(1,duration);
d(location) = 1;
% lp2hp.m
% Convert LPF to HPF by subtracting
% LPF impulse response from a unit
% impulse at its midpoint.
function hH = lp2hp(Lb,hL)
LhL = length(hL);
if rem(Lb,2) == 0
disp('Does not work for even length')
else
midpoint = (Lb+1)/2;
end
imp = delta(midpoint,LhL);
hH = imp - hL;
Exercise 3
Show, by adding further commands to the “mainFIR.m” program, that direct design of
a highpass filter using the fir1.m function gives rise to a better frequency response.
This shows that the fir1.m function definitely does not use the lowpass-to-highpass
transformation technique implemented in “lp2hp.m”.
Exercise 4
By adding further commands to the “mainFIR” program, generate and plot (in
additional figures) the phase versus frequency characteristics of each of the two
filters. Remember to modify the “housekeeping” commands at the top of the main
program, used to clear the figures before program execution. Are the filters linear
phase? If so, why?
It is not necessary, from the point of view of understanding how specially-written functions
work, to understand the detail of the use of the filter design (fir1.m) and fast Fourier
transform (fft.m and fftshift.m): the user can still run the main program and check that the
results make sense. For learning the basics of MATLAB , the main point is to notice how the
functions are called (indicated by red in the code given above), and how the syntax of each
function is defined in the function header.
On the other hand, for those who are also learning the digital signal processing applications
of MATLAB, this example is a good illustration of some fundamental DSP concepts: filter
design and evaluation in the time-domain and frequency-domain.
S. Sanei/K. Lever MATLAB - 11 October 2004
