MATLAB Signal Processing Toolbox Greg Reese Ph D

MATLAB Signal Processing Toolbox Greg Reese, Ph. D Research Computing Support Group Academic Technology Services Miami University October 2013

MATLAB Signal Processing Toolbox © 2013 Greg Reese. All rights reserved 2

Toolbox • Collection of code devoted to solving problems in one field of research • Can be purchased from MATLAB • Can be purchased from third parties • Can be obtained for free from third parties 3

Toolbox MATLAB Signal Processing Toolbox • Code for solving problems in signal processing (!) • Sold by MATLAB • Part of both Miami’s student and faculty license 4

Overview MATLAB divides Signal Processing Toolbox as follows • Waveforms – Pulses, modulated signals, peak-to-peak and RMS amplitude, rise time/fall time, overshoot/undershoot • Convolution and Correlation – Linear and circular convolution, autocorrelation, autocovariance, cross-correlation, cross-covariance • Transforms – Fourier transform, chirp z-transform, DCT, Hilbert transform, cepstrum, Walsh-Hadamard transform 5

Overview • Analog and Digital Filters – Analog filter design, frequency transformations, FIR and IIR filters, filter analysis, filter structures • Spectral Analysis – Nonparametric and parametric spectral estimation, high resolution spectral estimation • Parametric Modeling and Linear Prediction – Autoregressive (AR) models, linear predictive coding (LPC), Levinson-Durbin recursion • Multirate Signal Processing – Downsampling, upsampling, resampling, anti-aliasing filter, interpolation, decimation 6

Overview Will look very briefly at • Analog and Digital Filters • Spectral Analysis • Parametric Modeling and Linear Prediction • Multirate Signal Processing Will look in more depth at • Waveforms • Convolution and Correlation 7

Analog and Digital Filters Toolbox especially good for those serious about their filter design! Analog filters • Standard filters – Bessel, Butterworth, Chebyshev, Elliptic • Filter transforms – Low pass to: bandpass, bandstop, or highpass – Change cutoff frequency of lowpass • Analog to digital filter conversion – Bilinear transformation 8

Analog and Digital Filters Digital Filter Design with functions • Standard filters – Butterworth, Chebyshev, Elliptic • FIR and IIR design – Low pass to: bandpass, bandstop, or highpass – Change cutoff frequency of lowpass • Objects for specification of filters – Arbitrary, lowpass, highpass, bandpass, Hilbert 9

Analog and Digital Filters Digital Filter Design interactively (GUI) • Filterbuilder – specify desired characteristics first, then choose filter type – Butterworth, Chebyshev, Elliptic • FDATool (Filter Design and Analysis Tool) – Quickly design digital FIR or IIR filters – Quickly analyze filters, e. g. , magnitude/phase response, pole-zero plots 10

Analog and Digital Filters SPTool – composite of four tools 1. Signal Browser – analyze signals 2. FDATool 3. FVTool (Filter Visualization Tool) – analyze filter characteristics 4. Spectrum Viewer – spectral analysis 11

Analog and Digital Filters Digital Filter Analysis • Magnitude and phase response, impulse response, group delay, pole-zero plot Digital Filter Implementation • Filtering, direct form, lattice, biquad, statespace structures 12

Spectral Analysis Nonparametric Methods – Periodogram, Welch's overlapped segment averaging, multitaper, cross-spectrum, coherence, spectrogram Parametric Methods – Yule-Walker, Burg, covariance, modified covariance Subspace Methods – Multiple signal classification (MUSIC), eigenvectorestimator, pseudospectrum • Windows – Hamming, Blackman, Bartlett, Chebyshev, Taylor, Kaiser 13

Parametric Modeling and Linear Prediction Parametric Modeling • AR, ARMA, frequency response modeling Linear Predictive Waveforms • Linear predictive coefficients (LPC), line spectral frequencies (LSF), reflection coefficients (RC), Levinson-Durbin recursion 14

Multirate Signal Processing Multirate signal processing – Downsampling, upsampling, resampling, anti -aliasing filter, interpolation, decimation 15

