Digital Signal Processing - The Computer Laboratory

bunkietalentedAI and Robotics

Nov 24, 2013 (3 years and 10 months ago)

120 views

DigitalSignalProcessing
MarkusKuhn
ComputerLaboratory
http://www.cl.cam.ac.uk/teaching/0809/DSP/
Lent2009–PartII
Signals

flowofinformation

measuredquantitythatvarieswithtime(orposition)

electricalsignalreceivedfromatransducer
(microphone,thermometer,accelerometer,antenna,etc.)

electricalsignalthatcontrolsaprocess
Continuous-timesignals:voltage,current,temperature,speed,...
Discrete-timesignals:dailyminimum/maximumtemperature,
lapintervalsinraces,sampledcontinuoussignals,...
Electronics(unlikeoptics)canonlydealeasilywithtime-dependentsignals,thereforespatial
signals,suchasimages,aretypicallyfirstconvertedintoatimesignalwithascanningprocess
(TV,fax,etc.).
2
Signalprocessing
Signalsmayhavetobetransformedinorderto

amplifyorfilteroutembeddedinformation

detectpatterns

preparethesignaltosurviveatransmissionchannel

preventinterferencewithothersignalssharingamedium

undodistortionscontributedbyatransmissionchannel

compensateforsensordeficiencies

findinformationencodedinadifferentdomain
Todoso,wealsoneed

methodstomeasure,characterise,modelandsimulatetrans-
missionchannels

mathematicaltoolsthatsplitcommonchannelsandtransfor-
mationsintoeasilymanipulatedbuildingblocks
3
Analogelectronics
Passivenetworks(resistors,capacitors,
inductances,crystals,SAWfilters),
non-linearelements(diodes,...),
(roughly)linearoperationalamplifiers
Advantages:
•passivenetworksarehighlylinear
overaverylargedynamicrange
andlargebandwidths
•analogsignal-processingcircuits
requirelittleornopower
•analogcircuitscauselittleaddi-
tionalinterference
R
U
in
U
out
C L
0
ω(=2πf)
Uout
1/

LC
U
in

U
in
U
out
t
U
in
−U
out
R
=
1
L
Z
t
−∞
U
out
dτ+C
dU
out
dt
4
Digitalsignalprocessing
Analog/digitalanddigital/analogconverter,CPU,DSP,ASIC,FPGA.
Advantages:

noiseiseasytocontrolafterinitialquantization

highlylinear(withinlimiteddynamicrange)

complexalgorithmsfitintoasinglechip

flexibility,parameterscaneasilybevariedinsoftware

digitalprocessingisinsensitivetocomponenttolerances,aging,
environmentalconditions,electromagneticinterference
But:

discrete-timeprocessingartifacts(aliasing)

canrequiresignificantlymorepower(battery,cooling)

digitalclockandswitchingcauseinterference
5
TypicalDSPapplications

communicationsystems
modulation/demodulation,channel
equalization,echocancellation

consumerelectronics
perceptualcodingofaudioandvideo
onDVDs,speechsynthesis,speech
recognition

music
syntheticinstruments,audioeffects,
noisereduction

medicaldiagnostics
magnetic-resonanceandultrasonic
imaging,computertomography,
ECG,EEG,MEG,AED,audiology

geophysics
seismology,oilexploration

astronomy
VLBI,speckleinterferometry

experimentalphysics
sensor-dataevaluation

aviation
radar,radionavigation

security
steganography,digitalwatermarking,
biometricidentification,surveillance
systems,signalsintelligence,elec-
tronicwarfare

