A Matlab Tutorial

An Introduction to Matlab
(index)

Matlab stands for Matrix­laboratory. It is an interactive 
program that provides easy access to many of the most 
widely used functions for matrix analysis and digital 
signal processing. 

Due to its large number of functions for manipulating 
matrices,

and the digital image processing toolbox, Matlab is an 
excellent enviroment for developing prototype systems for 
digital image processing.

Matlab provides excellent tutorials that are accessible by 
typing

>>demo

The Basic matrix operations tutorial under the Matrices 
tutorial,

the Image Processing and Signal Processing tutorial under 
Toolboxes are highly recommended.
To get information on a particular function of Matlab, we 
type

>>help function_name

Defining Matrices and Functions in 
Matlab
(index)

In Matlab, we can define an array by simply listing its 
elements and separating out each row using ;

For example, the matrix A which is mathematically 
defined by

is described in Matlab by

» a=[1 2 3; 4 5 6; 7 8 9];

To view the contents of the Matrix, we simply type the
name of the Matrix, without inserting a ; at the end

» a
 a = 1 2 3

4 5 6

7 8 9

Comments in Matlab must be preceeded by %

>>% This is a comment in Matlab

To define a function in Matlab, we first create a function 
with the name of the function. We then define the function 
in the file. For example, the Matlab documentation 
recommends stat.m written as:

function [mean,stdev] = stat(x)

% This comment is printed out with help stat

n = length(x); % assumes that x is a vector.

mean = sum(x) / n; % sum adds up all elements.

stdev = sqrt(sum((x ­ mean).^2)/n);

Viewing Matlab Matrices as

One Dimensional Arrays
(index)
 Convert a two­dimensional to a one­dimensional array
by simply using array_name(:)

» a(:)

ans =

Note that the elements in the one­dimensional array are 
arranged column by column, not row by row. 

This is the same convention that is used for two 
dimensional arrays in Fortran, and this will also be the 
convention that we will adopt for representing two 
dimensional arrays in C++.

Displaying and Printing in Matlab
(index)

To plot an array or display an image, select the figure 
using:

» figure(1); % select the figure to display.

» clf % clear the figure from any previous commands.

and use plot for displaying a one­dimensional array:

» plot ([3 4 5 3 2 2]); % plots broken­line graph.

or use imagesc for displaying a two­dimensional array 
(image):

» A=[1 2 3; 4 5 6; 7 8 9]; % define the matrix to 
display.

» imagesc(A); % display the image. 
 » axis image % fix aspect ratio in display.

» colormap(gray) % use gray for black and white.

» axis off % turn­off display of each pixel.

We can then print the figure using

» print filename ­f1 

which saves Figure 1 in the postscript file 
filename.ps.

Matlab Colon Operations 
(index)

In Matlab, we use the colon : to specify different elements 
in an array. For example

» 1:2:7

ans = 1 3 5 7
For one­dimensional arrays, we have:

» a=[8 12 14 32 45];

» a(1:2:5)

ans = 8 14 45

We can also access the elements in the array randomly by 
using an index array to access the elements:

» a=[6 12 8 9 13 17 21 34];

» I=[5 6 3 2];

» a(I)

ans = 13 17 8 12

These basic methods for accessing elements in one­
dimensional arrays are compatible with BLAS 
specifications. Using colon operations, we avoid 
unnecessary loops and the Matlab programs are far more 
efficient. As we shall see, we also have similar colon 
operations for two dimensional arrays. Avoiding the use of 
for loops is the key to achieving great computational 
performance using Matlab.

Matlab Colon Operations
 for two­dimensional arrays
(index)

For two­dimensional arrays, we can use 

 one­dimensional sequential iterators

» a=[1 2 3 4; 5 6 7 8; 9 10 11 12];

» a(1:2:12)

ans = 1 9 6 3 11 8 ,

which is equivalent to converting a to its one­
dimensional form 

[1 5 9 2 6 10 3 7 11 4 8 12], 

and then returning the 1st, 3rd, ..., 11th elements.

 one­dimensional random access iterator (using an 
index array)

» I=[1 5 8 10]; % indices apply to one­dimensional array 
a(:)

» a(I)

ans =
1 6 7 4

 two dimensional iterators

» a(1:2,1:3)

ans = 1 2 3

5 6 7 .

Working with Memory in Matlab
(index)

To keep track of how memory is allocated, we use 
whos

» b=[34 54 56];

» a=[2 3 56];

» whos

Name Size Bytes Class

a 1x3 24 double array