Waveforms part of toolbox lets you create many commonly used signals, which you can use to study models programmed in MATLAB Uses of waveforms • Testing – E. g. , have simple waveform and can analytically determine model’s output. Use toolbox to create that waveform, run it through MATLAB model, confirm result 16

Waveforms Uses of waveforms • Simulation – Most of time can’t get analytical output – Make waveform of known characteristics and study model’s response • Modeling of real signals – Create waveform that looks like the real signal 17

Waveforms Time vectors Digital signals usually sampled from analog at fixed intervals Δt. Want time axis with N points 0 1Δt 2Δt … (N-2)Δt (N-1)Δt 18

For – – Waveforms t. Start: starting time t. End: ending time N: number of points delta. T: sampling interval • If have starting time, number of points, interval, (t. Start, N, delta. T): >> delta. T = 0. 1; >> N = 10; >> t 0 = 5; >> t = t 0 + delta. T * (0: N-1) t = 5. 0 5. 1 5. 2 5. 3 5. 4 5. 5 5. 6 5. 7 5. 8 5. 9 19

Waveforms • If have starting time, ending time, interval (t. Start, t. End, delta. T) >> t. Start = 5; >> t. End = 5. 9; >> delta. T = 0. 1; >> t = t. Start: delta. T: t. End t = 5. 0 5. 1 5. 2 5. 3 5. 4 5. 5 5. 6 5. 7 5. 8 5. 9 • If have starting time, ending time, number of points (t. Start, t. End, N) >> t. Start = 5; >> t. End = 5. 9; >> N = 10; >> t = linspace( t. Start, t. End, N ) t = 5. 0 5. 1 5. 2 5. 3 5. 4 5. 5 5. 6 5. 7 5. 8 5. 9 20

Waveforms In multichannel processing, a number of signals are gathered at the same time • Will assume all sampled at same time and same rate Signal processing toolbox, and MATLAB in general, treats each column of a matrix (2 D array) as an independent column vector 21

Waveforms Example >> M = [ 1: 3; 4: 6; 7: 9; 10: 12 ] M = 1 2 3 4 5 6 7 8 9 10 11 12 >> mean( M ) ans = 5. 5000 6. 5000 7. 5000 Result is mean of each column 22

Waveforms TIP Make time vector be a column vector • Any vectors created from time vector will also be column vectors and so can be processed more easily >> t = t t = 1 2 3 4 5 6 >> y = abs( t - 3 ) y = 2 1 0 1 Column vector begeteth column vector 2 3 23

TIP Waveforms repmat (replicate matrix) is general purpose function to make large matrix by replicating small one Trick - quick way to replicate column vector, i. e. , to make an m x n matrix T out of an m x 1 column vector v, is T = v(: , ones(1, n)) >> v = (1: 5) v = 1 2 3 4 5 >> v(: , ones(1, 3)) ans = 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 24

Waveforms TIP Trick - quick way to replicate row vector, i. e. , to make an m x n matrix T out of an 1 x n row vector v, is T = v(ones(m, 1), : ) >> v = 1: 3 v = 1 2 3 4 5 >> T = v( ones(6, 1), : ) ans = 1 2 3 1 2 3 25

Waveforms TIP Can use preceding two tips to make multichannel signal, e. g. , Simulate the multichannel signal sin(2πt), sin(2πt/2), sin(2πt/3), sin(2πt/4) sampled for one second at one-tenth second per sample 26

Waveforms TIP >> t = (0: 0. 1: 1)' t = 0 0. 1000 0. 2000 0. 3000 0. 4000 0. 5000 0. 6000 0. 7000 0. 8000 0. 9000 1. 0000 27

TIP Waveforms >> T = t(: , ones(1, 4)) T = 0 0 0. 1000 0. 2000 0. 3000 0. 4000 0. 5000 0. 6000 0. 7000 0. 8000 0. 9000 1. 0000 28

TIP Waveforms >> M = T. / C M = 0 0 0. 1000 0. 0500 0. 2000 0. 1000 0. 3000 0. 1500 0. 4000 0. 2000 0. 5000 0. 2500 0. 6000 0. 3000 0. 7000 0. 3500 0. 8000 0. 4000 0. 9000 0. 4500 1. 0000 0. 5000 0 0. 0333 0. 0667 0. 1000 0. 1333 0. 1667 0. 2000 0. 2333 0. 2667 0. 3000 0. 3333 0 0. 0250 0. 0500 0. 0750 0. 1000 0. 1250 0. 1500 0. 1750 0. 2000 0. 2250 0. 2500 29

