A Matlab Tutorial
An Introduction to Matlab
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
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
Convert a two-dimensional to a one-dimensional array by simply using array_name(:)
» a(:)
ans =
1
4
7
2
5
8
3
6
9
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
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
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
For two-dimensional arrays, we can use
» 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.
» I=[1 5 8 10]; % indices apply to one-dimensional array a(:)
» a(I)
ans =
1 6 7 4
» a(1:2,1:3)
ans = 1 2 3
5 6 7 .
Working with Memory in Matlab
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
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.
functions in Matlab (contd)
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
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)
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 =
0 3 4
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 =
0 1 1
1 1 0
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 =
0 0 0
0 0 1
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
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
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 can rewrite the Fourier Series expansion in the exponential form:
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
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 c_{n} 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
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 ):
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.
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?).
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.
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
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)
For , we have:
Similarly, for :
Hence, the only harmonics that we can sample are for positive n, where .
Sampling Constrains
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
Returning to the exponential form of the Fourier Series, we use up to the highest harmonic to approximate the original function:
where the c_{n} 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)
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
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
Cn=(an-jbn)/2, c-n=conj(cn)
» fft([1 3 4 3])
ans = 11 -3 -1 -3
All the coefficients are real.
» fft([0 3 4 0 -4 -3])
ans = 0 -12.1244i 1.7321i 0 -1.7321i 12.1244i
All the coefficients are imaginary.
» 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.
» 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
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)
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
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)
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