IntroductionThis lab will describe how to use Matlab for some basic signal representation and manipulation:
• Creating and importing signals
• Sampling and resampling
• Signal visualization
• Modeling noise
• Modulation
Discrete SignalsTime base: t = [0.0 0.1 0.2 0.3]
Signal data: x = [1.0 3.2 2.0 8.5]
The central data construct in Matlab is the numeric array, an ordered collection of real
or complex numeric data with one or more dimensions. The basic data objects of signal
processing (one-dimensional signals or sequences, multichannel signals, and two-dimensional
signals) are all naturally suited to array representation.
Matlab represents ordinary one-dimensional sampled data signals, or sequences, as vectors.
Vectors are 1-by-n or n-by-1 arrays, where n is the number of samples in the sequence.
One way to introduce a sequence into Matlab is to enter it as a list of elements at the command
prompt. The statement
x = [1 2 3 4 5]
creates a simple five-element real sequence in a row vector. It can be converted to a column
vector by taking the transpose:
x = [1 2 3 4 5]’
Column vectors extend naturally to the multichannel case, where each channel is represented
by a column of an array.
c 2006GM
Another method for creating vector data is to use the colon operator. Consider a 1-second
signal sampled at 1000 Hz. An appropriate time vector would be
t = 0:1e-3:1;
where the colon operator creates a 1001-element row vector representing time from zero to
one second in steps of one millisecond.
You can also use linspace to create vector data:
t = linspace(0,1,1e3);
creates a vector of 1000 linearly spaced points between 0 and 1.
Try:
t1 = [0 .1 .2 .3];
t2 = 0:0.1:0.3;
t3 = linspace(0, 0.3, 4);
T = [t1’ t2’ t3’];
X = sin(T)
Q: What does this code show?
Sampling SignalsAnalog signal sources include electromagnetic, audio, sonar, biomedical and others. Analog
signals must be sampled in order to be processed digitally.
Sampling
x(n) = xa(nTs)
x is a discrete signal sampled from the analog signal xa with a sample period of Ts and a
sample frequency of Fs = 1/Ts.
Try:
Fs = 100;
N = 1000;
stoptime = 9.99;
t1 = (0:N-1)/Fs;
t2 = 0:1/Fs:stoptime;
x1 = sin(2*pi*2*t1);
x2 = sin(2*pi*3*t2);
plot(x1)
figure, plot(x2)
An alternative to creating signals is to use a toolbox function. A variety of toolbox functions
generate waveforms . Each of them requires that you begin with a vector representing a time
base. Some of these functions will be described later in this lab.
AliasingDigital signals are often derived by sampling a continuous-time signal with an analog-todigital
(A/D) converter. If the continuous signal, xa(t), is bandlimited, meaning that it does
not contain any frequencies higher than a maximum frequency fM, the Shannon sampling
theorem says that it can be completely recovered from a set of samples if the sampling
frequency fs is greater than two times the maximum frequency of the signal to be sampled:
Fs > 2fM
This maximum frequency fM is known as the Nyquist frequency. If the sampling frequency is
not greater than two times the Nyquist frequency, the continuous signal cannot be uniquely
recovered and aliasing occurs. (You heard examples of aliased signals in Homework No.1).
fs > 2fM: Original signal and sampled signal have the same frequency.
fs 2fM: Sampled signal is aliased to half the original frequency.
Try:t = 0:0.001:2;
xa = sin(2*pi*5*t);
plot(t,xa)
hold on
fs = 15;
ts = 0:1/fs:2;
xs1 = sin(2*pi*5*ts);
plot(ts,xs1,’ro-’)
fs = 7.5;
ts = 0:1/fs:2;
xs2 = sin(2*pi*5*ts);
plot(ts,xs2,’ro-’)
hold off
Q: What is the frequency of xs2?// (There is aliasing here. We need sampling theory.
However can use the fft function on the signal to determine the frequency).
Signal Visualization
• View signal amplitude vs. time index
• Functions: plot, stem, stairs, strips
• Listen to data: sound
Note: the sound and soundsc commands will not work if your computer hardware isn’t set
up. If that is the case, view the signals instead of listening to them.
Try:t = [0.1 0.2 0.3 0.4];
x = [1.0 8.0 4.5 9.7];
plot(t,x)
figure, stem(t,x)
figure, stairs(t,x)
fs = 1000;
ts = 0:1/fs:2;
f = 250 + 240*sin(2*pi*ts);
x = sin(2*pi*f.*ts);
strips(x,0.25,fs)
sound(x,fs)
plot(ts,x)
plot(ts(1:200),x(1:200))
Q: What does the strips command do? (See ’help strips’.)
Q: What does the .* operator do?
Signal Processing Tool
The Signal Processing Toolbox application, SPTool, provides a rich graphical environment
for signal viewing, filter design, and spectral analysis.
You can use SPTool to analyze signals, design filters, analyze filters, filter signals, and analyze
signal spectra. You can accomplish these tasks using four GUIs that you access from
within SPTool:
• The Signal Browser is for analyzing signals. You can also play portions of signals using
your computer’s audio hardware.
• The Filter Designer is for designing or editing FIR and IIR digital filters. Note that the
FDATool is the preferred GUI to use for filter designs. FDATool is discussed in later labs.
• The Filter Viewer is for analyzing filter characteristics.
• The Spectrum Viewer is for spectral analysis.
Open SPTool by typing sptool at the command prompt.
Try:sptool
Look at the train signal, FIRbp filter, and trainse spectrum. (You see 3 panes
- Signals, Filters, Spectra. Filter Designer is available through File 7! Preferences. You
can play sounds using the LS icon. When viewing spectra, note that many methods of
determining spectra, including the fft, are available.)
Importing a Signal
You can use SPTool to analyze the signals, filters, or spectra that you create at the Matlab
command line.
You can import signals, filters, or spectra from the Matlab workspace into the SPTool
workspace using the Import item under the File menu.
Try:
fs = 1000;
ts = 0:1/fs:0.5;
f = 250 + 240*sin(2*pi*ts);
x = sin(2*pi*f.*ts);
Import these signals (f and x) into the SPTool and use the tool to examine them.
Q: What are the icons to use for horizontal zoom?
Try zooming in using the mouse.
Signal Browser
The Signal Browser tool is an interactive signal exploration environment. It provides a
graphical view of the signal object(s) currently selected in teh Signals list of SPTool.
Using the Signal Browser you can
• View and compare vector/array signals
• Zoom in on a range of signal data to examine it more closely
• Measure a variety of characteristics of signal data
• Play signal data on audio hardware
To open/activate the Signal Browser for the SPTool,
• Click one or more signals (use the Shift key for multiple selections) in the Signals list of
SPTool.
• Click the View button in the Signals list of SPTool.
Changing Sample Rates
To change the sample rate of a signal in SPTool,
1. Click a signal in the Signals list in SPTool.
2. Select the Sampling frequency item in the Edit menu.
3. Enter the desired sampling frequency and cliick OK.
Try changing the sampling rate of the imported signal.
Signal GenerationSignals
• Create a time base vector
t = [0:0.1:2];
• Create a signal as a function of time
x = sin(pi*t/2);
plot(t,x)
Useful Matlab functions
• Nonperiodic functions
ones, zeros
• Periodic functions
sin, cos, square, sawtooth
Nonperiodic Signals
t = linspace(0,1,11)
• Step:
y = ones(1,11);
stem(y)
• Impulse:
y = [1 zeros(1,10)];
stem(y)
• Ramp:
y = 2*t;
plot(y)
Useful Matlab functions
step, impulse, gensig
Try:
Step function:
fs = 10;
ts = [0:1/fs:5 5:1/fs:10];
x = [zeros(1,51) ones(1,51)];
stairs(ts,x)
Impulse function with width w:
fs = 10;
w = 0.1;
ts = [-1:1/fs:-w 0 w:1/fs:1];
x = [zeros(1,10) 1 zeros(1,10)];
plot(ts,x)
Delta function:
ts = 0:0.5:5;
x = [1 zeros(1,length(ts)-1)];
stem(ts,x)
axis([-1 6 0 2])
Sinusoids
Sinusoid parameters
• Amplitude, A
• Frequency, f
• Phase shift,
• Vertical offset, B
The general form of a sine wave is
y = Asin(2ft + ) + B
Example: generate a sine wave given the following specifications:
• A = 5
• f = 2 Hz
• = /8 radians
t = linspace(0,1,1001);
A = 5;
f = 2;
p = pi/8;
sinewave = A*sin(2*pi*f*t + p);
plot(t, sinewave)
Try:
edit sine_wave
sine_wave
edit sinfun
[A T] = sinfun(1,2,3,4)
Square Waves
Square wave generation is like sine wave generation, but you specify a duty cycle, which is
the percentage of the time over one period that the amplitude is high.
Example:
• duty cycle is 50% (the Matlab default)
• frequency is 4 Hz.
t = linspace(0,1,1001);
sqw1 = square(2*pi*4*t);
plot(t,sqw1)
axis([-0.1 1.1 -1.1 1.1])
Example:• duty cycle is 75%
• frequency is 4 Hz.
t = linspace(0,1,1001);
sqw2 = square(2*pi*4*t,75);
plot(t,sqw2)
axis([-0.1 1.1 -1.1 1.1])
Sawtooth Waves
Sawtooth waves are like square waves except that instead of specifying a duty cycle, you
specify the location of the peak of the sawtooth.
Example:
• peak at the end of the period (the Matlab default)
• frequency is 3 Hz.
t = linspace(0,1,1001);
saw1 = sawtooth(2*pi*3*t);
plot(t,saw1)
Example:• peak is halfway through the period
• frequency is 3 Hz.
t = linspace(0,1,1001);
saw2 = sawtooth(2*pi*3*t,1/2);
plot(t,saw2)
Complex Signals
Periodic signals can be represented by complex exponentials:
x(t) = ej2ft = cos(2ft) + jsin(2ft) = cos(
t) + jsin(
t)
If t is measured in seconds, then f will have units of sec−1, and
will have units of radians/
second.
In signal processing, we associate the unit circle with one sampling cycle, so that a sampling
frequency of Fs is associated with 2 radians, and the Nyquist frequency Fs/2 is associated
with radians. Values of
in the upper half-plane, in units of Hz, then correspond to
frequencies within the sampled signal.
In Matlab, type:x = exp(2*pi*j*f*t);
plot(x)
Matlab recognizes either j or i as the square root of -1, unless you have defined variables j
or i with different values.
Useful Matlab functions
real, imag, abs, angle
Try:
edit zsig
zsig(5)
Look at both figures and describe what you see.
Importing Data
An important component of the Matlab environment is the ability to read and write data
from/to external sources. Matlab has extensive capabilities for interfacing directly with data
from external programs and instrumentation.
In this lab, we concentrate on reading and writing data that has already been stored in
external files.
Files come in a variety of standard formats, and Matlab has specialized routines for working
with each of them. To see a list of supported file formats, type:
help fileformats
To see a list of associated I/O functions, type:
help iofun
Matlab provides a graphical user interface, the Import Wizard, to the various I/O functions.
You access the Wizard by choosing File ! Import Data or by typing:
uiimport
The Matlab command importdata is a programmatic version of the Wizard, accepting all of
the default choices without opening the graphical user interface. You can use importdata in
M-files to read in data from any of the supported file formats.
Matlab also has a large selection of low-level file I/O functions, modeled after those in the C
programming language. These allow you to work with unsupported formats by instructing
Matlab to open a file in memory, position itself within the file, read or write specific formatted
data, and then close the file.
Try:
help fileformats
help iofun
jan = textread(’all_temps.txt’,’%*u%u%*[^\n]’,’headerlines’,4);
[data text] = xlsread(’stockdata.xls’);
plot(data(:,2))
legend(text{1,3})
Explain how the colon operator works in the preceding plot command.
I = importdata(’eli.jpg’);
image(I)
which theme.wav
uiimport
Browse for:
theme.wav
soundsc(data,fs)
Save and LoadTwo data I/O functions are especially useful when working with Matlab variables.
• The save command writes workspace variables to a binary Matlab data file (MAT-file)
with a .mat extension. The file is placed in the current directory.
• The load command reads variables from a MAT-file back into the Matlab workspace.
Although quite specialized, save and load can be used for day-to-day management of your
Matlab computations.
Try:
doc save
doc load
t = 0:0.1:10;
x1 = sin(t);
x2 = sin(2*t);
x3 = sin(3*t);
save myvars
clear
load myvars t x3
Note the list of variables in the workspace tab in the upper left of the Matlab window.
Modeling Noise
To model signals in space, in the atmosphere, in sea water, or in any communications channel,
it is necessary to model noise.
Matlab has two functions for generating random numbers, which can be added to signals to
model noise.
Uniform random numbers
A = rand(m,n);
generates an mxn array of random numbers from the uniform distribution on the interval
[0,1]. To generate uniformly distributed random numbers from the interval [a,b], shift and
stretch:
A = a + (b-a)*rand(m,n);
Gaussian random numbers
A = randn(m,n);
generates an mxn array of random numbers from the standard normal distribution with
mean 0 and standard deviation 1. To generate random numbers from a normal distribution
with mean mu and standard deviation sigma, shift and stretch:
A = mu + sigma*rand(m,n);
Random numbers from other distributions
Random numbers from other distributions can be generated using the uniform random number
generator and knowledge of the distribution’s inverse cumulative distribution function.
Random number generators for several dozen common distributions are available in the
Statistics Toolbox.
Adding Noise to a Signal
noisy signal = signal + noise
y1 = x + rand(size(x)) % uniform noise
y2 = x + randn(size(x)) % Gaussian noise
Example:
Add Gaussian noise to middle C.
fs = 1e4;
t = 0:1/fs:5;
sw = sin(2*pi*262.62*t); % middle C
n = 0.1*randnsize(sw);
swn = sw + n:
Try:
edit noisyC
noisyC
strips(swn, .1,1e4)
Zoom in on the strips plot. (Note: you might have to cut and paste from the noisyC script
to generate swn.)
Pseudorandomness
This number:
0.95012928514718
is the first number produced by the Matlab uniform random number generator with its
default settings. Start up Matlab, set format long, type rand, and you get the number.
If all Matlab users, all around the world, all on different computers, keep getting this same
number, is it really “random”? No, it isn’t. Computers are deterministic machines and
should not exhibit random behavior. If your computer doesn’t access some external device,
like a gamma ray counter or a clock, then it must really be computing pseudorandom
numbers.
A working definition of randomness was given in 1951 by Berkeley professor D. H. Lehmer,
a pioneer in computing and, especially, computational number theory:
A random sequence is a vague notion ... in which each term is unpredictable
to the uninitiated and whose digits pass a certain number of tests traditional with
statisticians ...
Random number generators proceed deterministically from their current state. To view the
current state of rand, type:
s = rand(’state’)
This returns a 35-element vector containing the current state.
To change the state of rand:
rand(’state’,s) Sets the state to s.
rand(’state’,0) Resets the generator to its initial state.
rand(’state’, sum(100*clock)) Sets to a new state each time.
Commands for randn are analogous.
Try:
s = rand(’state’)
format long
rand
rand(’state’,sum(100*clock))
s = rand(’state’)
format long
rand
Resampling
The Signal Processing Toolbox provides a number of functions that resample a signal at a
higher or lower rate.
y = downsample(x,n)
decreases the effective sampling rate of x by keeping every nth sample starting with the first
sample. x can be a vector or a matrix. If x is a matrix, each column is considered a separate
sequence.
y = upsample(x,n)
increases the effective sampling rate of x by inserting n−1 zeros between samples. x can be
a vector or a matrix. If x is a matrix, each column is considered a separate sequence. The
upsampled y has x n samples.
y = resample(x,p,q)
resamples the sequence in vector x at p/q times the original sampling rate, using a polyphase
filter implementation. p and q must be positive integers. The length of y is equal to
ceil(length(x) p/q). If x is a matrix, resample works down the columns of x.
y = interp(x,r)
increases the sampling rate of x by a factor of r. The interpolated vector y is r times longer
than the original input x.
y = decimate(x,r)
reduces the sampling rate of x by a factor of r. The decimated vector y is r times shorter
in length than the input vector x. By default, decimate employs an eighth-order lowpass
Chebyshev Type I filter. It filters the input sequence in both the forward and reverse
directions to remove all phase distortion, effectively doubling the filter order.
Try:
load mtlb
sound(mtlb,Fs)
mtlb4 = downsample(mtlb,4)
mtlb8 = downsample(mtlb,8)
sound(mtlb8,fs/8)
What are the sizes of mtlb, mtlb4, and mtlb8?
(If sound doesn’t work, plot the signals.)
t = 0:0.00025:1;
x = sin(2*pi*30*t) + sin(2*pi*60*t);
y = decimate(x,4);
subplot(211), stem(x(1:120))
axis([0 120 -2 2])
title(’Original Signal’)
subplot(212), stem(y(1:30))
title(’Decimated Signal’)
Modulation and Demodulation
Modulation varies the amplitude, phase, or frequency of a carrier signal with reference to a
message signal.
The Matlab modulate function modulates a message signal with a specified modulation
method. The syntax is
y = modulate(x,fc,fs,’method’)
where:
• x is the message signal.
• fc is the carrier frequency.
• fs is the sampling frequency.
• method is a flag for the desired modulation method (see table below).
Method Description
amdsb-sc or am Amplitude modulation, double side-band, suppressed carrier
amdsb-tc Amplitude modulation, double side-band, transmitted carrier
amssb Amplitude modulation, single side-band
fm Frequency modulation
pm Phase modulation
ppm Pulse position modulation
pwm Pulse width modulation
qam Quadrature amplitude modulation
The demod function performs demodulation, that is, it obtains the original message signal
from the modulated signal. The syntax is:
x = demod(y,fs,fs,’method’)
demod uses any of the methods shown for modulate. The signal x is attenuated relative to
y because demodulation uses lowpass filtering.
Exercise: High and Low
1. Create a signal equal to the sum of two sine waves with the following characteristics:
• 3-second duration
• Sampling frequency = 2 kHz
• Sinusoid 1: frequency 50 Hz (low), amplitude 10, phase = 0
• Sinusoid 2: frequency 950 Hz (high), amplitude 1, phase = 0
2. View and listen to the signal using an M-file
3. Import the signal into SPTool and view it. Listen to the signal.
HAND IN:
ANSWERS TO ALL QUESTIONS STARTING WITH THE LETTER Q IN
FRONT OF IT.
(The material in this lab handout was put together by Paul Beliveau and derives principally from
the MathWorks training document “MATLAB for Signal Processing”, 2006.)
.