TIP Waveforms >> signal = sin( 2*pi*M signal = 0 0 0. 5878 0. 3090 0. 9511 0. 5878 0. 9511 0. 8090 0. 5878 0. 9511 0. 0000 1. 0000 -0. 5878 0. 9511 -0. 9511 0. 8090 -0. 9511 0. 5878 -0. 5878 0. 3090 -0. 0000 ) 0 0. 2079 0. 4067 0. 5878 0. 7431 0. 8660 0. 9511 0. 9945 0. 9511 0. 8660 0 0. 1564 0. 3090 0. 4540 0. 5878 0. 7071 0. 8090 0. 8910 0. 9511 0. 9877 1. 0000 30

Waveforms Impulse • Use to compute impulse response of linear, time-invariant system >> t = (0: 0. 1: 1. 9)'; >> impulse = zeros( size(t) ); >> impulse( 1 ) = 1; 31

Waveforms Step • Use to model switch turning on >> t = (-1: 0. 9)'; >> step = [ zeros(10, 1); ones(10, 1) ]; 32

Waveforms Ramp • Use to model something gradually turning on >> t = (0: 0. 1: 1. 9)'; >> ramp = t; 33

Autocorrelation – the detection of a delayed version of a signal • In temporal signal, delay often called “lag” • In spatial signal, delay often called “translation” or “offset” • Delayed signal may also be scaled Can also think of autocorrelation as similarity of a signal to itself as a function of lag 34

Autocorrelation and cross-correlation • Common in both signal processing and statistics • Definitions and uses are different • When looking for help on these topics, make sure you’re looking at a signal-processing source 35

Autocorrelation Examples of autocorrelation of digital signals • Radar – determine distance to object • Sonar – determine distance to object • Music – Determine tempo – Detect and estimate pitch • Astronomy – Find rotation frequency of pulsars 36

Autocorrelation Examples of spatial autocorrelation • Optical Character Recognition (OCR) – reading text from images of writing/printing • X-ray diffraction – helps recover the "Fourier phase information" on atom positions • Statistics – helps estimate mean value uncertainties when sampling a heterogeneous population • Astrophysics – used to study and characterize the spatial distribution of galaxies 37

Autocorrelation Examples of optical autocorrelation • Measurement of optical spectra and of very short-duration light pulses produced by lasers • Analysis of dynamic light scattering data to determine particle size distributions • The small-angle X-ray scattering intensity of some systems related to the spatial autocorrelation function of the electron density • In optics, normalized autocorrelations and cross-correlations give the degree of coherence of an electromagnetic field 38

Autocorrelation Typical use • Blip sent to object • Small blip reflected from object back to sender • Use autocorrelation to detect small blip at some lag • Know velocity of blip in medium so total distance blip traveled is distance = velocity * lag • Distance is round trip, so object distance/2 away 39

Autocorrelation • Multiply and sum. Result is autocorrelation at that point • Slide over one, multiply and sum 10 10 x x 0 0 0 1 1 1 0 0 10∙ 0 +=10∙ 0 0 =0 40

Autocorrelation • Repeat, sliding in both directions until have covered all positions • What happens when go past end? 10 10 10 0 0 1 1 1 x x x 0 0 ? 10∙ 0 + ? 41

Autocorrelation When go past end, have two options Zero – padding – put zeros on both end of both signals. Can either imagine they are there or actually extend arrays in memory and put in zeros 0 0 0 0 1 10 10 10 x x x 0 0 0 10∙ 0 + 10∙ 0 = 0 Will discuss second option later 42

Autocorrelation • 43

Autocorrelation • 44

Autocorrelation TRY IT At time n=0 a transmitter sends out a pulse of amplitude ten and duration 3. At time n=5 it gets back the reflected pulse with the same duration but one tenth the amplitude. What is the autocorrelation? 45