b 1x3 24 double array

Grand total is 6 elements using 48 bytes 
We may then get rid of unwanted variables using 
clear

» clear a; % removes the variable a from the memory.

We can use clear all to remove all variables from 
memory (and all functions and mex links) and close 
all to close all the figures and hence save a lot of 
memory.

Since Matlab supports dynamic memory allocation, it 
is possible to make inefficient use of the memory due 
to memory fragmentation problems. In this case, the 
pack command can be used to defragment the 
memory and help reclaim some of the memory.

To use memory efficiently, it is best to avoid changing
the number of elements in an array during execution. 
It is best to

preallocate arrays before using them. To this end we 
can use zeros and ones:

» a=zeros(3,4); % a two­dimensional 3x4 zero matrix.

» b=ones(1,5); % an array of 5 elements of value 1. 

Evaluating one­dimensional
 functions in Matlab
(index)

To evaluate one­dimensional functions in Matlab we 
need to specify the points at which we would like our 
function to be evaluated at.

We can use either the built­in colon notation to specify
the points » t=­1:0.5:1

t = ­1.0000 ­0.5000 0 0.5000 1.0000

or use linspace( ) 

» t=linspace(­2,4,5) % generates 5 points between ­2 
and 4

t = ­2.0000 ­0.5000 1.0000 2.5000 4.0000

or use logspace( )

» t=logspace(­3, 2, 6) % 6 points between 10^(­3), 
10^2.

t = 0.0010 0.0100 0.1000 1.0000 10.0000 100.0000

After specifying the points in t, we can apply any 
mathematical function to t. However, we must 
 remember that t represents a vector and all operations 
apply to each element in t » t=[2 3 5 2];

» t+2

ans = 4 5 7 4

When applying addition, substraction, multiplication and 
division between a scalar (that is a single number) and a 
vector we use +, ­, *, and / respectively. 

Evaluating one­dimensional

functions in Matlab (contd)
(index)

To apply element­by­element multiplication or 
division between vectors, we must use the period . in 
front of the operator .* and ./. as in: 

» t1=[1 2 3 4 5];

» t2=[4 5 6 4 5];

» t1.*t2
ans = 4 10 18 16 25

To apply element­by­element addition or substraction
between vectors, we can simply use + and ­.

In general, all numbers in Matlab are double. Complex
numbers are represented by multiplying the imaginary 
part by 1i. For example, 2+3j can be represented by 
1+1i*3 or 1+3i.

All the well­known mathematical functions are 
defined in Matlab, and but they apply to to each 
element in the vector directly. For example, cos([2 3 
4]) is equivalent to [cos(2) cos(3) cos(4)].

Evaluating two­dimensional

functions in Matlab
(index)

To evaluate two­dimensional functions in Matlab we 
use the built­in function meshgrid ( ) with the ranges 
for x and y; eg:

» [x,y]=meshgrid(0.1:0.1:0.3, 0.2:0.1:0.4)
yields the matrices x and y given by:

x =

0.1 0.2 0.3

0.1 0.2 0.3

0.1 0.2 0.3

y =

0.2 0.2 0.2

0.3 0.3 0.3

0.4 0.4 0.4

We may then apply any binary operation to x and y 
provided that we preceede each operation with a 
period. For example x.*y represents the matrix where 
each element of x multiplies each element of y. 
Similarly, exp(x.^2+y.^2) represents the matrix 
resulting from applying exp( ) to the sum of the 
squares of each of the elements of x and y. For our 
example, we have
 Evaluating two­dimensional 

functions in Matlab (contd)
(index)
Even if we need to use conditional statements (if-
statements), we may still be able to avoid for loops.
For example, the elements of a matrix can be clipped
between 0 and 15 using:
» a=[-3 3 4; 5 6 17];
» a=a.*(a<=15).*(a>0)+15*(a>15)
a=
034
5 6 15
To understand how the code works, we note that
conditional statements return 0 when false and 1 when true.
Hence:
» (a<=15).*(a>0)
ans =
011
110
returns 0 for each element of a that is non-positive, and
below 15. The complement of this condition is represented
by all members of a that are above 15:
» a>15
ans =
000
001
Each conditional statement returns 1 for two distinct parts
of the matrix. For each distinct part of the matrix the two
conditionals perform two distinct operations (setting
negative elements to zero and clipping values larger than
15 to 15).

A Fourier Series Example
(index)

In one­dimension, the Fourier Series expansion can be
written in the form of:

where  represents the period:  .
 The following Matlab code demonstrates how to plot 