engineering
controlsystems,featureextraction
forpatternrecognition
6
Syllabus
Signalsandsystems.Discretesequencesandsystems,theirtypesandproper-
ties.Lineartime-invariantsystems,convolution.Harmonicphasorsaretheeigen
functionsoflineartime-invariantsystems.Reviewofcomplexarithmetic.Some
examplesfromelectronics,opticsandacoustics.
MATLAB.UseofMATLABonPWFmachinestoperformnumericalexperiments
andvisualisetheresultsinhomeworkexercises.
Fouriertransform.Harmonicphasorsasorthogonalbasefunctions.Formsofthe
Fouriertransform,convolutiontheorem,Dirac’sdeltafunction,impulsecombsin
thetimeandfrequencydomain.
Discretesequencesandspectra.Periodicsamplingofcontinuoussignals,pe-
riodicsignals,aliasing,samplingandreconstructionoflow-passandband-pass
signals,spectralinversion.
DiscreteFouriertransform.ContinuousversusdiscreteFouriertransform,sym-
metry,linearity,reviewoftheFFT,real-valuedFFT.
Spectralestimation.Leakageandscallopingphenomena,windowing,zeropadding.
7
Finiteandinfiniteimpulse-responsefilters.Propertiesoffilters,implementa-
tionforms,window-basedFIRdesign,useoffrequency-inversiontoobtainhigh-
passfilters,useofmodulationtoobtainband-passfilters,FFT-basedconvolution,
polynomialrepresentation,z-transform,zerosandpoles,useofanalogIIRdesign
techniques(Butterworth,ChebyshevI/II,ellipticfilters).
Randomsequencesandnoise.Randomvariables,stationaryprocesses,autocor-
relation,crosscorrelation,deterministiccrosscorrelationsequences,filteredrandom
sequences,whitenoise,exponentialaveraging.
Correlationcoding.Randomvectors,dependenceversuscorrelation,covariance,
decorrelation,matrixdiagonalisation,eigendecomposition,Karhunen-Lo`evetrans-
form,principal/independentcomponentanalysis.Relationtoorthogonaltransform
codingusingfixedbasisvectors,suchasDCT.
Lossyversuslosslesscompression.Whatinformationisdiscardedbyhuman
sensesandcanbeeliminatedbyencoders?Perceptualscales,masking,spatial
resolution,colourcoordinates,somedemonstrationexperiments.
Quantization,imageandaudiocodingstandards.A/-lawcoding,deltacod-
ing,JPEGphotographicstill-imagecompression,motioncompensation,MPEG
videoencoding,MPEGaudioencoding.
Note:Thelastthreelecturesonaudio-visualcodingwerepreviouslypartofthecourse“Informa-
tionTheoryandCoding”.AbriefintroductiontoMATLABwasgivenin“UnixTools”.
8
Objectives
Bytheendofthecourse,youshouldbeableto

applybasicpropertiesoftime-invariantlinearsystems

understandsampling,aliasing,convolution,filtering,thepitfallsof
spectralestimation

explaintheaboveintimeandfrequencydomainrepresentations

usefilter-designsoftware

visualiseanddiscussdigitalfiltersinthez-domain

usetheFFTforconvolution,deconvolution,filtering

implement,applyandevaluatesimpleDSPapplicationsinMATLAB

applytransformsthatreducecorrelationbetweenseveralsignalsources

understandandexplainlimitsinhumanperceptionthatareex-
ploitedbylossycompressiontechniques

provideagoodoverviewoftheprinciplesandcharacteristicsofsev-
eralwidely-usedcompressiontechniquesandstandardsforaudio-
visualsignals
9
Textbooks

R.G.Lyons:Understandingdigitalsignalprocessing.Prentice-
Hall,2004.(£45)

A.V.Oppenheim,R.W.Schafer:Discrete-timesignalprocess-
ing.2nded.,Prentice-Hall,1999.(£47)

J.Stein:Digitalsignalprocessing–acomputerscienceper-
spective.Wiley,2000.(£74)

S.W.Smith:Digitalsignalprocessing–apracticalguidefor
engineersandscientists.Newness,2003.(£40)

K.Steiglitz:Adigitalsignalprocessingprimer–withappli-
cationstodigitalaudioandcomputermusic.Addison-Wesley,
1996.(£40)

SanjitK.Mitra:Digitalsignalprocessing–acomputer-based
approach.McGraw-Hill,2002.(£38)
10
Sequencesandsystems
Adiscretesequence{x
n
}

n=−∞
isasequenceofnumbers
...,x
−2
,x
−1
,x
0
,x
1
,x
2
,...
wherex
n
denotesthen-thnumberinthesequence(n∈Z).Adiscrete
sequencemapsintegernumbersontoreal(orcomplex)numbers.
Wenormallyabbreviate{x
n
}

n=−∞
to{x
n
},orto{x
n
}
n
iftherunningindexisnotobvious.
Thenotationisnotwellstandardized.Someauthorswritex[n]insteadofx
n
,othersx(n).
Whereadiscretesequence{x
n
}samplesacontinuousfunctionx(t)as
x
n
=x(t
s
n)=x(n/f
s
),
wecallt
s
thesamplingperiodandf
s
=1/t
s
thesamplingfrequency.
AdiscretesystemTreceivesasinputasequence{x
n
}andtransforms
itintoanoutputsequence{y
n
}=T{x
n
}:
...,x
2
,x
1
,x
0
,x
−1
,......,y
2
,y
1
,y
0
,y
−1
,...
discrete
systemT
11
Propertiesofsequences
Asequence{x
n
}is
absolutelysummable⇔

X
n=−∞
|x
n
|<∞
squaresummable⇔

X
n=−∞
|x
n
|
2
<∞
periodic⇔∃k>0:∀n∈Z:x
n
=x
n+k
Asquare-summablesequenceisalsocalledanenergysignal,and

X
n=−∞
|x
n
|
2
isitsenergy.ThisterminologyreflectsthatifUisavoltagesuppliedtoaload
resistorR,thenP=UI=U
2
/Risthepowerconsumed,and
R
P(t)dttheenergy.
Soevenwherewedropphysicalunits(e.g.,volts)forsimplicityincalculations,it
isstillcustomarytorefertothesquaredvaluesofasequenceaspowerandtoits
sumorintegralovertimeasenergy.
12
Anon-square-summablesequenceisapowersignalifitsaveragepower
lim
k→∞
1
1+2k
k
X
n=−k
|x
n
|
2
exists.
Specialsequences
Unit-stepsequence:
u
n
=

0,n<0
1,n≥0
Impulsesequence:
δ
n
=

1,n=0
0,n6=0
=u
n
−u
n−1
13
Typesofdiscretesystems
Acausalsystemcannotlookintothefuture:
y
n
=f(x
n
,x
n−1
,x
n−2
,...)
Amemory-lesssystemdependsonlyonthecurrentinputvalue:
y
n
=f(x
n
)
Adelaysystemshiftsasequenceintime:
y
n
=x
n−d
Tisatime-invariantsystemifforanyd
{y
n
}=T{x
n
}⇐⇒{y
n−d
}=T{x
n−d
}.
Tisalinearsystemifforanypairofsequences{x
n
}and{x

n
}
T{ax
n
+bx

n
}=aT{x
n
}+bT{x

n
}.
14
Examples:
Theaccumulatorsystem
y
n
=
n
X
k=−∞
x
k
isacausal,linear,time-invariantsystemwithmemory,asaretheback-
warddifferencesystem
y
n
=x
n
−x
n−1
,
theM-pointmovingaveragesystem
y
n
=
1
M
M−1
X
k=0
x
n−k
=
x
n−M+1
++x
n−1
+x
n
M
andtheexponentialaveragingsystem
y
n
=αx
n
+(1−α)y
n−1


X
k=0
(1−α)
k
x
n−k
.
15
Examplesfortime-invariantnon-linearmemory-lesssystems:
y
n
=x
2
n
,y
n
=log
2
x
n
,y
n
=max{min{⌊256x
n
⌋,255},0}
Examplesforlinearbutnottime-invariantsystems:
y
n
=

x
n
,n≥0
0,n<0
=x
n
u
n
y
n
=x
⌊n/4⌋
y
n
=x
n
ℜ(e
ωjn
)
Examplesforlineartime-invariantnon-causalsystems:
y
n
=
1
2
(x
n−1
+x
n+1
)
y
n
=
9
X
k=−9
x
n+k

sin(πkω)
πkω
[0.5+0.5cos(πk/10)]
16
Constant-coefficientdifferenceequations
Ofparticularpracticalinterestarecausallineartime-invariantsystems
oftheform
y
n
=b
0
x
n

N
X
k=1
a
k
y
n−k
z
−1
z
−1
z
−1
x
n
b
0
y
n−1
y
n−2
y
n−3
−a
1
−a
2
−a
3
y
n
Blockdiagramrepresentation
ofsequenceoperations:
z
−1
x
n
x
n
x
n
x

n
x
n−1
ax
n
a
x
n
+x

n
Delay:
Addition:
Multiplication
byconstant:
Thea
k
andb
m
are
constantcoefficients.
17
or
y
n
=
M
X
m=0
b
m
x
n−m
z
−1
z
−1
z
−1
x
n
y
n
b
0
b
1
b
2
b
3
x
n−1
x
n−2
x
n−3
orthecombinationofboth:
N
X
k=0
a
k
y
n−k
=
M
X
m=0
b
m
x
n−m
z
−1
z
−1
z
−1
z
−1
z
−1
z
−1
b
0
y
n−1
y
n−2
y
n−3
x
n
a
−1
0
b
1
b
2
b
3
x
n−1
x
n−2
x
n−3
−a
1
−a
2
−a
3
y
n
TheMATLABfunctionfilterisanefficientimplementationofthelastvariant.
18
Convolution
Alllineartime-invariant(LTI)systemscanberepresentedintheform
y
n
=

X
k=−∞
a
k
x
n−k
where{a
k
}isasuitablychosensequenceofcoefficients.
Thisoperationoversequencesiscalledconvolutionanddefinedas
{p
n
}∗{q
n
}={r
n
}⇐⇒∀n∈Z:r
n
=

X
k=−∞
p
k
q
n−k
.
If{y
n
}={a
n
}∗{x
n
}isarepresentationofanLTIsystemT,with
{y
n
}=T{x
n
},thenwecallthesequence{a
n
}theimpulseresponse
ofT,because{a
n
}=T{δ
n
}.
19
Convolutionexamples
A
B
C
D
E
F
AB
AC
CA
AE
DE
AF
20
Propertiesofconvolution
Forarbitrarysequences{p
n
},{q
n
},{r
n
}andscalarsa,b:

Convolutionisassociative
({p
n
}∗{q
n
})∗{r
n
}={p
n
}∗({q
n
}∗{r
n
})

Convolutioniscommutative
{p
n
}∗{q
n
}={q
n
}∗{p
n
}

Convolutionislinear
{p
n
}∗{aq
n
+br
n
}=a({p
n
}∗{q
n
})+b({p
n
}∗{r
n
})

Theimpulsesequence(slide13)isneutralunderconvolution
{p
n
}∗{δ
n
}={δ
n
}∗{p
n
}={p
n
}

Sequenceshiftingisequivalenttoconvolvingwithashifted
impulse
{p
n−d
}={p
n
}∗{δ
n−d
}
21
CanallLTIsystemsberepresentedbyconvolution?
Anysequence{x
n
}canbedecomposedintoaweightedsumofshifted
impulsesequences:
{x
n
}=

X
k=−∞
x
k
{δ
n−k
}
Let’sseewhathappensifweapplyalinear
(∗)
time-invariant
(∗∗)
system
Ttosuchadecomposedsequence:
T{x
n
}=T


X
k=−∞
x
k
{δ
n−k
}
!
(∗)
=

X
k=−∞
x
k
T{δ
n−k
}
(∗∗)
=

X
k=−∞
x
k
{δ
n−k
}∗T{δ
n
}=


X
k=−∞
x
k
{δ
n−k
}
!
∗T{δ
n
}
={x
n
}∗T{δ
n
}q.e.d.
⇒TheimpulseresponseT{δ
n
}fullycharacterizesanLTIsystem.
22
Exercise1Whattypeofdiscretesystem(linear/non-linear,time-invariant/
non-time-invariant,causal/non-causal,causal,memory-less,etc.)is:
(a)y
n
=|x
n
|
(b)y
n
=−x
n−1
+2x
n
−x
n+1
(c)y
n
=
8
Y
i=0
x
n−i
(d)y
n
=
1
2
(x
2n
+x
2n+1
)
(e)y
n
=
3x
n−1
+x
n−2
x
n−3
(f)y
n
=x
n
e
n/14
(g)y
n
=x
n
u
n
(h)y
n
=

X
i=−∞
x
i

i−n+2
Exercise2
Provethatconvolutionis(a)commutativeand(b)associative.
23
Exercise3Afinite-lengthsequenceisnon-zeroonlyatafinitenumberof
positions.Ifmandnarethefirstandlastnon-zeropositions,respectively,
thenwecalln−m+1thelengthofthatsequence.Whatmaximumlength
cantheresultofconvolvingtwosequencesoflengthkandlhave?
Exercise4Thelength-3sequencea
0
=−3,a
1
=2,a
2
=1isconvolved
withasecondsequence{b
n
}oflength5.
(a)Writedownthislinearoperationasamatrixmultiplicationinvolvinga
matrixA,avector
~
b∈R
5
,andaresultvector~c.
(b)UseMATLABtomultiplyyourmatrixbythevector
~
b=(1,0,0,2,2)
andcomparetheresultwiththatofusingtheconvfunction.
(c)UsetheMATLABfacilitiesforsolvingsystemsoflinearequationsto
undotheaboveconvolutionstep.
Exercise5(a)Findapairofsequences{a
n
}and{b
n
},whereeachone
containsatleastthreedifferentvaluesandwheretheconvolution{a
n
}∗{b
n
}
resultsinanall-zerosequence.
(b)DoeseveryLTIsystemThaveaninverseLTIsystemT
−1
suchthat
{x
n
}=T
−1
T{x
n
}forallsequences{x
n
}?Why?
24
DirectformIandIIimplementations
z
−1
z
−1
z
−1
z
−1
z
−1
z
−1
b
0
b
1
b
2
b
3
a
−1
0
−a
1
−a
2
−a
3
x
n−1
x
n−2
x
n−3
x
n
y
n−3
y
n−2
y
n−1
y
n
=
z
−1
z
−1
z
−1
a
−1
0
−a
1
−a
2
−a
3
x
n
b
3
b
0
b
1
b
2
y
n
Theblockdiagramrepresentationoftheconstant-coefficientdifference
equationonslide18iscalledthedirectformIimplementation.
Thenumberofdelayelementscanbehalvedbyusingthecommuta-
tivityofconvolutiontoswapthetwofeedbackloops,leadingtothe
directformIIimplementationofthesameLTIsystem.
Thesetwoformsareonlyequivalentwithidealarithmetic(noroundingerrorsandrangelimits).
25
Convolution:opticsexample
Ifaprojectivelensisoutoffocus,theblurredimageisequaltothe
originalimageconvolvedwiththeapertureshape(e.g.,afilledcircle):

=
Point-spreadfunctionh(disk,r=
as
2f
):
h(x,y)=

1
r
2
π
,x
2
+y
2
≤r
2
0,x
2
+y
2
>r
2
OriginalimageI,blurredimageB=I∗h,i.e.
B(x,y)=
ZZ
I(x−x

,y−y

)h(x

,y

)dx

dy

a
f
imageplane
s
focalplane
26
Convolution:electronicsexample
R
U
in
CU
out

U
in
U
out
t
0
0
ω(=2πf)
Uout
1/RC
U
in
U
in

2
Anypassivenetwork(R,L,C)convolvesitsinputvoltageU
in
withan
impulseresponsefunctionh,leadingtoU
out
=U
in
∗h,thatis
U
out
(t)=
Z

−∞
U
in
(t−τ)h(τ)dτ
Inthisexample:
U
in
−U
out
R
=C
dU
out
dt
,h(t)=

1
RC
e
−t
RC
,t≥0
0,t<0
27
Whyaresinewavesuseful?
1)Addingtogethersinewavesofequalfrequency,butarbitraryampli-
tudeandphase,resultsinanothersinewaveofthesamefrequency:
A
1
sin(ωt+ϕ
1
)+A
2
sin(ωt+ϕ
2
)=Asin(ωt+ϕ)
with
A=
q
A
2
1
+A
2
2
+2A
1
A
2
cos(ϕ
2
−ϕ
1
)
tanϕ=
A
1
sinϕ
1
+A
2
sinϕ
2
A
1
cosϕ
1
+A
2
cosϕ
2
ωt
A
2
A
A
1
ϕ
2
ϕ
ϕ
1
A
1
sin(ϕ
1
)
A
2
sin(ϕ
2
)
A
2
cos(ϕ
2
)
A
1
cos(ϕ
1
)
Sinewavesofanyphasecanbe
formedfromsinandcosalone:
Asin(ωt+ϕ)=
asin(ωt)+bcos(ωt)
witha=Acos(ϕ),b=Asin(ϕ)andA=

a
2
+b
2
,tanϕ=
b
a
.
28
Note:Convolutionofadiscretesequence{x
n
}withanothersequence
{y
n
}isnothingbutaddingtogetherscaledanddelayedcopiesof{x
n
}.
(Thinkof{y
n
}decomposedintoasumofimpulses.)
If{x
n
}isasampledsinewaveoffrequencyf,sois{x
n
}∗{y
n
}!
=⇒Sine-wavesequencesformafamilyofdiscretesequences
thatisclosedunderconvolutionwitharbitrarysequences.
Thesameappliesforcontinuoussinewavesandconvolution.
2)Sinewavesareorthogonaltoeachother:
Z

−∞
sin(ω
1
t+ϕ
1
)sin(ω
2
t+ϕ
2
)dt“=”0
⇐⇒ω
1
6=ω
2
∨ϕ
1
−ϕ
2
=(2k+1)π/2(k∈Z)
Theycanbeusedtoformanorthogonalfunctionbasisforatransform.
Theterm“orthogonal”isusedhereinthecontextofan(infinitelydimensional)vectorspace,
wherethe“vectors”arefunctionsoftheformf:R→R(orf:R→C)andthescalarproduct
isdefinedasfg=
R

−∞
f(t)g(t)dt.
29
Whyareexponentialfunctionsuseful?
Addingtogethertwoexponentialfunctionswiththesamebasez,but
differentscalefactorandoffset,resultsinanotherexponentialfunction
withthesamebase:
A
1
z
t+ϕ
1
+A
2
z
t+ϕ
2
=A
1
z
t
z
ϕ
1
+A
2
z
t
z
ϕ
2
=(A
1
z
ϕ
1
+A
2
z
ϕ
2
)z
t
=Az
t
Likewise,ifweconvolveasequence{x
n
}ofvalues
...,z
−3
,z
−2
,z
−1
,1,z,z
2
,z
3
,...
x
n
=z
n
withanarbitrarysequence{h
n
},weget{y
n
}={z
n
}∗{h
n
},
y
n
=

X
k=−∞
x
n−k
h
k
=

X
k=−∞
z
n−k
h
k
=z
n


X
k=−∞
z
−k
h
k
=z
n
H(z)
whereH(z)isindependentofn.
Exponentialsequencesareclosedunderconvolutionwith
arbitrarysequences.Thesameappliesinthecontinuouscase.
30
Whyarecomplexnumberssouseful?
1)Theygiveusallnsolutions(“roots”)ofequationsinvolvingpoly-
nomialsuptodegreen(the“

−1=j”story).
2)Theygiveusthe“greatunifyingtheory”thatcombinessineand
exponentialfunctions:
cos(ωt)=
1
2

e
jωt
+e
−jωt

sin(ωt)=
1
2j

e
jωt
−e
−jωt

or
cos(ωt+ϕ)=
1
2

e
jωt+ϕ
+e
−jωt−ϕ

or
cos(ωn+ϕ)=ℜ(e
jωn+ϕ
)=ℜ[(e

)
n
e

]
sin(ωn+ϕ)=ℑ(e
jωn+ϕ
)=ℑ[(e

)
n
e

]
Notation:ℜ(a+jb):=aandℑ(a+jb):=bwherej
2
=−1anda,b∈R.
31
Wecannowrepresentsinewavesasprojectionsofarotatingcomplex
vector.Thisallowsustorepresentsine-wavesequencesasexponential
sequenceswithbasise

.
Aphaseshiftinsuchasequencecorrespondstoarotationofacomplex
vector.
3)Complexmultiplicationallowsustomodifytheamplitudeandphase
ofacomplexrotatingvectorusingasingleoperationandvalue.
Rotationofa2Dvectorin(x,y)-formisnotationallyslightlymessy,
butfortunatelyj
2
=−1doesexactlywhatisrequiredhere:

x
3
y
3

=

x
2
−y
2
y
2
x
2



x
1
y
1

=

x
1
x
2
−y
1
y
2
x
1
y
2
+x
2
y
1

z
1
=x
1
+jy
1
,z
2
=x
2
+jy
2
z
1
z
2
=x
1
x
2
−y
1
y
2
+j(x
1
y
2
+x
2
y
1
)
(x
2
,y
2
)
(x
1
,y
1
)
(x
3
,y
3
)
(−y
2
,x
2
)
32
Complexphasors
Amplitudeandphasearetwodistinctcharacteristicsofasinefunction
thatareinconvenienttokeepseparatenotationally.
Complexfunctions(anddiscretesequences)oftheform
Ae
j(ωt+ϕ)
=A[cos(ωt+ϕ)+jsin(ωt+ϕ)]
(wherej
2
=−1)areabletorepresentbothamplitudeandphasein
onesinglealgebraicobject.
Thankstocomplexmultiplication,wecanalsoincorporateinonesingle
factorbothamultiplicativechangeofamplitudeandanadditivechange
ofphaseofsuchafunction.Thismakesdiscretesequencesoftheform
x
n
=e
jωn
eigensequenceswithrespecttoanLTIsystemT,becauseforeachω,
thereisacomplexnumber(eigenvalue)H(ω)suchthat
T{x
n
}=H(ω){x
n
}
Inthenotationofslide30,wheretheargumentofHisthebase,wewouldwriteH(e

).
33
Recall:Fouriertransform
TheFourierintegraltransformanditsinversearedefinedas
F{g(t)}(ω)=G(ω)=α
Z

−∞
g(t)e
∓jωt
dt
F
−1
{G(ω)}(t)=g(t)=β
Z

−∞
G(ω)e
±jωt

whereαandβareconstantschosensuchthatαβ=1/(2π).
ManyequivalentformsoftheFouriertransformareusedintheliterature.Thereisnostrong
consensusonwhethertheforwardtransformusese
−jωt
andthebackwardstransforme
jωt
,or
viceversa.Wefollowherethoseauthorswhosetα=1andβ=1/(2π),tokeeptheconvolution
theoremfreeofaconstantprefactor;otherspreferα=β=1/

2π,intheinterestofsymmetry.
Thesubstitutionω=2πfleadstoaformwithoutprefactors:
F{h(t)}(f)=H(f)=
Z

−∞
h(t)e
∓2πjft
dt
F
−1
{H(f)}(t)=h(t)=
Z

−∞
H(f)e
±2πjft
df
34
Anothernotationisinthecontinuouscase
F{h(t)}(ω)=H(e

)=
Z

−∞
h(t)e
−jωt
dt
F
−1
{H(e

)}(t)=h(t)=
1

Z

−∞
H(e

)e
jωt

andinthediscrete-sequencecase
F{h
n
}(ω)=H(e

)=

X
n=−∞
h
n
e
−jωn
F
−1
{H(e

)}(t)=h
n
=
1

Z
π
−π
H(e

)e
jωn

whichtreatstheFouriertransformasaspecialcaseofthez-transform
(tobeintroducedshortly).
35
PropertiesoftheFouriertransform
If
x(t)•−◦X(f)andy(t)•−◦Y(f)
arepairsoffunctionsthataremappedontoeachotherbytheFourier
transform,thensoarethefollowingpairs.
Linearity:
ax(t)+by(t)•−◦aX(f)+bY(f)
Timescaling:
x(at)•−◦
1
|a|
X

f
a

Frequencyscaling:
1
|a|
x

t
a

•−◦X(af)
36
Timeshifting:
x(t−Δt)•−◦X(f)e
−2πjfΔt
Frequencyshifting:
x(t)e
2πjΔft
•−◦X(f−Δf)
Parseval’stheorem(totalpower):
Z

−∞
|x(t)|
2
dt=
Z

−∞
|X(f)|
2
df
37
Fouriertransformexample:rectandsinc
TheFouriertransformofthe“rectangularfunction”
rect(t)=





1if|t|<
1
2
1
2
if|t|=
1
2
0otherwise
isthe“(normalized)sincfunction”
F{rect(t)}(f)=
Z
1
2

1
2
e
−2πjft
dt=
sinπf
πf
=sinc(f)
andviceversa
F{sinc(t)}(f)=rect(f).
Somenoteworthypropertiesofthesefunctions:

R

−∞
sinc(t)dt=1=
R

−∞
rect(t)dt
•sinc(0)=1=rect(0)
•∀n∈Z\{0}:sinc(k)=0
38
Convolutiontheorem
Continuousform:
F{(f∗g)(t)}=F{f(t)}F{g(t)}
F{f(t)g(t)}=F{f(t)}∗F{g(t)}
Discreteform:
{x
n
}∗{y
n
}={z
n
}⇐⇒X(e

)Y(e

)=Z(e

)
Convolutioninthetimedomainisequivalentto(complex)scalarmul-
tiplicationinthefrequencydomain.
Convolutioninthefrequencydomaincorrespondstoscalarmultiplica-
tioninthetimedomain.
Proof:z(r)=
R
s
x(s)y(r−s)ds⇐⇒
R
r
z(r)e
−jωr
dr=
R
r
R
s
x(s)y(r−s)e
−jωr
dsdr=
R
s
x(s)
R
r
y(r−s)e
−jωr
drds=
R
s
x(s)e
−jωs
R
r
y(r−s)e
−jω(r−s)
drds
t:=r−s
=
R
s
x(s)e
−jωs
R
t
y(t)e
−jωt
dtds=
R
s
x(s)e
−jωs
ds
R
t
y(t)e
−jωt
dt.(Samefor
P
insteadof
R
.)
39
Dirac’sdeltafunction
Thecontinuousequivalentoftheimpulsesequence{δ
n
}isknownas
Dirac’sdeltafunctionδ(x).Itisageneralizedfunction,definedsuch
that
δ(x)=

0,x6=0
∞,x=0
Z

−∞
δ(x)dx=1
0x
1
andcanbethoughtofasthelimitoffunctionsequencessuchas
δ(x)=lim
n→∞

0,|x|≥1/n
n/2,|x|<1/n
or
δ(x)=lim
n→∞
n

π
e
−n
2
x
2
Thedeltafunctionismathematicallyspeakingnotafunction,butadistribution,thatisan
expressionthatisonlydefinedwhenintegrated.
40
SomepropertiesofDirac’sdeltafunction:
Z

−∞
f(x)δ(x−a)dx=f(a)
Z

−∞
e
±2πjft
df=δ(t)
1

Z

−∞
e
±jωt
dω=δ(t)
Fouriertransform:
F{δ(t)}(ω)=
Z

−∞
δ(t)e
−jωt
dt=e
0
=1
F
−1
{1}(t)=
1

Z

−∞
1e
jωt
dω=δ(t)
http://mathworld.wolfram.com/DeltaFunction.html
41
Sineandcosineinthefrequencydomain
cos(2πf
0
t)=
1
2
e
2πjf
0
t
+
1
2
e
−2πjf
0
t
sin(2πf
0
t)=
1
2j
e
2πjf
0
t

1
2j
e
−2πjf
0
t
F{cos(2πf
0
t)}(f)=
1
2
δ(f−f
0
)+
1
2
δ(f+f
0
)
F{sin(2πf
0
t)}(f)=−
j
2
δ(f−f
0
)+
j
2
δ(f+f
0
)
ℑℑ
ℜℜ
1
2
1
2
1
2
j
1
2
j
f f f
0
−f
0
−f
0
f
0
Asanyreal-valuedsignalx(t)canberepresentedasacombinationofsineandcosinefunctions,
thespectrumofanyreal-valuedsignalwillshowthesymmetryX(e

)=[X(e
−jω
)]

,where

denotesthecomplexconjugate(i.e.,negatedimaginarypart).
42
Fouriertransformsymmetries
Wecallafunctionx(t)
oddifx(−t)=−x(t)
evenifx(−t)=x(t)
and

isthecomplexconjugate,suchthat(a+jb)

=(a−jb).
Then
x(t)isreal⇔X(−f)=[X(f)]

x(t)isimaginary⇔X(−f)=−[X(f)]

x(t)iseven⇔X(f)iseven
x(t)isodd⇔X(f)isodd
x(t)isrealandeven⇔X(f)isrealandeven
x(t)isrealandodd⇔X(f)isimaginaryandodd
x(t)isimaginaryandeven⇔X(f)isimaginaryandeven
x(t)isimaginaryandodd⇔X(f)isrealandodd
43
Example:amplitudemodulation
Communicationchannelsusuallypermitonlytheuseofagivenfre-
quencyinterval,suchas300–3400Hzfortheanalogphonenetworkor
590–598MHzforTVchannel36.Modulationwithacarrierfrequency
f
c
shiftsthespectrumofasignalx(t)intothedesiredband.
Amplitudemodulation(AM):
y(t)=Acos(2πtf
c
)x(t)
00 fff f
l
f
c
−f
l
−f
c
∗=
−f
c
f
c
X(f)Y(f)
Thespectrumofthebasebandsignalintheinterval−f
l
<f<f
l
is
shiftedbythemodulationtotheintervals±f
c
−f
l
<f<±f
c
+f
l
.
Howcansuchasignalbedemodulated?
44
SamplingusingaDiraccomb
Thelossofinformationinthesamplingprocessthatconvertsacon-
tinuousfunctionx(t)intoadiscretesequence{x
n
}definedby
x
n
=x(t
s
n)=x(n/f
s
)
canbemodelledthroughmultiplyingx(t)byacombofDiracimpulses
s(t)=t
s


X
n=−∞
δ(t−t
s
n)
toobtainthesampledfunction
ˆx(t)=x(t)s(t)
Thefunctionˆx(t)nowcontainsexactlythesameinformationasthe
discretesequence{x
n
},butisstillinaformthatcanbeanalysedusing
theFouriertransformoncontinuousfunctions.
45
TheFouriertransformofaDiraccomb
s(t)=t
s


X
n=−∞
δ(t−t
s
n)
=

X
n=−∞
e
2πjnt/t
s
isanotherDiraccomb
S(f)=F
(
t
s


X
n=−∞
δ(t−t
s
n)
)
(f)=
t
s


Z
−∞

X
n=−∞
δ(t−t
s
n)e
2πjft
dt=

X
n=−∞
δ

f−
n
t
s

.
t
s
s(t)S(f)
f
s
−2t
s
−t
s
2t
s
−2f
s
−f
s
2f
s
00f t
46
Samplingandaliasing
0


sample
cos(2 tf)
cos(2 t(k× f
s
± f))
Sampledatfrequencyf
s
,thefunctioncos(2πtf)cannotbedistin-
guishedfromcos[2πt(kf
s
±f)]foranyk∈Z.
47
Frequency-domainviewofsampling
x(t)
ttt
X(f)
fff
00
0
=
............
−1/f
s
1/f
s
1/f
s
0 −1/f
s
s(t)

∗=
−f
s
f
s
0f
s
−f
s
......
S(f)
ˆx(t)
ˆ
X(f)
......
Samplingasignalinthetimedomaincorrespondsinthefrequency
domaintoconvolvingitsspectrumwithaDiraccomb.Theresulting
copiesoftheoriginalsignalspectruminthespectrumofthesampled
signalarecalled“images”.
48
Nyquistlimitandanti-aliasingfilters
Ifthe(double-sided)bandwidthofasignaltobesampledislargerthan
thesamplingfrequencyf
s
,theimagesofthesignalthatemergeduring
samplingmayoverlapwiththeoriginalspectrum.
Suchanoverlapwillhinderreconstructionoftheoriginalcontinuous
signalbyremovingthealiasingfrequencieswithareconstructionfilter.
Therefore,itisadvisabletolimitthebandwidthoftheinputsignalto
thesamplingfrequencyf
s
beforesampling,usingananti-aliasingfilter.
Inthecommoncaseofareal-valuedbase-bandsignal(withfrequency
contentdownto0Hz),allfrequenciesfthatoccurinthesignalwith
non-zeropowershouldbelimitedtotheinterval−f
s
/2<f<f
s
/2.
Theupperlimitf
s
/2forthesingle-sidedbandwidthofabaseband
signalisknownasthe“Nyquistlimit”.
49
Nyquistlimitandanti-aliasingfilters
f f
s
−2f
s
−f
s
02f
s
f f
s
−2f
s
−f
s
02f
s
f −f
s
0 f 0f
s
Withanti-aliasingfilter
X(f)
ˆ
X(f)
X(f)
ˆ
X(f)
Withoutanti-aliasingfilter
double-sidedbandwidth
bandwidth
single-sidedNyquist
limit=f
s
/2
reconstructionfilter
anti-aliasingfilter
Anti-aliasingandreconstructionfiltersbothsuppressfrequenciesoutside|f|<f
s
/2.
50
Reconstructionofacontinuous
band-limitedwaveform
Theidealanti-aliasingfilterforeliminatinganyfrequencycontentabove
f
s
/2beforesamplingwithafrequencyoff
s
hastheFouriertransform
H(f)=
(
1if|f|<
f
s
2
0if|f|>
f
s
2
=rect(t
s
f).
Thisleads,afteraninverseFouriertransform,totheimpulseresponse
h(t)=f
s

sinπtf
s
πtf
s
=
1
t
s
sinc

t
t
s

.
Theoriginalband-limitedsignalcanbereconstructedbyconvolving
thiswiththesampledsignalˆx(t),whicheliminatestheperiodicityof
thefrequencydomainintroducedbythesamplingprocess:
x(t)=h(t)∗ˆx(t)
Notethatsamplingh(t)givestheimpulsefunction:h(t)s(t)=δ(t).
51
Impulseresponseofideallow-passfilterwithcut-offfrequencyf
s
/2:
-3
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
3
0
t× f
s
52
Reconstructionfilterexample
1
2
3
4
5


sampled signal
interpolation result
scaled/shifted sin(x)/x pulses
53
Reconstructionfilters
Themathematicallyidealformofareconstructionfilterforsuppressing
aliasingfrequenciesinterpolatesthesampledsignalx
n
=x(t
s
n)back
intothecontinuouswaveform
x(t)=

X
n=−∞
x
n

sinπ(t−t
s
n)
π(t−t
s
n)
.
Choiceofsamplingfrequency
Duetocausalityandeconomicconstraints,practicalanalogfilterscanonlyapprox-
imatesuchanideallow-passfilter.Insteadofasharptransitionbetweenthe“pass
band”(<f
s
/2)andthe“stopband”(>f
s
/2),theyfeaturea“transitionband”
inwhichtheirsignalattenuationgraduallyincreases.
Thesamplingfrequencyisthereforeusuallychosensomewhathigherthantwice
thehighestfrequencyofinterestinthecontinuoussignal(e.g.,4×).Ontheother
hand,thehigherthesamplingfrequency,thehigherareCPU,powerandmemory
requirements.Therefore,thechoiceofsamplingfrequencyisatradeoffbetween
signalquality,analogfiltercostanddigitalsubsystemexpenses.
54
Exercise6Digital-to-analogconverterscannotoutputDiracpulses.In-
stead,foreachsample,theyholdtheoutputvoltage(approximately)con-
stant,untilthenextsamplearrives.Howcanthisbehaviourbemodeled
mathematicallyasalineartime-invariantsystem,andhowdoesitaffectthe
spectrumoftheoutputsignal?
Exercise7ManyDSPsystemsuse“oversampling”tolessentherequire-
mentsonthedesignofananalogreconstructionfilter.Theyuse(afinite
approximationof)thesinc-interpolationformulatomultiplythesampling
frequencyf
s
oftheinitialsampledsignalbyafactorNbeforepassingitto
thedigital-to-analogconverter.WhilethisrequiresmoreCPUoperations
andafasterD/Aconverter,therequirementsonthesubsequentlyapplied
analogreconstructionfilteraremuchlessstringent.Explainwhy,anddraw
schematicrepresentationsofthesignalspectrumbeforeandafterallthe
relevantsignal-processingsteps.
Exercise8Similarly,explainhowoversamplingcanbeappliedtolessen
therequirementsonthedesignofananaloganti-aliasingfilter.
55
Band-passsignalsampling
Sampledsignalscanalsobereconstructediftheirspectralcomponents
remainentirelywithintheintervalnf
s
/2<|f|<(n+1)f
s
/2for
somen∈N.(Thebasebandcasediscussedsofarisjustn=0.)
Inthiscase,thealiasingcopiesofthepositiveandthenegativefrequencieswillinterleaveinstead
ofoverlap,andcanthereforeberemovedagainwithareconstructionfilterwiththeimpulse
response
h(t)=f
s
sinπtf
s
/2
πtf
s
/2
cos

2πtf
s
2n+1
4
«
=(n+1)f
s
sinπt(n+1)f
s
πt(n+1)f
s
−nf
s
sinπtnf
s
πtnf
s
.
f 0f 0
ˆ
X(f) X(f)anti-aliasingfilterreconstructionfilter

5
4
f
s
f
s
−f
s
−f
s
2
f
s
2
5
4
f
s
n=2
56
Exercise9Reconstructingasampledbasebandsignal:
•GenerateaonesecondlongGaussiannoisesequence{r
n
}(using
MATLABfunctionrandn)withasamplingrateof300Hz.
•Usethefir1(50,45/150)functiontodesignafiniteimpulsere-
sponselow-passfilterwithacut-offfrequencyof45Hz.Usethe
filtfiltfunctioninordertoapplythatfiltertothegeneratednoise
signal,resultinginthefilterednoisesignal{x
n
}.
•Thensample{x
n
}at100Hzbysettingallbuteverythirdsample
valuetozero,resultinginsequence{y
n
}.
•Generateanotherlow-passfilterwithacut-offfrequencyof50Hz
andapplyitto{y
n
},inordertointerpolatethereconstructedfiltered
noisesignal{z
n
}.Multiplytheresultbythree,tocompensatethe
energylostduringsampling.
•Plot{x
n
},{y
n
},and{z
n
}.Finallycompare{x
n
}and{z
n
}.
Whyshouldthefirstfilterhavealowercut-offfrequencythanthesecond?
57
Exercise10Reconstructingasampledband-passsignal:
•Generatea1snoisesequence{r
n
},asinexercise9,butthistime
useasamplingfrequencyof3kHz.
•Applytothataband-passfilterthatattenuatesfrequenciesoutside
theinterval31–44Hz,whichtheMATLABSignalProcessingToolbox
functioncheby2(3,30,[3144]/1500)willdesignforyou.
•Thensampletheresultingsignalat30Hzbysettingallbutevery
100-thsamplevaluetozero.
•Generatewithcheby2(3,20,[3045]/1500)anotherband-pass
filterfortheinterval30–45Hzandapplyittotheabove30-Hz-
sampledsignal,toreconstructtheoriginalsignal.(You’llhaveto
multiplyitby100,tocompensatetheenergylostduringsampling.)
•Plotalltheproducedsequencesandcomparetheoriginalband-pass
signalandthatreconstructedafterbeingsampledat30Hz.
Whydoesthereconstructedwaveformdiffermuchmorefromtheoriginal
ifyoureducethecut-offfrequenciesofbothband-passfiltersby5Hz?
58
Spectrumofaperiodicsignal
Asignalx(t)thatisperiodicwithfrequencyf
p
canbefactoredintoa
singleperiod˙x(t)convolvedwithanimpulsecombp(t).Thiscorre-
spondsinthefrequencydomaintothemultiplicationofthespectrum
ofthesingleperiodwithacombofimpulsesspacedf
p
apart.
=
x(t)
ttt
=∗

X(f)
fff
p(t) ˙x(t)
˙
X(f)P(f)
............
......
−1/f
p
1/f
p
0−1/f
p
1/f
p
0
0f
p
−f
p
0f
p
−f
p
59
Spectrumofasampledsignal
Asignalx(t)thatissampledwithfrequencyf
s
hasaspectrumthatis
periodicwithaperiodoff
s
.
x(t)
ttt
X(f)
fff
00
0
=
............
−1/f
s
1/f
s
1/f
s
0 −1/f
s
s(t)

∗=
−f
s
f
s
0f
s
−f
s
............
S(f)
ˆx(t)
ˆ
X(f)
60
ContinuousvsdiscreteFouriertransform
•Samplingacontinuoussignalmakesitsspectrumperiodic
•Aperiodicsignalhasasampledspectrum
Wesampleasignalx(t)withf
s
,gettingˆx(t).Wetakenconsecutive
samplesofˆx(t)andrepeattheseperiodically,gettinganewsignal¨x(t)
withperiodn/f
s
.Itsspectrum
¨
X(f)issampled(i.e.,hasnon-zero
value)atfrequencyintervalsf
s
/nandrepeatsitselfwithaperiodf
s
.
Nowboth¨x(t)anditsspectrum
¨
X(f)arefinitevectorsoflengthn.
f t
............
f
−1
s
f
−1
s
0 −n/f
s
n/f
s
0f
s
f
s
/n −f
s
/n −f
s
¨x(t)
¨
X(f)
61
DiscreteFourierTransform(DFT)
X
k
=
n−1
X
i=0
x
i
e
−2πj
ik
n
x
k
=
1
n

n−1
X
i=0
X
i
e
2πj
ik
n
Then-pointDFTmultipliesavectorwithann×nmatrix
F
n
=










11111
1e
−2πj
1
n
e
−2πj
2
n
e
−2πj
3
n
e
−2πj
n−1
n
1e
−2πj
2
n
e
−2πj
4
n
e
−2πj
6
n
e
−2πj
2(n−1)
n
1e
−2πj
3
n
e
−2πj
6
n
e
−2πj
9
n
e
−2πj
3(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1e
−2πj
n−1
n
e
−2πj
2(n−1)
n
e
−2πj
3(n−1)
n
e
−2πj
(n−1)(n−1)
n










F
n








x
0
x
1
x
2
.
.
.
x
n−1







=







X
0
X
1
X
2
.
.
.
X
n−1







,
1
n
F

n








X
0
X
1
X
2
.
.
.
X
n−1







=







x
0
x
1
x
2
.
.
.
x
n−1







62
DiscreteFourierTransformvisualized











































x
0
x
1
x
2
x
3
x
4
x
5
x
6
x
7














=














X
0
X
1
X
2
X
3
X
4
X
5
X
6
X
7














Then-pointDFTofasignal{x
i
}sampledatfrequencyf
s
containsin
theelementsX
0
toX
n/2
oftheresultingfrequency-domainvectorthe
frequencycomponents0,f
s
/n,2f
s
/n,3f
s
/n,...,f
s
/2,andcontains
inX
n−1
downtoX
n/2
thecorrespondingnegativefrequencies.Note
thatforareal-valuedinputvector,bothX
0
andX
n/2
willbereal,too.
Whyistherenophaseinformationrecoveredatf
s
/2?
63
InverseDFTvisualized
1
8












































X
0
X
1
X
2
X
3
X
4
X
5
X
6
X
7














=














x
0
x
1
x
2
x
3
x
4
x
5
x
6
x
7














64
FastFourierTransform(FFT)

F
n
{x
i
}
n−1
i=0

k
=
n−1
X
i=0
x
i
e
−2πj
ik
n
=
n
2
−1
X
i=0
x
2i
e
−2πj
ik
n/2
+e
−2πj
k
n
n
2
−1
X
i=0
x
2i+1
e
−2πj
ik
n/2
=








F
n
2
{x
2i
}
n
2
−1
i=0

k
+e
−2πj
k
n


F
n
2
{x
2i+1
}
n
2
−1
i=0

k
,k<
n
2

F
n
2
{x
2i
}
n
2
−1
i=0

k−
n
2
+e
−2πj
k
n


F
n
2
{x
2i+1
}
n
2
−1
i=0

k−
n
2
,k≥
n
2
TheDFTovern-elementvectorscanbereducedtotwoDFTsover
n/2-elementvectorsplusnmultiplicationsandnadditions,leadingto
log
2
nroundsandnlog
2
nadditionsandmultiplicationsoverall,com-
paredton
2
fortheequivalentmatrixmultiplication.
Ahigh-performanceFFTimplementationinCwithmanyprocessor-specificoptimizationsand
supportfornon-power-of-2sizesisavailableathttp://www.fftw.org/.
65
Efficientreal-valuedFFT
ThesymmetrypropertiesoftheFouriertransformappliedtothediscrete
Fouriertransform{X
i
}
n−1
i=0
=F
n
{x
i
}
n−1
i=0
havetheform
∀i:x
i
=ℜ(x
i
)⇐⇒∀i:X
n−i
=X

i
∀i:x
i
=jℑ(x
i
)⇐⇒∀i:X
n−i
=−X

i
Thesetwosymmetries,combinedwiththelinearityoftheDFT,allowsus
tocalculatetworeal-valuedn-pointDFTs
{X

i
}
n−1
i=0
=F
n
{x

i
}
n−1
i=0
{X
′′
i
}
n−1
i=0
=F
n
{x
′′
i
}
n−1
i=0
simultaneouslyinasinglecomplex-valuedn-pointDFT,bycomposingits
inputas
x
i
=x

i
+jx
′′
i
anddecomposingitsoutputas
X

i
=
1
2
(X
i
+X

n−i
)X
′′
i
=
1
2
(X
i
−X

n−i
)
Tooptimizethecalculationofasinglereal-valuedFFT,usethistricktocalculatethetwohalf-size
real-valueFFTsthatoccurinthefirstround.
66
Fastcomplexmultiplication
Calculatingtheproductoftwocomplexnumbersas
(a+jb)(c+jd)=(ac−bd)+j(ad+bc)
involvesfour(real-valued)multiplicationsandtwoadditions.
Thealternativecalculation
(a+jb)(c+jd)=(α−β)+j(α+γ)with
α=a(c+d)
β=d(a+b)
γ=c(b−a)
providesthesameresultwiththreemultiplicationsandfiveadditions.
ThelattermayperformfasteronCPUswheremultiplicationstakethree
ormoretimeslongerthanadditions.
Thistrickismosthelpfulonsimplermicrocontrollers.Specializedsignal-processingCPUs(DSPs)
feature1-clock-cyclemultipliers.High-enddesktopprocessorsusepipelinedmultipliersthatstall
whereoperationsdependoneachother.
67
FFT-basedconvolution
Calculatingtheconvolutionoftwofinitesequences{x
i
}
m−1
i=0
and{y
i
}
n−1
i=0
oflengthsmandnvia
z
i
=
min{m−1,i}
X
j=max{0,i−(n−1)}
x
j
y
i−j
,0≤i<m+n−1
takesmnmultiplications.
CanweapplytheFFTandtheconvolutiontheoremtocalculatethe
convolutionfaster,injustO(mlogm+nlogn)multiplications?
{z
i
}=F
−1
(F{x
i
}F{y
i
})
Thereisobviouslynoproblemifthisconditionisfulfilled:
{x
i
}and{y
i
}areperiodic,withequalperiodlengths
Inthiscase,thefactthattheDFTinterpretsitsinputasasingleperiod
ofaperiodicsignalwilldoexactlywhatisneeded,andtheFFTand
inverseFFTcanbeapplieddirectlyasabove.
68
Inthegeneralcase,measureshavetobetakentopreventawrap-over:
A
B
F
-1
[F(A)×F(B)]
A'
B'
F
-1
[F(A')×F(B')]
Bothsequencesarepaddedwithzerovaluestoalengthofatleastm+n−1.
Thisensuresthatthestartandendoftheresultingsequencedonotoverlap.
69
Zeropaddingisusuallyappliedtoextendbothsequencelengthstothe
nexthigherpoweroftwo(2
⌈log
2
(m+n−1)⌉
),whichfacilitatestheFFT.
Withacausalsequence,simplyappendthepaddingzerosattheend.
Withanon-causalsequence,valueswithanegativeindexnumberare
wrappedaroundtheDFTblockboundariesandappearattheright
end.Inthiscase,zero-paddingisappliedinthecenteroftheblock,
betweenthelastandfirstelementofthesequence.
ThankstotheperiodicnatureoftheDFT,zeropaddingatbothends
hasthesameeffectaspaddingonlyatoneend.
IfbothsequencescanbeloadedentirelyintoRAM,theFFTcanbeap-
pliedtotheminonestep.However,oneofthesequencesmightbetoo
largeforthat.Itcouldalsobearealtimewaveform(e.g.,atelephone
signal)thatcannotbedelayeduntiltheendofthetransmission.
Insuchcases,thesequencehastobesplitintoshorterblocksthatare
separatelyconvolvedandthenaddedtogetherwithasuitableoverlap.
70
Eachblockiszero-paddedatbothendsandthenconvolvedasbefore:
===
∗∗∗
Theregionsoriginallyaddedaszeropaddingare,afterconvolution,aligned
tooverlapwiththeunpaddedendsoftheirrespectiveneighbourblocks.
Theoverlappingpartsoftheblocksarethenaddedtogether.
71
Deconvolution
Asignalu(t)wasdistortedbyconvolutionwithaknownimpulsere-
sponseh(t)(e.g.,throughatransmissionchannelorasensorproblem).
The“smeared”results(t)wasrecorded.
Canweundothedamageandrestore(oratleastestimate)u(t)?

=

Theconvolutiontheoremturnstheproblemintooneofmultiplication:
s(t)=
Z
u(t−τ)h(τ)dτ
s=u∗h
F{s}=F{u}F{h}
F{u}=F{s}/F{h}
u=F
−1
{F{s}/F{h}}
Inpractice,wealsorecordsomenoisen(t)(quantization,etc.):
c(t)=s(t)+n(t)=
Z
u(t−τ)h(τ)dτ+n(t)
Problem–AtfrequenciesfwhereF{h}(f)approacheszero,the
noisewillbeamplified(potentiallyenormously)duringdeconvolution:
˜u=F
−1
{F{c}/F{h}}=u+F
−1
{F{n}/F{h}}
73
Typicalworkarounds:

ModifytheFouriertransformoftheimpulseresponse,suchthat
|F{h}(f)|>ǫforsomeexperimentallychosenthresholdǫ.

Ifestimatesofthesignalspectrum|F{s}(f)|andthenoise
spectrum|F{n}(f)|canbeobtained,thenwecanapplythe
“Wienerfilter”(“optimalfilter”)
W(f)=
|F{s}(f)|
2
|F{s}(f)|
2
+|F{n}(f)|
2
beforedeconvolution:
˜u=F
−1
{WF{c}/F{h}}
Exercise11UseMATLABtodeconvolvetheblurredstarsfromslide26.
Thefilesstars-blurred.pngwiththeblurred-starsimageandstars-psf.pngwiththeimpulse
response(point-spreadfunction)areavailableonthecourse-materialwebpage.Youmayfind
theMATLABfunctionsimread,double,imagesc,circshift,fft2,ifft2ofuse.
Trydifferentwaystocontrolthenoise(seeabove)anddistortionsnearthemargins(window-
ing).[TheMATLABimageprocessingtoolboxprovidesready-made“professional”functions
deconvwnr,deconvreg,deconvlucy,edgetaper,forsuchtasks.Donotusethese,exceptper-
hapstocomparetheiroutputswiththeresultsofyourownattempts.]
74
Spectralestimation
0
10
20
30
-1
0
1
Sine wave 4f
s
/32
0
10
20
30
0
5
10
15
Discrete Fourier Transform
0
10
20
30
-1
0
1
Sine wave 4.61f
s
/32
0
10
20
30
0
5
10
15
Discrete Fourier Transform
75
WeintroducedtheDFTasaspecialcaseofthecontinuousFourier
transform,wheretheinputissampledandperiodic.
Iftheinputissampled,butnotperiodic,theDFTcanstillbeused
tocalculateanapproximationoftheFouriertransformoftheoriginal
continuoussignal.However,therearetwoeffectstoconsider.They
areparticularlyvisiblewhenanalysingpuresinewaves.
Sinewaveswhosefrequencyisamultipleofthebasefrequency(f
s
/n)
oftheDFTareidenticaltotheirperiodicextensionbeyondthesize
oftheDFT.Theyare,therefore,representedexactlybyasinglesharp
peakintheDFT.Alltheirenergyfallsintoonesinglefrequency“bin”
intheDFTresult.
Sinewaveswithotherfrequencies,whichdonotmatchexactlyoneof
theoutputfrequencybinsoftheDFT,arestillrepresentedbyapeak
attheoutputbinthatrepresentsthenearestintegermultipleofthe
DFT’sbasefrequency.However,suchapeakisdistortedintwoways:

Itsamplitudeislower(downto63.7%).

Muchsignalenergyhas“leaked”tootherfrequencies.
76
0
5
10
15
20
25
30
15
15.5
16
0
5
10
15
20
25
30
35
input freq.
DFT index
Theleakageofenergytootherfrequencybinsnotonlyblurstheestimatedspec-
trum.Thepeakamplitudealsochangessignificantlyasthefrequencyofatone
changesfromthatassociatedwithoneoutputbintothenext,aphenomenon
knownasscalloping.Intheabovegraphic,aninputsinewavegraduallychanges
fromthefrequencyofbin15tothatofbin16(onlypositivefrequenciesshown).
77
Windowing
0
200
400
-1
0
1
Sine wave
0
200
400
0
100
200
300
Discrete Fourier Transform
0
200
400
-1
0
1
Sine wave multiplied with window function
0
200
400
0
50
100
Discrete Fourier Transform
78
Thereasonfortheleakageandscallopinglossesiseasytovisualizewiththe
helpoftheconvolutiontheorem:
TheoperationofcuttingasequenceofthesizeoftheDFTinputvectorout
ofalongeroriginalsignal(theonewhosecontinuousFourierspectrumwe
trytoestimate)isequivalenttomultiplyingthissignalwitharectangular
function.Thisdestroysallinformationandcontinuityoutsidethe“window”
thatisfedintotheDFT.
MultiplicationwitharectangularwindowoflengthTinthetimedomainis
equivalenttoconvolutionwithsin(πfT)/(πfT)inthefrequencydomain.
Thesubsequentinterpretationofthiswindowasaperiodicsequenceby
theDFTleadstosamplingofthisconvolutionresult(samplingmeaning
multiplicationwithaDiraccombwhoseimpulsesarespacedf
s
/napart).
Wherethewindowlengthwasanexactmultipleoftheoriginalsignalperiod,
samplingofthesin(πfT)/(πfT)curveleadstoasingleDiracpulse,and
thewindowingcausesnodistortion.Inallothercases,theeffectsofthecon-
volutionbecomevisibleinthefrequencydomainasleakageandscalloping
losses.
79
Somebetterwindowfunctions
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
1


Rectangular window
Triangular window
Hann window
Hamming window
Allthesefunctionsare0outsidetheinterval[0,1].
80
0
0.5
1
-60
-40
-20
0
20
Normalized Frequency ( rad/sample)
Magnitude (dB)
Rectangular window (64-point)
0
0.5
1
-60
-40
-20
0
20
Normalized Frequency ( rad/sample)
Magnitude (dB)
Triangular window
0
0.5
1
-60
-40
-20
0
20
Normalized Frequency ( rad/sample)
Magnitude (dB)
Hann window
0
0.5
1
-60
-40
-20
0
20
Normalized Frequency ( rad/sample)
Magnitude (dB)
Hamming window
81
Numerousalternativestotherectangularwindowhavebeenproposed
thatreduceleakageandscallopinginspectralestimation.Theseare
vectorsmultipliedelement-wisewiththeinputvectorbeforeapplying
theDFTtoit.Theyallforcethesignalamplitudesmoothlydownto
zeroattheedgeofthewindow,therebyavoidingtheintroductionof
sharpjumpsinthesignalwhenitisextendedperiodicallybytheDFT.
Threeexamplesofsuchwindowvectors{w
i
}
n−1
i=0
are:
Triangularwindow(Bartlettwindow):
w
i
=1−




1−
i
n/2




Hannwindow(raised-cosinewindow,Hanningwindow):
w
i
=0.5−0.5×cos


i
n−1

Hammingwindow:
w
i
=0.54−0.46×cos


i
n−1

82
ZeropaddingincreasesDFTresolution
Thetwofiguresbelowshowtwospectraofthe16-elementsequence
s
i
=cos(2π3i/16)+cos(2π4i/16),i∈{0,...,15}.
TheleftplotshowstheDFTofthewindowedsequence
x
i
=s
i
w
i
,i∈{0,...,15}
andtherightplotshowstheDFTofthezero-paddedwindowedsequence
x

i
=

s
i
w
i
,i∈{0,...,15}
0,i∈{16,...,63}
wherew
i
=0.54−0.46×cos(2πi/15)istheHammingwindow.
0
5
10
15
0
2
4
DFT without zero padding
0
20
40
60
0
2
4
DFT with 48 zeros appended to window
83
ApplyingthediscreteFouriertransformtoann-elementlongreal-
valuedsequenceleadstoaspectrumconsistingofonlyn/2+1discrete
frequencies.
Sincetheresultingspectrumhasalreadybeendistortedbymultiplying
the(hypotheticallylonger)signalwithawindowingfunctionthatlimits
itslengthtonnon-zerovaluesandforcesthewaveformsmoothlydown
tozeroatthewindowboundaries,appendingfurtherzerosoutsidethe
windowwillnotdistortthesignalfurther.
ThefrequencyresolutionoftheDFTisthesamplingfrequencydivided
bytheblocksizeoftheDFT.Zeropaddingcanthereforebeusedto
increasethefrequencyresolutionoftheDFT.
Notethatzeropaddingdoesnotaddanyadditionalinformationtothe
signal.Thespectrumhasalreadybeen“low-passfiltered”bybeing
convolvedwiththespectrumofthewindowingfunction.Zeropadding
inthetimedomainmerelysamplesthisspectrumblurredbythewin-
dowingstepatahigherresolution,therebymakingiteasiertovisually
distinguishspectrallinesandtolocatetheirpeakmoreprecisely.
84
Frequencyinversion
InordertoturnthespectrumX(f)ofareal-valuedsignalx
i
sampledatf
s
intoaninvertedspectrumX

(f)=X(f
s
/2−f),wemerelyhavetoshift
theperiodicspectrumbyf
s
/2:
=∗
00 fff
X(f)
−f
s
f
s
0 −f
s
f
s
X

(f)
f
s
2

f
s
2
............
Thiscanbeaccomplishedbymultiplyingthesampledsequencex
i
withy
i
=
cosπf
s
t=cosπi,whichisnothingbutmultiplicationwiththesequence
...,1,−1,1,−1,1,−1,1,−1,...
Soinordertodesignadiscretehigh-passfilterthatattenuatesallfrequencies
foutsidetherangef
c
<|f|<f
s
/2,wemerelyhavetodesignalow-pass
filterthatattenuatesallfrequenciesoutsidetherange−f
c
<f<f
c
,and
thenmultiplyeverysecondvalueofitsimpulseresponsewith−1.
85
Window-baseddesignofFIRfilters
Recallthattheidealcontinuouslow-passfilterwithcut-offfrequency
f
c
hasthefrequencycharacteristic
H(f)=

1if|f|<f
c
0if|f|>f
c
=rect

f
2f
c

andtheimpulseresponse
h(t)=2f
c
sin2πtf
c
2πtf
c
=2f
c
sinc(2f
c
t).
Samplingthisimpulseresponsewiththesamplingfrequencyf
s
ofthe
signaltobeprocessedwillleadtoaperiodicfrequencycharacteristic,
thatmatchestheperiodicspectrumofthesampledsignal.
Therearetwoproblemsthough:

theimpulseresponseisinfinitelylong

thisfilterisnotcausal,thatish(t)6=0fort<0
86
Solutions:

Maketheimpulseresponsefinitebymultiplyingthesampled
h(t)withawindowingfunction

Maketheimpulseresponsecausalbyaddingadelayofhalfthe
windowsize
Theimpulseresponseofann-thorderlow-passfilteristhenchosenas
h
i
=2f
c
/f
s

sin[2π(i−n/2)f
c
/f
s
]
2π(i−n/2)f
c
/f
s
w
i
where{w
i
}isawindowingsequence,suchastheHammingwindow
w
i
=0.54−0.46×cos(2πi/n)
withw
i
=0fori<0andi>n.
Notethatforf
c
=f
s
/4,wehaveh
i
=0forallevenvaluesofi.Therefore,thisspecialcase
requiresonlyhalfthenumberofmultiplicationsduringtheconvolution.Such“half-band”FIR
filtersareused,forexample,asanti-aliasingfilterswhereverasamplingrateneedstobehalved.
87
FIRlow-passfilterdesignexample
-1
0
1
-1
-0.5
0
0.5
1
30
Real Part
Imaginary Part
0
10
20
30
-0.1
0
0.1
0.2
0.3
n (samples)
Amplitude
Impulse Response
0
0.5
1
-60
-40
-20
0
Normalized Frequency ( rad/sample)
Magnitude (dB)
0
0.5
1
-1500
-1000
-500
0
Normalized Frequency ( rad/sample)
Phase (degrees)
order:n=30,cutofffrequency(−6dB):f
c
=0.25×f
s
/2,window:Hamming
88
Wetruncatetheideal,infinitely-longimpulseresponsebymultiplicationwithawindowsequence.
Inthefrequencydomain,thiswillconvolvetherectangularfrequencyresponseoftheideallow-
passfilterwiththefrequencycharacteristicofthewindow.Thewidthofthemainlobedetermines
thewidthofthetransitionband,andthesidelobescauseripplesinthepassbandandstopband.
Convertingalow-passintoaband-passfilter
Toobtainaband-passfilterthatattenuatesallfrequenciesfoutside
therangef
l
<f<f
h
,wefirstdesignalow-passfilterwithacut-off
frequency(f
h
−f
l
)/2andmultiplyitsimpulseresponsewithasine
waveoffrequency(f
h
+f
l
)/2,beforeapplyingtheusualwindowing:
h
i
=(f
h
−f
l
)/f
s

sin[π(i−n/2)(f
h
−f
l
)/f
s
]
π(i−n/2)(f
h
−f
l
)/f
s
sin[π(f
h
+f
l
)]w
i
=∗
00 fff f
h
f
l
H(f)
f
h
+f
l
2
−f
h
−f
l

f
h
−f
l
2
f
h
−f
l
2

f
h
+f
l
2
89
Exercise12ExplainthedifferencebetweentheDFT,FFT,andFFTW.
Exercise13Push-buttontelephonesuseacombinationoftwosinetones
tosignal,whichbuttoniscurrentlybeingpressed:
1209Hz
1336Hz
1477Hz
1633Hz
697Hz
1
2
3
A
770Hz
4
5
6
B
852Hz
7
8
9
C
941Hz
*
0
#
D
(a)Youreceiveadigitaltelephonesignalwithasamplingfrequencyof
8kHz.Youcuta256-samplewindowoutofthissequence,multiplyitwitha
windowingfunctionandapplya256-pointDFT.Whataretheindiceswhere
theresultingvector(X
0
,X
1
,...,X
255
)willshowthehighestamplitudeif
button
9
waspushedatthetimeoftherecording?
(b)UseMATLABtodetermine,whichbuttonsequencewastypedinthe
touchtonesrecordedinthefiletouchtone.wavonthecourse-materialweb
page.
90
Polynomialrepresentationofsequences
Wecanrepresentsequences{x
n
}aspolynomials:
X(v)=

X
n=−∞
x
n
v
n
Exampleofpolynomialmultiplication:
(
1
+
2
v+
3
v
2
)(
2
+
1
v)
2+4v+6v
2
+1v+2v
2
+3v
3
=
2
+
5
v+
8
v
2
+
3
v
3
Comparethiswiththeconvolutionoftwosequences(inMATLAB):
conv([
123
],[
21
])equals[
2583
]
91
Convolutionofsequencesisequivalenttopolynomialmultiplication:
{h
n
}∗{x
n
}={y
n
}⇒y
n
=

X
k=−∞
h
k
x
n−k
↓↓
H(v)X(v)=


X
n=−∞
h
n
v
n
!



X
n=−∞
x
n
v
n
!
=

X
n=−∞

X
k=−∞
h
k
x
n−k
v
n
NotehowtheFouriertransformofasequencecanbeaccessedeasily
fromitspolynomialform:
X(e
−jω
)=

X
n=−∞
x
n
e
−jωn
92
Exampleofpolynomialdivision:
1
1−av
=1+av+a
2
v
2
+a
3
v
3
+=

X
n=0
a
n
v
n
1+av+a
2
v
2
+
1−av
1
1−av
av
av−a
2
v
2
a
2
v
2
a
2
v
2
−a
3
v
3

Rationalfunctions(quotientsoftwopolynomials)canprovideacon-
venientclosed-formrepresentationsforinfinitely-longexponentialse-
quences,inparticulartheimpulseresponsesofIIRfilters.
93
Thez-transform
Thez-transformofasequence{x
n
}isdefinedas:
X(z)=

X
n=−∞
x
n
z
−n
Notethatisdiffersonlyinthesignoftheexponentfromthepolynomialrepresentationdiscussed
onthepreceedingslides.
RecallthattheaboveX(z)isexactlythefactorwithwhichanexpo-
nentialsequence{z
n
}ismultiplied,ifitisconvolvedwith{x
n
}:
{z
n
}∗{x
n
}={y
n
}
⇒y
n
=

X
k=−∞
z
n−k
x
k
=z
n


X
k=−∞
z
−k
x
k
=z
n
X(z)
94
Thez-transformdefinesforeachsequenceacontinuouscomplex-valued
surfaceoverthecomplexplaneC.Forfinitesequences,itsvalueisal-
waysdefinedacrosstheentirecomplexplane.
Forinfinitesequences,itcanbeshownthatthez-transformconverges
onlyfortheregion
lim
n→∞




x
n+1
x
n




<|z|<lim
n→−∞




x
n+1
x
n




Thez-transformidentifiesasequenceunambiguouslyonlyinconjunctionwithagivenregionof
convergence.Inotherwords,thereexistdifferentsequences,thathavethesameexpressionas
theirz-transform,butthatconvergefordifferentamplitudesofz.
Thez-transformisageneralizationoftheFouriertransform,whichit
containsonthecomplexunitcircle(|z|=1):
F{x
n
}(ω)=X(e

)=

X
n=−∞
x
n
e
−jωn
95
Thez-transformoftheimpulse
response{h
n
}ofthecausalLTI
systemdefinedby
k
X
l=0
a
l
y
n−l
=
m
X
l=0
b
l
x
n−l
with{y
n
}={h
n
}∗{x
n
}isthe
rationalfunction
z
−1
z
−1
z
−1
z
−1
z
−1
z
−1
b
0
b
1
a
−1
0
−a
1
x
n−1
x
n
y
n−1
y
n




y
n−k
−a
k
b
m
x
n−m
H(z)=
b
0
+b
1
z
−1
+b
2
z
−2
++b
m
z
−m
a
0
+a
1
z
−1
+a
2
z
−2
++a
k
z
−k
(b
m
6=0,a
k
6=0)whichcanalsobewrittenas
H(z)=
z
k
P
m
l=0
b
l
z
m−l
z
m
P
k
l=0
a
l
z
k−l
.
H(z)hasmzerosandkpolesatnon-zerolocationsinthezplane,
plusk−mzeros(ifk>m)orm−kpoles(ifm>k)atz=0.
96
Thisfunctioncanbeconvertedintotheform
H(z)=
b
0
a
0

m
Y
l=1
(1−c
l
z
−1
)
k
Y
l=1
(1−d
l
z
−1
)
=
b
0
a
0
z
k−m

m
Y
l=1
(z−c
l
)
k
Y
l=1
(z−d
l
)
wherethec
l
arethenon-zeropositionsofzeros(H(c
l
)=0)andthed
l
arethenon-zeropositionsofthepoles(i.e.,z→d
l
⇒|H(z)|→∞)
ofH(z).Exceptforaconstantfactor,H(z)isentirelycharacterized
bythepositionofthesezerosandpoles.
AswiththeFouriertransform,convolutioninthetimedomaincorre-
spondstocomplexmultiplicationinthez-domain:
{x
n
}•−◦X(z),{y
n
}•−◦Y(z)⇒{x
n
}∗{y
n
}•−◦X(z)Y(z)
Delayingasequencebyonecorrespondsinthez-domaintomultipli-
cationwithz
−1
:
{x
n−Δn
}•−◦X(z)z
−Δn
97
-1
-0.5
0
0.5
1
-1
-0.5
0
0.5
1
0
0.25
0.5
0.75
1
1.25
1.5
1.75
2
real
imaginary
|H(z)|
Thisexampleisanamplitudeplotof
H(z)=
0.8
1−0.2z
−1
=
0.8z
z−0.2
whichfeaturesazeroat0andapoleat0.2.
z
−1
y
n
x
n
y
n−1
0.8
0.2
98
H(z)=
z
z−0.7
=
1
1−0.7z
−1
-1
0
1
-1
0
1
Real Part
Imaginary Part
z Plane
0
10
20
30
0
0.5
1
n (samples)
Amplitude
Impulse Response
H(z)=
z
z−0.9
=
1
1−0.9z
−1
-1
0
1
-1
0
1
Real Part
Imaginary Part
z Plane
0
10
20
30
0
0.5
1
n (samples)
Amplitude
Impulse Response
99
H(z)=
z
z−1
=
1
1−z
−1
-1
0
1
-1
0
1
Real Part
Imaginary Part
z Plane
0
10
20
30
0
0.5
1
n (samples)
Amplitude
Impulse Response
H(z)=
z
z−1.1
=
1
1−1.1z
−1
-1
0
1
-1
0
1
Real Part
Imaginary Part
z Plane
0
10
20
30
0
10
20
n (samples)
Amplitude
Impulse Response
100
H(z)=
z
2
(z−0.9e
jπ/6
)(z−0.9e
−jπ/6
)
=
1
1−1.8cos(π/6)z
−1
+0.9
2
z
−2
-1
0
1
-1
0
1
2
Real Part
Imaginary Part
z Plane
0
10
20
30
-2
0
2
n (samples)
Amplitude
Impulse Response
H(z)=
z
2
(z−e
jπ/6
)(z−e
−jπ/6
)
=
1
1−2cos(π/6)z
−1
+z
−2
-1
0
1
-1
0
1
2
Real Part
Imaginary Part
z Plane
0
10
20
30
-5
0
5
n (samples)
Amplitude
Impulse Response
101
H(z)=
z
2
(z−0.9e
jπ/2
)(z−0.9e
−jπ/2
)
=
1
1−1.8cos(π/2)z
−1
+0.9
2
z
−2
=
1
1+0.9
2
z
−2
-1
0
1
-1
0
1
2
Real Part
Imaginary Part
z Plane
0
10
20
30
-1
0
1
n (samples)
Amplitude
Impulse Response
H(z)=
z
z+1
=
1
1+z
−1
-1
0
1
-1
0
1
Real Part
Imaginary Part
z Plane
0
10
20
30
-1
0
1
n (samples)
Amplitude
Impulse Response
102
IIRFilterdesigntechniques
Thedesignofafilterstartswithspecifyingthedesiredparameters:

Thepassbandisthefrequencyrangewherewewanttoapprox-
imateagainofone.

Thestopbandisthefrequencyrangewherewewanttoapprox-
imateagainofzero.

Theorderofafilteristhenumberofpolesitusesinthe
z-domain,andequivalentlythenumberofdelayelementsnec-
essarytoimplementit.

Bothpassbandandstopbandwillinpracticenothavegains
ofexactlyoneandzero,respectively,butmayshowseveral
deviationsfromtheseidealvalues,andtheseripplesmayhave
aspecifiedmaximumquotientbetweenthehighestandlowest
gain.
103

Therewillinpracticenotbeanabruptchangeofgainbetween
passbandandstopband,butatransitionbandwherethefre-
quencyresponsewillgraduallychangefromitspassbandtoits
stopbandvalue.
Thedesignercanthentradeoffconflictinggoalssuchasasmalltran-
sitionband,aloworder,alowrippleamplitude,orevenanabsenceof
ripples.
Designtechniquesformakingthesetradeoffsforanalogfilters(involv-
ingcapacitors,resistors,coils)canalsobeusedtodesigndigitalIIR
filters:
Butterworthfilters
Havenoripples,gainfallsmonotonicallyacrossthepassandtransition
band.Withinthepassband,thegaindropsslowlydownto1−
p
1/2
(−3dB).Outsidethepassband,itdropsasymptoticallybyafactor2
N
peroctave(N20dB/decade).
104
ChebyshevtypeIfilters
Distributethegainerroruniformlythroughoutthepassband(equirip-
ples)anddropoffmonotonicallyoutside.
ChebyshevtypeIIfilters
Distributethegainerroruniformlythroughoutthestopband(equirip-
ples)anddropoffmonotonicallyinthepassband.
Ellipticfilters(Cauerfilters)
Distributethegainerrorasequiripplesbothinthepassbandandstop-
band.Thistypeoffilterisoptimalintermsofthecombinationofthe
passband-gaintolerance,stopband-gaintolerance,andtransition-band
widththatcanbeachievedatagivenfilterorder.
AllthesefilterdesigntechniquesareimplementedintheMATLABSignalProcessingToolboxin
thefunctionsbutter,cheby1,cheby2,andellip,whichoutputthecoefficientsa
n
andb
n
ofthe
differenceequationthatdescribesthefilter.Thesecanbeappliedwithfiltertoasequence,or
canbevisualizedwithzplaneaspoles/zerosinthez-domain,withimpzasanimpulseresponse,
andwithfreqzasanamplitudeandphasespectrum.Thecommandssptoolandfdatool
provideinteractiveGUIstodesigndigitalfilters.
105
Butterworthfilterdesignexample
-1
0
1
-1
-0.5
0
0.5
1
Real Part
Imaginary Part
0
10
20
30
0
0.2
0.4
0.6
0.8
n (samples)
Amplitude
Impulse Response
0
0.5
1
-60
-40
-20
0
Normalized Frequency ( rad/sample)
Magnitude (dB)
0
0.5
1
-100
-50
0
Normalized Frequency ( rad/sample)
Phase (degrees)
order:1,cutofffrequency(−3dB):0.25×f
s
/2
106
Butterworthfilterdesignexample
-1
0
1
-1
-0.5
0
0.5
1
Real Part
Imaginary Part
0
10
20
30
-0.1
0
0.1
0.2
0.3
n (samples)
Amplitude
Impulse Response
0
0.5
1
-60
-40
-20
0
Normalized Frequency ( rad/sample)
Magnitude (dB)
0
0.5
1
-600
-400
-200
0
Normalized Frequency ( rad/sample)
Phase (degrees)
order:5,cutofffrequency(−3dB):0.25×f
s
/2
107
ChebyshevtypeIfilterdesignexample
-1
0
1
-1
-0.5
0
0.5
1
Real Part
Imaginary Part
0
10
20
30
-0.2
0
0.2
0.4
0.6
n (samples)
Amplitude
Impulse Response
0
0.5
1
-60
-40
-20
0
Normalized Frequency ( rad/sample)
Magnitude (dB)
0
0.5
1
-600
-400
-200
0
Normalized Frequency ( rad/sample)
Phase (degrees)
order:5,cutofffrequency:0.5×f
s
/2,pass-bandripple:−3dB
108
ChebyshevtypeIIfilterdesignexample
-1
0
1
-1
-0.5
0
0.5
1
Real Part
Imaginary Part
0
10
20
30
-0.2
0
0.2
0.4
0.6
n (samples)
Amplitude
Impulse Response
0
0.5
1
-60
-40
-20
0
Normalized Frequency ( rad/sample)
Magnitude (dB)
0
0.5
1
-300
-200
-100
0
100
Normalized Frequency ( rad/sample)
Phase (degrees)
order:5,cutofffrequency:0.5×f
s
/2,stop-bandripple:−20dB
109
Ellipticfilterdesignexample
-1
0
1
-1
-0.5
0
0.5
1
Real Part
Imaginary Part
0
10
20
30
-0.2
0
0.2
0.4
0.6
n (samples)
Amplitude
Impulse Response
0
0.5
1
-60
-40
-20
0
Normalized Frequency ( rad/sample)
Magnitude (dB)
0
0.5
1
-400
-300
-200
-100
0
Normalized Frequency ( rad/sample)
Phase (degrees)
order:5,cutofffrequency:0.5×f
s
/2,pass-bandripple:−3dB,stop-bandripple:−20dB
110
Exercise14DrawthedirectformIIblockdiagramsofthecausalinfinite-
impulseresponsefiltersdescribedbythefollowingz-transformsandwrite
downaformuladescribingtheirtime-domainimpulseresponses:
(a)H(z)=
1
1−
1
2
z
−1
(b)H

(z)=
1−
1
4
4
z
−4
1−
1
4
z
−1
(c)H
′′
(z)=
1
2
+
1
4
z
−1
+
1
2
z
−2
Exercise15(a)Performthepolynomialdivisionoftherationalfunction
giveninexercise14(a)untilyouhavefoundthecoefficientofz
−5
inthe
result.
(b)Performthepolynomialdivisionoftherationalfunctiongiveninexercise
14(b)untilyouhavefoundthecoefficientofz
−10
intheresult.
(c)Useitsz-transformtoshowthatthefilterinexercise14(b)hasactually
afiniteimpulseresponseanddrawthecorrespondingblockdiagram.
111
Exercise16Considerthesystemh:{x
n
}→{y
n
}withy
n
+y
n−1
=
x
n
−x
n−4
.
(a)DrawthedirectformIblockdiagramofadigitalfilterthatrealisesh.
(b)Whatistheimpulseresponseofh?
(c)Whatisthestepresponseofh(i.e.,h∗u)?
(d)Applythez-transformto(theimpulseresponseof)htoexpressitasa
rationalfunctionH(z).
(e)Canyoueliminateacommonfactorfromnumeratoranddenominator?
Whatdoesthismean?
(f)Forwhatvaluesz∈CisH(z)=0?
(g)HowmanypolesdoesHhaveinthecomplexplane?
(h)WriteHasafractionusingthepositionofitspolesandzerosanddraw
theirlocationinrelationtothecomplexunitcircle.
(i)Ifhisappliedtoasoundfilewithasamplingfrequencyof8000Hz,
sinewavesofwhatfrequencywillbeeliminatedandsinewavesofwhat
frequencywillbequadrupledintheiramplitude?
112
Randomsequencesandnoise
Adiscreterandomsequence{x
n
}isasequenceofnumbers
...,x
−2
,x
−1
,x
0
,x
1
,x
2
,...
whereeachvaluex
n
istheoutcomeofarandomvariablex
n
ina
correspondingsequenceofrandomvariables
...,x
−2
,x
−1
,x
0
,x
1
,x
2
,...
Suchacollectionofrandomvariablesiscalledarandomprocess.Each
individualrandomvariablex
n
ischaracterizedbyitsprobabilitydistri-
butionfunction
P
x
n
(a)=Prob(x
n
≤a)
andtheentirerandomprocessischaracterizedcompletelybyalljoint
probabilitydistributionfunctions
P
x
n
1
,...,x
n
k
(a
1
,...,a
k
)=Prob(x
n
1
≤a
1
∧...∧x
n
k
≤a
k
)
forallpossiblesets{x
n
1
,...,x
n
k
}.
113
Tworandomvariablesx
n
andx
m
arecalledindependentif
P
x
n
,x
m
(a,b)=P
x
n
(a)P
x
m
(b)
andarandomprocessiscalledstationaryif
P
x
n
1
+l
,...,x
n
k
+l
(a
1
,...,a
k
)=P
x
n
1
,...,x
n
k
(a
1
,...,a
k
)
foralll,thatis,iftheprobabilitydistributionsaretimeinvariant.
Thederivativep
x
n
(a)=P

x
n
(a)iscalledtheprobabilitydensityfunc-
tion,andhelpsustodefinequantitiessuchasthe

expectedvalueE(x
n
)=
R
ap
x
n
(a)da

mean-squarevalue(averagepower)E(|x
n
|
2
)=
R
|a|
2
p
x
n
(a)da

varianceVar(x
n
)=E[|x
n
−E(x
n
)|
2
]=E(|x
n
|
2
)−|E(x
n
)|
2

correlationCor(x
n
,x
m
)=E(x
n
x

m
)
RememberthatE()islinear,thatisE(ax)=aE(x)andE(x+y)=E(x)+E(y).Also,
Var(ax)=a
2
Var(x)and,ifxandyareindependent,Var(x+y)=Var(x)+Var(y).
114
Astationaryrandomprocess{x
n
}canbecharacterizedbyitsmean
value
m
x
=E(x
n
),
itsvariance
σ
2
x
=E(|x
n
−m
x
|
2
)=γ
xx
(0)

x
isalsocalledstandarddeviation),itsautocorrelationsequence
φ
xx
(k)=E(x
n+k
x

n
)
anditsautocovariancesequence
γ
xx
(k)=E[(x
n+k
−m
x
)(x
n
−m
x
)

]=φ
xx
(k)−|m
x
|
2
Apairofstationaryrandomprocesses{x
n
}and{y
n
}can,inaddition,
becharacterizedbyitscrosscorrelationsequence
φ
xy
(k)=E(x
n+k
y

n
)
anditscrosscovariancesequence
γ
xy
(k)=E[(x
n+k
−m
x
)(y
n
−m
y
)

]=φ
xy
(k)−m
x
m

y
115
Deterministiccrosscorrelationsequence
Fordeterministicsequences{x
n
}and{y
n
},thecrosscorrelationsequence
is
c
xy
(k)=

X
i=−∞
x
i+k
y
i
.
Afterdividingthroughtheoverlappinglengthofthefinitesequencesinvolved,c
xy
(k)canbe
usedtoestimate,fromafinitesampleofastationaryrandomsequence,theunderlyingφ
xy
(k).
MATLAB’sxcorrfunctiondoesthatwithoptionunbiased.
If{x
n
}issimilarto{y
n
},butlagslelementsbehind(x
n
≈y
n−l
),then
c
xy
(l)willbeapeakinthecrosscorrelationsequence.Itisthereforewidely
calculatedtolocateshiftedversionsofaknownsequenceinanotherone.
Thedeterministiccrosscorrelationsequenceisaclosecousinoftheconvo-
lution,withjustthesecondinputsequencemirrored:
{c
xy
(n)}={x
n
}∗{y
−n
}
ItcanthereforebecalculatedequallyeasilyviatheFouriertransform:
C
xy
(f)=X(f)Y

(f)
Swappingtheinputsequencesmirrorstheoutputsequence:c
xy
(k)=c
yx
(−k).
116
Equivalently,wedefinethedeterministicautocorrelationsequencein
thetimedomainas
c
xx
(k)=

X
i=−∞
x
i+k
x
i
.
whichcorrespondsinthefrequencydomainto
C
xx
(f)=X(f)X

(f)=|X(f)|
2
.
Inotherwords,theFouriertransformC
xx
(f)oftheautocorrelation
sequence{c
xx
(n)}ofasequence{x
n
}isidenticaltothesquaredam-
plitudesoftheFouriertransform,orpowerspectrum,of{x
n
}.
Thissuggests,thattheFouriertransformoftheautocorrelationse-
quenceofarandomprocessmightbeasuitablewayfordefiningthe
powerspectrumofthatrandomprocess.
WhatcanwesayaboutthephaseintheFourierspectrumofatime-invariantrandomprocess?
117
Filteredrandomsequences
Let{x
n
}bearandomsequencefromastationaryrandomprocess.
Theoutput
y
n
=

X
k=−∞
h
k
x
n−k
=

X
k=−∞
h
n−k
x
k
ofanLTIappliedtoitwillthenbeanotherrandomsequence,charac-
terizedby
m
y
=m
x

X
k=−∞
h
k
and
φ
yy
(k)=

X
i=−∞
φ
xx
(k−i)c
hh
(i),where
φ
xx
(k)=E(x
n+k
x

n
)
c
hh
(k)=
P

i=−∞
h
i+k
h
i
.
118
Inotherwords:
{y
n
}={h
n
}∗{x
n
}⇒

yy
(n)}={c
hh
(n)}∗{φ
xx
(n)}
Φ
yy
(f)=|H(f)|
2

xx
(f)
Similarly:
{y
n
}={h
n
}∗{x
n
}⇒

yx
(n)}={h
n
}∗{φ
xx
(n)}
Φ
yx
(f)=H(f)Φ
xx
(f)
Whitenoise
Arandomsequence{x
n
}isawhitenoisesignal,ifm
x
=0and
φ
xx
(k)=σ
2
x
δ
k
.
Thepowerspectrumofawhitenoisesignalisflat:
Φ
xx
(f)=σ
2
x
.
119
Applicationexample:
WhereanLTI{y
n
}={h
n
}∗{x
n
}canbeobservedtooperateon
whitenoise{x
n
}withφ
xx
(k)=σ
2
x
δ
k
,thecrosscorrelationbetween
inputandoutputwillrevealtheimpulseresponseofthesystem:
φ
yx
(k)=σ
2
x
h
k
whereφ
yx
(k)=φ
xy
(−k)=E(y
n+k
x

n
).
120
DFTaveraging
Theabovediagramsshowdifferenttypesofspectralestimatesofasequence
x
i
=sin(2πj×8/64)+sin(2πj×14.32/64)+n
i
withφ
nn
(i)=4δ
i
.
Leftisasingle64-elementDFTof{x
i
}(withrectangularwindow).The
flatspectrumofwhitenoiseisonlyanexpectedvalue.Inasinglediscrete
Fouriertransformofsuchasequence,thesignificantvarianceofthenoise
spectrumbecomesvisible.Italmostdrownsthetwopeaksfromsinewaves.
Aftercutting{x
i
}into1000windowsof64elementseach,calculatingtheir
DFT,andplottingtheaverageoftheirabsolutevalues,thecentrefigure
showsanapproximationoftheexpectedvalueoftheamplitudespectrum,
withaflatnoisefloor.Takingtheabsolutevaluebeforespectralaveraging
iscalledincoherentaveraging,asthephaseinformationisthrownaway.
121
Therightmostfigurewasgeneratedfromthesamesetof1000windows,
butthistimethecomplexvaluesoftheDFTswereaveragedbeforethe
absolutevaluewastaken.Thisiscalledcoherentaveragingand,because
ofthelinearityoftheDFT,identicaltofirstaveragingthe1000windows
andthenapplyingasingleDFTandtakingitsabsolutevalue.Thewindows
start64samplesapart.Onlyperiodicwaveformswithaperiodthatdivides
64arenotaveragedaway.Thisperiodicaveragingstepsuppressesboththe
noiseandthesecondsinewave.
Periodicaveraging
Ifazero-meansignal{x
i
}hasaperiodiccomponentwithperiodp,the
periodiccomponentcanbeisolatedbyperiodicaveraging:
¯x
i
=lim
k→∞
1
2k+1
k
X
n=−k
x
i+pn
Periodicaveragingcorrespondsinthetimedomaintoconvolutionwitha
Diraccomb
P
n
δ
i−pn
.Inthefrequencydomain,thismeansmultiplication
withaDiraccombthateliminatesallfrequenciesbutmultiplesof1/p.
122
Image,videoandaudiocompression
Structureofmodernaudiovisualcommunicationsystems:
signal
sensor+
sampling
perceptual
coding
entropy
coding
channel
coding
noise
channel
human
senses
display
perceptual
decoding
entropy
decoding
channel
decoding











123
Audio-visuallossycodingtodaytypicallyconsistsofthesesteps:

Atransducerconvertstheoriginalstimulusintoavoltage.

Thisanalogsignalisthensampledandquantized.
Thedigitizationparameters(samplingfrequency,quantizationlevels)arepreferably
chosengenerouslybeyondtheabilityofhumansensesoroutputdevices.

Thedigitizedsensor-domainsignalisthentransformedintoa
perceptualdomain.
Thisstepoftenmimicssomeofthefirstneuralprocessingstepsinhumans.

Thissignalisquantizedagain,basedonaperceptualmodelofwhat
levelofquantization-noisehumanscanstillsense.

Theresultingquantizedlevelsmaystillbehighlystatisticallyde-
pendent.Apredictionordecorrelationtransformexploitsthisand
producesalessdependentsymbolsequenceoflowerentropy.

Anentropycoderturnsthatintoanapparently-randombitstring,
whoselengthapproximatestheremainingentropy.
Thefirstneuralprocessingstepsinhumansareineffectoftenakindofdecorrelationtransform;
oureyesandearswereoptimizedlikeanyotherAVcommunicationssystem.Thisallowsusto
usethesametransformfordecorrelatingandtransformingintoaperceptuallyrelevantdomain.
124
Outlineoftheremaininglectures

Quickreviewofentropycoding

Transformcoding:techniquesforconvertingsequencesofhighly-
dependentsymbolsintoless-dependentlower-entropysequences.

run-lengthcoding

decorrelation,Karhunen-Lo`evetransform(PCA)

otherorthogonaltransforms(especiallyDCT)

Introductiontosomecharacteristicsandlimitsofhumansenses

perceptualscalesandsensitivitylimits

colourvision

humanhearinglimits,criticalbands,audiomasking

Quantizationtechniquestoremoveinformationthatisirrelevantto
humansenses
125

Imageandaudiocodingstandards

A/-lawcoding(digitaltelephonenetwork)

JPEG

MPEGvideo

MPEGaudio
Literature

D.Salomon:Aguidetodatacompressionmethods.
ISBN0387952608,2002.

L.Gulick,G.Gescheider,R.Frisina:Hearing.ISBN0195043073,
1989.

H.Schiffman:Sensationandperception.ISBN0471082082,1982.
126
Entropycodingreview–Huffman
Entropy:H=