Autocorrelation TRY IT Sent Received 10 10 10 0 0 0 10 10 10 0 0 1 1 1 0 0 Rxy(0) = ? 0 0 1 1 1 0 0 Rxy(0)=0 1 1 0 0 Rxy(1)=0 1 1 0 0 Rxy(2)=0 Rxy(1) = ? 0 0 0 1 Rxy(2) = ? 10 0 0 0 1 46

Autocorrelation TRY IT Rxy(-1)=0 0 10 10 0 0 1 10 10 10 Rxy(-2) = ? 0 0 0 1 1 10 10 10 Rxy(-3) = ? 0 0 1 1 10 10 10 0 1 10 10 10 1 1 1 Rxy(-2)=0 0 0 Rxy(-3)=10 0 Rxy(-4)=20 0 Rxy(-1) = ? 10 0 Rxy(-5)=30 0 0 1 1 0 0 0 Rxy(-4) = ? 0 0 Rxy(-5) = ? 0 0 47

Autocorrelation TRY IT Rxy(-6)=20 Rxy(-6) = ? 10 10 10 0 0 1 1 0 0 Rxy(-7) = ? 10 10 10 0 10 10 0 0 Rxy(-7)=10 0 0 1 Rxy(-8)=0 0 0 1 1 Rxy(-8) = ? 0 0 0 1 Rxy(-9)=0 0 0 1 1 1 Rxy(-9) = ? 10 1 0 10 10 10 48

Autocorrelation TRY IT Put it together >> m = 2: -1: -9; >> R = [ 0 0 0 10 20 30 20 10 0 0 ]; >> [~, max. Index] = max( R ) max. Index = 8 >> m(max. Index) ans = -5 % max when shifted right by 5 0 0 0 1 1 0 1 2 3 4 5 6 1 7 0 8 0 9 49

Autocorrelation TRY IT >> plot( -m, R, ‘o’ ) Note shape is a triangle, not a rectangle, which is shape of pulse. Autocorrelation detects signal of given shape – it does not replicate signal 50

Autocorrelation MATLAB considers what we’re doing to be cross-correlation • Concept is same as what described here for autocorrelation • If one array shorter than another, MATLAB pads shorter one with zeros until both same length 51

Autocorrelation To compute cross-correlation of vectors x and y in MATLAB, use c = xcorr( x, y ) where • c is a vector with 2 N-1 elements • N is length of longer of x and y If m is lag as previously defined, c(k) is autocorrelation for lag m = k - N 52

Autocorrelation TRY IT Let’s do previous graphical autocorrelation with MATLAB 10 10 10 0 0 1 1 1 0 0 >> x = [ 10 10 10 ]; >> y = [ 0 0 0 1 1 1 0 0 ]; >> c = xcorr( x, y ); >> [~, max. Cix] = max( c ) max. Cix = 5 >> m = max. Cix - length( y ) m = -5 % Move 5 to right from element 1 53

Autocorrelation TRY IT Example of finding signal buried in noise 1. Make a sine wave with a period of 20 and amplitude of 100 >> wave = 100 * sin( 2*pi*(0: 19)/20 ); 2. Reset random number generator (so we all get the same random numbers) >> rng default 3. Make 500 points of noise with randn and variance 75% of wave amplitude noisy. Wave = 75 * randn( 1, 500 ); 54

Autocorrelation TRY IT 4. Pick a random spot to place the wave, ensuring that the whole wave fits in ix = randi( [ 1, 481 ] ); 5. Add the wave to the noisy. Wave(ix: ix+19)=noisy. Wave(ix: ix+19)+wave; Plot wave in noise. Is wave visible? plot( noisy. Wave ) 6. Compute autocorrelation >> c = xcorr( wave, noisy. Wave ); 55

Autocorrelation TRY IT 7. Find max of autocorrelation and calculate lag from that >> [~, max. Ix] = max( c ) max. Ix = 269 >> m = max. Ix – 500 m = -231 % shift 231 to right 8. Show random spot where wave added to noise. Match? >> ix ix = 231 Very close match! m should be -230 56