the sum of the first two harmonics of a function (no 
DC):

» T0 = 4; % This is the period.

» t = linspace (­2*T0, 2*T0, 100); % 100 points in 4 
periods.

» h1 = 3*cos(2*pi/T0*t) + sin(2*pi/T0*t); % first 
harmonic.

» h2 = 6*cos(4*pi/T0*t) + 7*sin(4*pi/T0*t); % second
harmonic.

» S2 = h1+h2; % sum of the first two harmonics.

» plot(t,S2); % plots the sum of the harmonics against 
time.

» set(gca,'FontSize', 18); % set axis fonts large for 
presentations!

which produces the following figure:
 A Review of the Discrete Fourier
Series 

for Real Functions 
(index)

For a real function f, recall the Fourier Series 
expansion:

where  represents the period:  .

By defining  by
 ,

(where  means complex conjugate of the number  ), we 

where the  can be evaluated using the integrals:

Unfortunately, it is not always easy to evaluate the 
continuous­time integrals for  .

We would like to consider discrete­approximations to these
discrete­time integrals.

Periodicity due to Sampling the
Period 
(index)
A very simple approximation to the integral would be to 
sample time at N points within the period:

substitute  for t, and then multiply by the time­distance 
between the points (the Riemann approximation to the 
integral):

Notice that despite the fact that we are only using only N 
points over a single period, after computing expressions for
the cn coefficients, we will have an approximate expression 
of the function for all time:

From the periodicity  , and the fact that
, we get that:

but also:
Hence, our sampling over a single period can be extended 
from:

, 

to:

Implications On Symmetry
(index)

It is important to note that the periodicity that is inherent in 
the sampling:

also imposes sampling constrains in order to preserve the 
symmetry. 
The piecewise­constant approximation must have the 
same symmetries (constant within
):

 for even symmetry: 

Therefore, the first half of the array will also determine the 
second

half of the array through:  . For 
example, 

the array [1 3 4 3] has even symmetry.

 for odd symmetry: 

Again, the first half of the array will also determine the 
second

half of the array through:  . For 
example,

the array [0 3 4 0 ­4 ­3] has odd symmetry (note that first 
and

middle terms are zero, why?).
  for half­wave symmetry: 

Again, the first half of the array will also determine the 
second

half of the array through:  . For

example, the array [2 3 5 ­2 ­3 ­5] has half­wave symmetry.

 for quarter­wave symmetry: 

This symmetry requires half­wave symmetry, and 
also that the function is even within each half­period. 
An example is the array

[0 3 5 3 0 ­3 ­5 ­3]. 

Sampling the Harmonics

(
index)

Recall the approximation to the Fourier Series coefficients:

Hence, the DC and the "discrete" harmonics that we can 
represent:
Can we go forever? Can we sample an infinite number of 
harmonics?

NO, we cannot. 

First, we will establish that we cannot sample any 
harmonic corresponding to any  .

Second, we will establish that any harmonic for
is actually equivalent to a harmonic 
corresponding to a "negative one":  where  .

Together, these observations are the key to 
understanding how to use the FFT correctly.

Sampling the Harmonics (contd)
(index)
For  , we have: 

Similarly, for  :

Hence, the only harmonics that we can sample are for 
positive n, where  .

Sampling Constrains
(index)

We can also conclude that the fastest possible harmonic is 
the one for n=N/2 by using the sampling theorem. 
From the sampling theorem, we know that in order to 
avoid aliasing, we must take at­least two samples from 
each harmonic.

The first sample (for m=0) is always 1/N, and the fastest 
harmonic will come back to 1/N on the third sample. 
Hence, the second sample is the middle of the period: ­1/N, 
and the samples of the highest frequency harmonic are:

This is the  harmonic since:

giving 1, ­1, 1, ... for m=0, 1, ..., N­1.

The Fast Fourier Transform
(index)

Returning to the exponential form of the Fourier Series, we 
use up to the highest harmonic to approximate the original 
function:

where the cn are approximately given by:

 
 Recall that we have 
shown that:

and this leads to:

by simple inspection of the expression for  .

Furthermore, we note the approximation:

which we have shown
to be true for integer values of t.

The Fast Fourier Transform (contd)
(index)

We have that:
At the original signal samples, our approximation becomes:

(*)

Our approximations are closely related to the FFT.

In Matlab, the FFT of  , is given by  defined by:

which is computed using: F=fft (f); while:

can be computed using f=ifft (F);
Hence,  , where we are
using  . Furthermore, the sum in (*) gives 
back  while the second term is an error term.

Plotting the FFT in Matlab
(index)

Prepare a 
harmonic:

» N=16;

» t=0:N­1;

» s=sin(2*pi/N*(19/9)*t);

» plot(s)
» set (gca, 'FontSize', 18);

and plot its FFT:

» stem(abs(fft(s)))

» set (gca, 'FontSize', 18);

» xlabel('Harmonic number');

» ylabel('FFT magnitude');
 

FFT for signals with symmetry
(index)

Cn=(an­jbn)/2, c­n=conj(cn)

 Even symmetry: 
» fft([1 3 4 3])

ans = 11 ­3 ­1 ­3

All the coefficients are real.

 Odd symmetry: 

» fft([0 3 4 0 ­4 ­3])

ans = 0 ­12.1244i 1.7321i 0 ­1.7321i 12.1244i

All the coefficients are imaginary.

 Half­wave symmetry: 

» fft([2 3 5 ­2 ­3 ­5])

ans = 

Columns 1 through 4 

0 2.0000 ­13.8564i 0 8.0000 + 0.0000i

Columns 5 through 6 

0 2.0000 +13.8564i

Coefficients for even harmonics are zero.

 Quarter­wave symmetry: half­wave + possibly 
even/odd.
 » fft([0 3 5 3 0 ­3 ­5 ­3])

ans = 0 ­18.4853i 0 1.5147i 0 ­1.5147i 0 18.4853i 

Coefficients for even harmonics are zero and all other

coefficients are imaginary.

An Overview of

One Dimensional Digital Filtering
(index)

Let  be the input array. Using the arrays a[] and b[], we 
compute the output array y[] using:

We call the a[] to be the array of Finite Impulse Response 
(FIR) filter coefficients, while b[] represents the array of 
Infinite Impulse Response (IIR) filter coefficients.

Examples

Let p=1, q=0. This is the simplest possible filter. It 
multiplies the input array by a constant and stores it into the
output array:
 

Let p=2, q=0. This is the first real FIR filter of interest. It 
is given by: 

There is a problem: y[0] is defined in terms of x[­1], and we
do not know what that value is. We can specify x[­1] to be 
any value we want (called an initial condition), but we will
take it to be zero. Then, we can apply the filter definition to
compute the output array:

We can also assume x[len]=0 to compute y[len].

An Overview of One­Dimensional DSP
(contd)
(index)
Let p=q=1. This is the simplest example of an IIR filter. It 
is defined by:

Again, we have a problem defining y[0] and y[len] since 
y[­1] and x[len] are not known. We will assume that both 
unknowns are zero and proceed to write down expressions 
for y:

Note that if we assign zero to all values of x beyond x[len­
1], we will be able to compute an infinite number of values 
for the output array y. Thus, we have to decide when to stop
this recursion.

It is customary to keep the number of entries in the output 
array to be the same as the number of entries in the input 
array:

.
 One Dimensional Digital Filtering in
Matlab
(index)

In Matlab, digital filters are implemented using the filter 
command. Matlab implements a transposed form. In our 
notation:

implements:

In this form, there is a minus sign in front of the second 
sum, and a scaling by b[0].

For example, a simple implementation of a digital 
derivative is implemented in:

» y = filter ([1 ­1], [1], [3 4 5 6 8 9 10])

y = 3 1 1 1 2 1 1

where the input array is [3 4 5 6 8 9 10], it is convolved 
with the FIR filter [1 ­1]. Note that Matlab provides diff ( ) 
for computing the difference between adjacent samples: 
» diff ([3 4 5 6 8 9 10])

ans = 1 1 1 2 1 1

The difference from our filter example is that the first 3 is 
missing.

For FIR filters only, Matlab provides the conv( ) function 

as demonstrated in:

» conv([3 4 5 6 8 9 10], [1 ­1])

ans = 3 1 1 1 2 1 1 ­10

Notice that the returned array assumes that the input is
extended by zeros.

One Dimensional Digital Filtering

in Matlab (contd)
(index)

Similarly, we can implement a very simple digital 
integrator using an IIR filter:

» y = filter ([1], [1 ­1], [3 4 5 6 8 9 10])
y = 3 7 12 18 26 35 45 

Notice that we have to put a negative sign in the second 1 
in 

[1 ­1].

Matlab also provides cumsum ( ) for implementing the 
same filter:

» cumsum([3 4 5 6 8 9 10])

ans = 3 7 12 18 26 35 45