DSP LAB-1

LAB 1. Signals in Matlab

Introduction

This 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 Signals

Time 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 Signals

Analog 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.

Aliasing

Digital signals are often derived by sampling a continuous-time signal with an analog-to digital (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 Generation

Signals

• 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(2_ft + _) + 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) = ej2_ft = cos(2_ft) + jsin(2_ft) = 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 Load

Two 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.)

Try out these examples