Autocorrelation TRY IT 9. For grins, plot autocorrelation >> lags = -499: 499; >> plot( lags, c ) Why is almost all of right size zero? • Right side corresponds to shifting left and once shift wave more than 20, rest of wave is zeros 57

Correlation Questions? 58

Convolution convolution • Uses – Polynomial multiplication – LTI response – Joint PDF • Linear and circular – Explain, show when equivalent (padding), good for computing convolution with fft. Do example with tic, toc, time-domain convolution vs. fft, ifft, see cconv 59

Convolution Applications of convolution • Acoustics – reverberation is the convolution of the original sound with echoes from objects surrounding the sound source • Computational fluid dynamics – the large eddy simulation (LES) turbulence model uses convolution to lower the range of length scales necessary in computation and thereby reducing computational cost • Probability – probability distribution of the sum of two independent random variables is the convolution of their individual distributions 60

Convolution Applications of convolution • Spectroscopy – line broadening can be due to the Doppler effect and/or collision broadening. The effect due to both is the convolution of the two effects • Electronic music – imposition of a rhythmic structure on a sound done by convolution • Image processing – blurring, sharpening, and edge enhancement done by convolution • Numerical computation – can multiply polynomials quickly with convolution 61

Convolution finds many applications because it is central to linear, timeinvariant systems and many things can be modeled by such systems 62

Convolution Linear – a linear system obeys two principles – Principle of superposition – the output to a sum of inputs is equal to the sum of the outputs to the individual inputs – Scaling – the output to the product of an input and a constant is the product of the constant and the output to the input alone – In other words, for a linear system L, a L{ x(t) } + b. L{ x(t) } = L{ ax(t) + bx(t) } 63

Convolution Suppose you put some input into a system and get some output. If you put in the same input at a later time, if the system is time invariant, the output will be the same as the original output except it will occur at that later time – In other words, for a time-invariant system S, If y(t) = S{ x(t) }, then y(t-t 0) = S{ x(t-t 0) } Time-invariance and linearity are independent. A linear system can be time-invariant or not. A timeinvariant system can be linear or not. 64

Convolution Example Change machine at a laundromat. Put in dollar bills, press button, get out quarters Linear? 1. Put $1 in, press button, get 4 quarters out 2. Put $2 in, press button, get 8 quarters out 3. (output from $1) + (output from $2) = 12 quarters 4. Put $3, press button, get 12 quarters out 5. Two outputs equal, so system linear 65

Convolution Example Time invariant? 1. Put $1 in, press button, get 4 quarters out 2. Put $2 in, press button, get 8 quarters out An hour later 1. Put $1 in, press button, get 4 quarters out 2. Put $2 in, press button, get 8 quarters out Outputs identical except for same delay as input, so system is time invariant 66

Convolution • 67

Convolution Major fact The output of a linear, time-invariant (LTI) system to any input is the convolution of that input with the system’s impulse response Other words: • The impulse response of an LTI system completely characterizes that system • The impulse response of an LTI system specifies that system 68

Convolution Graphical view of convolving two signals • Pick one signal 3 2 1 • Flip it 180° around left edge 3 2 1 1 2 3 • Position right element of flipped signal over left element of unflipped signal 1 2 3 1 0 2 2 5 -2 -2 0 0 1 69

Convolution Graphical view • Multiply corresponding elements and sum 1 2 3 x 1 0 2 2 5 -2 -2 0 0 1 3∙ 1 3 70

Convolution Graphical view • Slide right, multiply, sum 1 2 3 x x 1 0 2 2 5 -2 -2 0 0 1 2∙ 1 + 3∙ 0 3 2 71

Convolution Graphical view • Repeat until “fall off” right side 1 2 3 x x x 1 0 2 2 5 -2 -2 0 0 1 1∙ 1 + 2∙ 0 + 3∙ 2 3 2 7 72

Convolution Graphical view 1 1 2 3 x x x 0 2 2 5 -2 -2 0 0 1 1∙ 0 + 2∙ 2 + 3∙ 2 3 2 7 10 73

Convolution Graphical view 1 0 1 2 3 x x x 2 2 5 -2 -2 0 0 1 1∙ 2 + 2∙ 2 + 3∙ 5 3 2 7 10 21 74

Convolution Graphical view 1 0 2 1 2 3 x x x 2 5 -2 -2 0 0 1 1∙ 2 + 2∙ 5 + 3∙(-2) 3 2 7 10 21 6 75

Convolution Graphical view 1 0 2 2 1 2 3 x x x 5 -2 -2 0 0 1 1∙ 5 + 2∙(-2) + 3∙(-2) 3 2 7 10 21 6 -5 76

Convolution Graphical view 1 0 2 2 5 1 2 3 x x x -2 -2 0 0 1 1∙(-2) + 2∙(-2) + 3∙ 0 3 2 7 10 21 6 -5 -6 77

Convolution Graphical view 1 0 2 2 5 -2 1 2 3 x x x -2 0 0 1 1∙(-2) + 2∙ 0 + 3∙ 0 3 2 7 10 21 6 -5 -6 -2 78

Convolution Graphical view 1 0 2 2 5 -2 -2 1 2 3 x x x 0 0 1 1∙ 0 + 2∙ 0 + 3∙ 1 3 2 7 10 21 6 -5 -6 -2 3 79

Convolution Graphical view • As with correlation, ignore elements that have fallen off or pad bottom array with zeros 1 0 2 2 5 -2 -2 0 1 2 x x 0 1 3 1∙ 0 + 2∙ 1 3 2 7 10 21 6 -5 -6 -2 3 2 80

Convolution Graphical view 1 2 3 2 1 x 1 0 2 2 5 -2 -2 0 0 1 1∙ 1 3 2 7 10 21 6 -5 -6 -2 3 81

Convolution So convolution of 3 2 1 with 1 0 2 2 5 -2 -2 0 0 1 gives 3 2 7 10 21 6 -5 -6 -2 3 2 1 82

Convolution 3 2 1 Note that when we convolved with 1 0 2 2 5 -2 -2 0 0 1 the shorter signal was not completely on top of the longer one for the first two and last two elements of the longer. For these cases, the multiply-and-add convolution computation was missing either one or two terms, so those four calculations are not valid and should be ignored. Bad values at the left and right of a convolution are known as edge effects. 83

Convolution In general, if the shorter of two signals in a convolution has M elements, you should ignore the first and last M-1 elements in the result 84

Convolution • 85

Convolution Compute a convolution in MATLAB with w = conv( u, v ) where u and v are vectors. The output vector w has length(u)+length(v)-1 86

Convolution TRY IT We figured out graphically that convolved with 1 0 2 2 5 gave 3 2 7 10 21 6 -5 -6 3 2 1 -2 -2 0 0 -2 3 2 1 1 Try it in MATLAB >> u = >> v = >> w = 3 [ 3 2 [ 1 0 conv( 2 7 1 ]; 2 2 5 -2 -2 0 0 1 ]; u, v ) 10 21 6 -5 -6 -2 87

Convolution Typically convolution involves a data signal and another signal. The second signal, called a kernel, is • Also called a filter • Usually much shorter than the data signal • Designed to produce a desired effect on the data 88

Convolution TRY IT To introduce a time lag of m units into the data signal, i. e. , to shift it to the right by m, use a kernel of m zeros followed by a one. Example 1 0 2 2 5 -2 -2 0 0 1 Introduce a lag of two into the signal >> kernel = [ 0 0 1 ]; >> w = conv( kernel, v ) w = 0 0 1 0 2 2 5 -2 -2 0 0 1 89

Convolution TRY IT Approximate the derivative by replacing a point with the difference between itself and the previous point Example Approximate the derivative of 1 0 2 2 5 -2 >> u = [ 1 -1 ]; >> w = conv( u, v ) w = 1 -1 2 0 3 -7 -2 0 0 0 2 1 0 1 -1 90

Convolution TRY IT Reduce noise in a signal by replacing a point with a weighted average of itself and the previous two points Example 1 0 2 2 5 -2 -2 0 0 1 >> kernel = (1/9)*[ 5 3 1 ]; >> w = conv( u, v ); w*9 = 5 3 11 16 33 7 -11 -8 -2 5 3 1 91

Convolution Questions? 92

Signal Processing Toolbox Questions? 93

The End 94