Digital Signal Processing Markus Kuhn
Computer Laboratory http://www.cl.cam.ac.uk/teaching/0910/DSP/
Michaelmas 2009 – Part II
Signals
→ flow of information → measured quantity that varies with time (or position) signal received from a transducer → electrical (microphone (micr ophone,, thermo thermometer, meter, accel accelerome erometer, ter, antenn antenna, a, etc.) → electrical signal that controls a process voltag tage, e, current current,, tem temperat perature ure,, spee speed, d, . . . Continuous-time signals: vol
Discrete-time signals: daily minimu minimum/maxi m/maximum mum temperatur temperature, e, lap interva intervals ls in rac races, es, sampled sampled continu continuous ous signals signals,, . . . Electroni Electr onics cs (un (unlik likee opt optics ics)) can onl onlyy dea deall eas easily ily wit with h tim time-d e-depen ependen dentt sig signal nals, s, the theref refor oree spa spatia tiall signals, such as images, are typically first converted into a time signal with a scanning process (TV, fax, etc.). 2
Signal processing Signals may have to be transformed in order to
→ amplify or filter out embedded information → detect patterns → prepare the signal to survive a transmission channel → prevent interference with other signals sharing a medium → undo distortions contributed by a transmission channel → compensate for sensor deficiencies → find information encoded in a different domain To do so, we also need to measure, characterise, model and simulate trans→ methods mission missio n chann channels els tools that split common channels and transfor→ mathematical mations into easily manipulated building blocks 3
Analog electronics Passive net Passive netwo works rks (resi (resisto stors, rs, capa capacito citors, rs, inductances, crystals, SAW filters), non-li non -linea nearr element elementss (diodes, (diodes, . . . ), (roughly) linear operational amplifiers Advantages:
•
passive networks are highly linear over a very large dynamic range and large bandwidths
•
analog signa analog signal-p l-processin rocessingg circ circuits uits require little or no power
•
analog circuits cause little additional interference
R
U iin n
L
C
U out out
U in in U in in t u o
U
U out out 0
U in in
√
1 / LC
ω (= 2 f )
− U out 1 out = R
t
t
dU out out U out τ C d + out L −∞ dt
4
Digital signal processing Analog/digital and digital/analog converter, CPU, DSP, ASIC, FPGA.
Advantages:
→ noise is easy to control after initial quantization → highly linear (within limited dynamic range) → complex algorithms fit into a single chip → flexibility, parameters can easily be varied in software processing is insensitive to component tolerances, aging, → digital environmental conditions, electromagnetic interference But:
→ discrete-time processing artifacts (aliasing) → can require significantly more power (battery, cooling) → digital clock and switching cause interference
5
Typical DSP applications
→ communication systems
modulation/demodulati modulation/de modulation, on, chann channel el equalization, echo cancellation
→ consumer electronics
perceptual coding of audio and video on DVDs, speech synthesis, speech recognition
→ music
synthetic instruments, audio effects, noise reduction
→ medical diagnostics
magnetic-resonan magnetic-re sonance ce and ultraso ultrasonic nic imag im agin ing, g, co comp mput uter er to tomo mogr grap aphy hy,, ECG, EEG, MEG, AED, audiology
→ geophysics
seismology, oil exploration
→ astronomy experimental ntal physic physicss → experime → aviation → security
VLBI, speckle interfe interferomet rometry ry
sensor-data evaluation
radar, radio navigation
steganography, digital waterma steganography, watermarking, rking, biometric biome tric identi identificati fication, on, survei surveillanc llancee system sys tems, s, sig signal nalss inte intellig lligenc ence, e, ele elecctronic warfare
→ engineering
controll sys contro system tems, s, fea feature ture ext extrac ractio tion n for pattern recognition 6
Syllabus Signals and systems. Discrete sequences and systems, their types and properties. Linea Linearr timetime-inva invariant riant systems, systems, conv convoluti olution. on. Ha Harmoni rmonicc phasors are the eigen functions funct ions of linear time-inva time-invariant riant systems. systems. Revi Review ew of com complex plex arithme arithmetic. tic. Some examples from electronics, optics and acoustics. Fourier transform. Harmonic phasors as orthogonal base functions. Forms of the Fourier transform, convolution theorem, Dirac’s delta function, impulse combs in the time and frequency domain. Discrete sequences and spectra. Periodic sampling of continuous signals, periodic riod ic sig signal nals, s, ali aliasi asing, ng, sam sampli pling ng and rec recons onstru tructi ction on of lo low-p w-pass ass and ban band-p d-pass ass signals, spectral inversion. Discrete Fourier transform. Continuous versus discrete Fourier transform, symmetry, linearity, review of the FFT, real-valued FFT. Spectral estimation. Leakage and scalloping phenomena, windowing, zero padding. MATLAB: Some of the most important exercises in this course require writing small programs, preferably in MATLAB (or a similar tool), which is available on PWF computers. A brief MATLAB introdu int roducti ction on wa wass giv given en in Pa Part rt IB “Un “Unix ix Tools ools”. ”. Rev Review iew that befo before re the first exe exerci rcise se and also read the “Getting Started” section in MATLAB’s built-in manual. 7
Finite and infinite impulse-response filters. Properties of filters, implementation forms, window-based FIR design, use of frequency-inversion to obtain highpass filters, use of modulation to obtain band-pass filters, FFT-based convolution, polynomial representation, z -transform, -transform, zeros and poles, use of analog IIR design techniques (Butterworth, Chebyshev I/II, elliptic filters). Random sequences and noise. Random variables, stationary processes, autocorrelation, crosscorrelation, deterministic crosscorrelation sequences, filtered random sequences, white noise, exponential averaging. Correlation coding. Random vectors, dependence versus correlation, covariance, decorrelation,, matrix diagonal decorrelation diagonalisation, isation, eigen decompos dec omposition, ition, Karhunen-L Karhunen-Lo` o` eve transeve form, principal/independent component analysis. Relation to orthogonal transform coding using fixed basis vectors, such as DCT. Lossy versus lossless compression. What information is discarded by human senses sen ses and can be eli elimin minate ated d by encoders encoders?? Pe Perce rceptu ptual al sca scales les,, mas maskin king, g, spa spatia tiall resolution, colour coordinates, some demonstration experiments. Quantization, image and audio coding standards. A/µ-law coding, delta coding, JPEG photographic photographic still still-imag -imagee comp compressi ression, on, moti motion on compen compensatio sation, n, MPEG video encoding, MPEG audio encoding.
8
Objectives By the end of the course, you should be able to
→ → → → → → → → → →
apply basic properties of time-invariant linear systems understand sampling, aliasing, convolution, filtering, the pitfalls of spectral estim estimation ation explain the above in time and frequency domain representations use filter-design software visualise and discuss digital filters in the z -domain use the FFT for convolution, deconvolution, filtering implement, apply and evaluate simple DSP applications in MATLAB apply transforms that reduce correlation between several signal sources understand and explain limits in human perception that are exploited by lossy compression techniques provide a good overview of the principles and characteristics of several widely-used compression techniques and standards for audiovisual signals 9
Textbooks
→ R.G. Lyons: PrenticeHall, 2004. ( 45) Oppenheim, R.W. Schafer: → A.V.2nd ed., Prentice-Hall, 1999. ( 47) → J. Stein: Wiley, 2000. ( 74) → S.W. Smith: Newness, 2003. ( 40) → K. Steiglitz: Addison-Wesley, Understanding digital signal processing. £
Discrete-time signal process-
ing.
£
Digital signal processing – a computer science per-
spective.
£
Digital signal processing – a practical guide for engineers and scientists. £ A digital signal processing primer – with applications to digital audio and computer music.
1996. (£40)
→ Sanjit K. Mitra: McGraw-Hill, 2002. ( 38)
Digital signal processing – a computer-based
approach.
£
10
Sequences and systems A discrete sequence xn ∞ n=−∞ is a sequence of numbers
{ }
. . . , x−2 , x−1 , x0 , x1 , x2 , . . . where xn denotes the n-th number in the sequence (n ( n Z). A discrete sequence maps integer numbers onto real (or complex) numbers.
∈
We normally abbreviate xn ∞ n=−∞ to xn , or to xn n if the running index is not obvious. The notation is not well standardized. Some authors write x[n] instead of xn , others x(n).
{ }
{ }
{ }
Where a discrete sequence xn samples a continuous function x(t) as
{ }
xn = x(ts n) = x(n/f s ),
·
we call ts the sampling period and f s = 1/ts the sampling frequency . A discrete system T receives as input a sequence xn and transforms it into an output sequence yn = T xn :
{ }
. . . , x2 , x1 , x0 , x−1 , . . .
{ }
{ }
discrete system T
. . . , y2 , y1 , y0 , y−1 , . . . 11
Properties of sequences A sequence xn is
{ }
absolutely summable
square summable
⇔
∞
| |
xn <
| ∞
n=
−∞ ∞
2
⇔ x | <∞ −∞ ⇔ ∃k > 0 : ∀ n ∈ Z : x n
n=
periodic
n
= xn+k
A square-summable sequence is also called an energy signal , and ∞ xn 2 n=−∞
|
|
is its energy energy. Thi Thiss ter termin minolo ology gy reflects reflects tha thatt if U is a voltage supplied to a load resistor R, then P = U I = U 2 /R is the power consumed, and P (t) dt the energy. So even where we drop physical units (e.g., volts) for simplicity in calculations, it is still customary to refer to the squared values of a sequence as power and to its sum or integral over time as energy .
12
A non-square-summable sequence is a power signal if its average power 1 lim k→∞ 1 + 2 k
k
|
xn
n= k
−
|2
exists.
Special sequences Unit-step sequence: un =
0, n < 0 1, n 0
≥
Impulse sequence: δn
1, n = 0 = 0, n = 0 = un un−1
−
13
Types of discrete systems A causal system cannot look into the future: yn = f f ((xn , xn−1 , xn−2 , . . .) .) A memory-less system depends only on the current input value: yn = f f ((xn ) A delay system shifts a sequence in time: yn = xn−d T is a time-invariant system if for any d
{y } = T {x } ⇐⇒ {y − } = T {x − }. T is a linear system if for any pair of sequences {x } and {x′ } T {a · x + b · x′ } = a · T {x } + b · T {x′ }. n
n
n d
n d
n
n
n
n
n
n
14
Examples: The accumulator system n
yn =
xk
k=
−∞
is a causal, linear, time-invariant system with memory, as are the backward difference system
yn = xn
−x − , n 1
the M-point moving average system 1 yn = M
M 1
−
k =0
xn−k =
xn−M +1 +1 +
··· + x −
n 1
+ xn
M
and the exponential averaging system yn = α xn + (1
·
− α) · y −
n 1
=α
∞
k=0
(1
k
− α) · x − . n k
15
Examples for time-invariant non-linear memory-less systems: yn = x2n ,
yn = log2 xn ,
yn = max min 256 256x xn , 255 , 0
{ {⌊
⌋
} }
Examples for linear but not time-invariant systems: yn =
xn , n 0 = xn un 0, n < 0
≥
yn = x⌊n/4⌋ yn = xn (eω jn )
·
·ℜ
Examples for linear time-invariant non-causal systems: yn yn
1 = (xn−1 + xn+1 ) 2 9 sin(πk sin( πkω ω) = xn+k [0.5 + 0. [0. 0 .5 cos(π k/10)] k/ 10)] πkω πk ω k=−9
·
·
·
16
Constant-coefficient difference equations Of particular practical interest are causal linear time-invariant systems of the form xn
b0
yn
N
yn = b0 xn
·
− k =1
ak yn−k
·
−a1 −a2
Block diagram representation of sequence operations: x′n Addition: Multiplication by constant: Delay:
−a3
xn + x′n
xn xn xn
a
axn
z −1
xn−1
z −1 yn−1 z −1 yn−2 z −1 yn−3
The ak and bm are constant consta nt coefficien coefficients. ts. 17
or xn M
yn =
m=0
bm xn−m
·
z −1
xn−1 −1 xn−2 −1 xn−3 z z
b0
b1
xn
b0 z −1
N
M
k =0
·
b3 yn
or the combination of both:
ak yn−k =
b2
m=0
bm xn−m
·
b1
xn−1 z −1
b2
xn−2 z −1 xn−3
b3
1 a− 0
−a1 −a2 −a3
yn z −1 yn−1 z −1 yn−2 z −1 yn−3
The MATLAB function filter is an efficient implementation of the last variant. 18
Convolution All linear time-invariant (LTI) systems can be represented in the form yn =
∞
k=
−∞
ak xn−k
·
where ak is a suitably chosen sequence of coefficients. This operation over sequences is called convolution and defined as
{ }
{ p } ∗ {q } = {r } ⇐⇒ ∀n ∈ Z : r n
n
n
n
=
∞
k=
pk qn−k .
−∞
·
If yn = an xn is a representation of an LTI system T T ,, with yn = T xn , then we call the sequence an the impulse response of T T ,, because an = T δn .
{ } { }∗{ } { } { } { } { }
{ }
19
Convolution examples A
B
C
D
E
F
A∗B
A∗C
C∗A
A∗E
D∗E
A∗F
20
Properties of convolution For arbitrary sequences pn , qn , rn and scalars a, b:
{ }{ }{ }
→ Convolution is associative ({ p } ∗ {q }) ∗ {r } = { p } ∗ ({q } ∗ {r }) → Convolution is commutative { p } ∗ {q } = {q } ∗ { p } → Convolution is linear { p } ∗ {a · q + b · r } = a · ({ p } ∗ {q }) + b · ({ p } ∗ {r }) → The impulse sequence (slide 13) is neutral under convolution { p } ∗ {δ } = {δ } ∗ { p } = { p } shifting is equivalent to convolving with a shifted → Sequence impulse n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
{ p − } = { p } ∗ {δ − } n d
n
n d
21
Can all LTI systems be represented by convolution? Any sequence xn can be decomposed into a weighted sum of shifted impulse sequences:
{ }
∞
{x } = n
xk
k=
−∞
· {δ − } n k
Let’s see what happens if we apply a linear (∗) time-invariant(∗∗) system T to such a decomposed sequence: =
T xn
{ }
∞ ∗ ∞ ·{ − } · −∞ −∞ ∞ · { − } ∗ { } ∞
T
xk
( )
δn
k
k=
(
∗∗) =
xk
k=
=
−∞
=
k=
δn
{xn} ∗ T {δn}
xk T δn−k
T δn =
k
k=
q.e.d.
−∞
{
xk
}
· {δn−k } ∗ T {δn}
⇒ The impulse response T {δ } fully characterizes an LTI system. n
22
Exercise 1 What type of discrete system (linear/non-linear, time-invariant/ non-time-in non-t ime-invar variant, iant, causal/ causal/non-ca non-causal, usal, causal, memo memory-les ry-less, s, etc.) is: 3xn−1 + xn−2 xn−3
(a) yn = xn
| |
(e) yn =
(b) yn =
−xn−1 + 2xn − xn+1
(f)) yn = xn en/14 (f
(c) yn =
8 i=0
(d) yn =
xn−i
1 2 (x2n
+ x2n+1 )
· (g) yn = xn · un (h) yn =
∞
i=
−∞
xi δi−n+2
·
Exercise 2 Prove that convolution is (a) commutative and (b) associative.
23
Exercise 3 A finite-length sequence is non-zero only at a finite number of positions. If m and n are the first and last non-zero positions, respectively, then we call n m + 1 the the length of that sequence. What maximum length can the result of convolving two sequences of length k and l have?
−
Exercise 4 The length-3 sequence a0 = 3, a1 = 2, a2 = 1 is convolved with a second sequence bn of length 5. (a) Write down this linear operation as a matrix multiplication involving a b R5 , and a result vector c. matrix A, a vector b = (1, 0, 0, 2, 2) (b) Use MATLAB to multiply your matrix by the vector and compare the result with that of using the conv function. (c) Use the MATLAB facilities for solving systems of linear equations to undo the above convolution step.
{ }
−
∈
Exercise 5 (a) Find a pair of sequences an and bn , where each one bn contains at least three different values and where the convolution an results in an all-zero sequence. (b) Does every LTI system T have an inverse LTI system T −1 such that xn = T −1 T xn for all sequences xn ? Why?
{ }
{ }
{ }
{ }
{ }
{ }∗{ } 24
Direct form I and II implementations xn
1 a− 0
b0
z −1
b1
−a1
xn−1 z −1
b2
−a2
xn−2 z −1
b3
−a3
xn−3
yn
xn
z −1 yn−1 z −1 yn−2 z −1 yn−3
1 a− 0
−a1 =
−a2 −a3
b0 z −1
z −1
z −1
yn
b1
b2
b3
The block diagram representation of the constant-coefficient difference equation on slide 18 is called the direct form I implementation. The number of delay elements can be halved by using the commutativity of convolution to swap the two feedback loops, leading to the direct form II implementation of the same LTI system. These two forms are only equivalent with ideal arithmetic (no rounding errors and range limits). 25
Convolutio Conv olution: n: optic opticss example If a projective lens is out of focus, the blurred image is equal to the original image convolved with the aperture shape (e.g., a filled circle):
∗ Point-spread function h (disk, r =
h(x, y ) =
1 , r2 π
0,
=
as ): 2f
image plane
focal plane
x2 + y 2 r 2 x2 + y 2 > r 2
≤
a
Original image I , blurred image B = I h, i.e.
∗
B (x, y ) =
Z Z
I (x
s
− x′ , y − y′ ) · h(x′ , y′ ) · dx′ dy′ f
26
Convolution Convo lution:: elect electroni ronics cs example U in in
R
U in in
U in
√
2
U iin n
C
t u o
U out out
U
U out out 0 0
t
1 /RC /RC
ω (= 2 f ) f )
Any passive network (R,L,C (R,L,C ) convolves its input voltage U in in with an impulse response function h, leading to U out out = U in in h, that is U out out (t) =
∞
∗
U in in (t
−∞
− τ τ )) · h(τ τ )) · dτ
In this example: U in in
− U
out out
R
dU out out = C , dt
·
h(t) =
1 RC
·e
−t RC
, t 0 t<0
≥
0,
27
Why are sine waves useful? 1) Adding together sine waves of equal frequency, but arbitrary amplitude and phase, results in another sine wave of the same frequency: sin(ωt sin(ωt sin(ωt A1 sin( ωt + ϕ1 ) + A2 sin( ωt + ϕ2 ) = A sin( ωt + ϕ)
·
with
A =
·
·
A21 + A22 + 2A 2A1 A2 cos( cos(ϕ ϕ2
−ϕ ) 1
A1 sin ϕ1 + A2 sin ϕ2 tan ϕ = A1 cos ϕ1 + A2 cos ϕ2 Sine waves of any phase can be formed from sin and cos alone: A sin(ωt sin(ωt + ϕ) = a sin( sin(ωt ωt)) + b cos( cos(ωt ωt))
·
·
·
ωt
A2 sin(ϕ2 )
A2
A
·
A1
ϕ2 A2 cos(ϕ2 )
ϕ
A1 sin(ϕ1 )
·
ϕ1 A1 cos(ϕ1 )
·
·
√ with a = A · cos( cos(ϕ ϕ), b = A · sin( sin(ϕ ϕ) and A = a
2
+ b2 , tan ϕ = ab . 28
Note: Convolution of a discrete sequence xn with another sequence yn is nothing but adding together scaled and delayed copies of xn . (Think of yn decomposed into a sum of impulses.) If xn is a sampled sine wave of frequency f f ,, so is xn yn ! = Sine-wave sequences form a family of discrete sequences that is closed under convolution with arbitrary sequences.
{ }
{ }
{ }
{ }
{ } ⇒
{ }∗{ }
The same applies for continuous sine waves and convolution.
2) Sine waves are orthogonal to each other:
∞
sin(ω sin( ω1 t + ϕ1 ) sin( sin(ω ω2 t + ϕ2 ) dt “=” 0
·
−∞
⇐⇒
ω1 = ω2
∨
ϕ1
−ϕ
2
= (2 (2k k + 1)π /2 (k
∈ Z)
They can be used to form an orthogonal function basis for a transform. The term “orthogonal” is used here in the context of an (infinitely dimensional) vector space, where the “vectors” are functions of the form f : R R (or f : R C) and the scalar product ∞ is defined as f g = −∞ f (t) g (t) dt.
·
R
→
·
→
29
Why are exponential functions useful? Adding together two exponential functions with the same base z , but different differ ent scale facto factorr and offse offset, t, resul results ts in anothe anotherr exponenti exponential al functi function on with the same base: A1 z t+ϕ1 + A2 z t+ϕ2 = A1 z t z ϕ1 + A2 z t z ϕ2 = (A1 z ϕ1 + A2 z ϕ2 ) z t = A z t
·
·
· · · · · · · Likewise, if we convolve a sequence {x } of values
·
n
. . . , z −3 , z −2 , z −1 , 1, z , z 2 , z 3 , . . .
xn = z n with an arbitrary sequence hn , we get yn = z n yn =
∞
k=
−∞
xn−k hk =
·
∞
k=
−∞
{ }
z n−k hk = z n
·
∞
· k=
{ } { } ∗ { h }, z − · h = z · H (z ) k
k
n
n
−∞
where H(z) is independent of n. Exponential sequences are closed under convolution with arbitrary sequences. The same applies in the continuous case. 30
Why are complex numbers so useful? 1) They give us all n solutions (“roots”) of equations involving polynomials up to degree n (the “ 1 = j ” story). 2) They give us the “great unifying theory” that combines sine and exponential exponenti al funct functions: ions:
√−
1 jωt cos(ωt cos( e + e− jωt ωt)) = 2 1 jωt sin(ωt sin( ωt)) = e e− jωt 2j
or
−
1 jωt+ϕ cos(ωt cos( ωt + ϕ) = e + e− jωt−ϕ 2
or
cos(ωn + ϕ) = cos(ωn sin(ωn sin( ωn + ϕ) = Notation:
jωn+ϕ
ℜ(e ℑ(e
) = jωn+ϕ ) =
jω n
) e jϕ ] jω n ) e jϕ ]
ℜ[(e ℑ[(e
· ·
ℜ(a + jb) := a and ℑ(a + jb) := b where j2 = −1 and a, b ∈ R.
31
We can now represent sine waves as projections of a rotating complex vector. This allows us to represent sine-wave sequences as exponential sequences with basis e jω . A phase shift in such a sequence corresponds to a rotation of a complex vector. 3) Complex multiplication allows us to modify the amplitude and phase of a complex rotating vector using a single operation and value. Rotation of a 2D vector in (x, ( x, y)-form is notationally slightly messy, but fortunately j2 = 1 does exactly what is required here:
−
x3 y3
= =
x2 y2
·
x1 y1
x1 y2 + x2 y1
z1 = x1 + jy1 , z1 z2 = x1 x2
−y2 · x2 x1 x2 − y1 y2
z2 = x2 + jy2
(x3 , y3 ) ( y2 , x2 )
−
(x2 , y2 ) (x1 , y1 )
− y1y2 + j(x1 y2 + x2y1) 32
Complex phasors Amplitude and phase are two distinct characteristics of a sine function that are inconvenient to keep separate notationally. Complex functions (and discrete sequences) of the form [cos(ωt sin(ωt A e j(ωt+ϕ) = A [cos( ωt + ϕ) + j sin( ωt + ϕ)]
· · · = −1) are able to represent both amplitude and phase in
(where j2 one single algebraic object. Thanks to complex multiplication, we can also incorporate in one single factor both a multiplicative change of amplitude and an additive change of phase of such a function. This makes discrete sequences of the form xn = e jωn T ,, because for each ω , eigensequences with respect to an LTI system T there is a complex number (eigenvalue) H (ω ) such that T xn = H (ω )
{ }
· {x } n
In the notation of slide 30, where the argument of H is the base, we would write H (e (e jω ). 33
Recall:: Fouri Recall ourier er transform transform The Fourier integral transform and its inverse are defined as
F {g(t)}(ω)
= G(ω ) = α
1
=
F − {G(ω)}(t)
g (t )
∞
= β
−∞ ∞
g (t) e∓ jωt dt
·
G(ω ) e± jωt dω
·
−∞
where α and β are constants chosen such that αβ = 1/(2π ). ). Many equ Many equiva ivalen lentt fo forms rms of the Fourie Fourierr tra transfo nsform rm ar aree use used d in the literatu literature. re. The There re is no str strong ong j j − ωt ωt consensus on whether the forward transform uses e and the backwards transform e , or vice versa. We follow here those authors who set α = 1 and β = 1/(2π ), ), to keep the convolution theorem free of a constant prefactor; others prefer α = β = 1/ 2π , in the interest of symmetry.
√
The substitution ω = 2π f f leads to a form without prefactors:
F {h(t)}(f f ))
= H (f f )) =
1
=
F − {H (f f ))}(t)
h(t)
=
∞
−∞ ∞ −∞
h(t) e∓2
·
H (f f )) e±2
·
jf t
dt
jf t
df
π
π
34
Another notation is in the continuous case
F {h(t)}(ω) 1
jω
F − {H (e(e )}(t)
∞
= H (e (e jω ) =
jωt
· e−
dt
−∞ ∞ 1 = H (e (e jω ) · e jωt dω 2π −∞
h(t)
=
h(t)
and in the discrete-sequence case
F {h }(ω) n
1
jω
F − {H (e(e )}(t)
= H (e (e jω ) =
∞
hn e− jωn
·
n=
−∞
1 = H (e (e jω ) e jωn dω 2π −
hn
=
π
π
·
which treats the Fourier transform as a special case of the z -transform (to be introduced shortly). 35
Properties of the Fourier transform If x(t)
•−◦
X (f f ))
and
y (t)
•−◦
Y ((f Y f ))
are pairs of functions that are mapped onto each other by the Fourier transform, then so are the following pairs. Linearity: ax((t) + by ax by((t) aX (f f )) + bY bY ((f f ))
•−◦
Time scaling: x(at at))
1 f X a a
•−◦ | |
Frequency scaling:
1 t x a a
||
•−◦
X (af af ))
36
Time shifting: x(t
X (f f )) e−2
− ∆t) •−◦
·
jf ∆t
π
Frequency shifting: x(t) e2
·
j∆f t
•−◦
π
X (f
− ∆f f ))
Parseval’s theorem (total power):
∞
| −∞
x(t) 2 dt
|
=
∞
|
X (f f )) 2 df
|
−∞
37
Four ourier ier tra transf nsfo orm exa exampl mple: e: re rect ct and sin sincc The Fourier transform of the “rectangular function” rect(tt) = rect(
1 if t < 12 1 1 if t = 2 2 0 ot othe herw rwis isee
|| ||
is the “(normalized) sinc function” rect(tt)}(f f )) = F {rect(
1 2
− 12
sin π f f 2 jf t − e dt = = sinc(f sinc(f )) π
f f
π
and vice versa
F {sinc( sinc(tt)}(f f )) = rect(f rect(f )). Some noteworthy properties of these functions:
R ∞ rect(t) dt ∞ sinc(t) dt = 1 = −∞ • R −∞ • sinc(0) = 1 = rect(0) • ∀n ∈ Z \ {0} : sinc(k) = 0 38
Convolution theorem Continuous form:
F {(f ∗ g)()(tt)} F {f f ((t) · g(t)}
F {f f ((t)} · F {g(t)} F {f f ((t)} ∗ F {g(t)}
= =
Discrete form:
{x } ∗ { y } n
=
n
{z } ⇐⇒ n
X (e (e jω ) Y Y (e (e jω ) = Z (e (e jω )
·
Convolution in the time domain is equivalent to (complex) scalar multiplication in the frequency domain. Convolution in the frequency domain corresponds to scalar multiplication in the time domain. − jωr
− jωr
R x(s)y(r − s)ds ⇐⇒ R z(r)e dr = R R x(s)y(r − s)e dsdr = R x(s) R y(r − s)e drds = R x(s)e R y(r − s)e dr ds = R x(s)e R y(t)e dtds = R x(s)e ds · R y(t)e dt. (Same for P instead of R .) Proof: z (r ) = s s
s
− jωr
r
− jωs
r
− jωt
t
− jωs
s
s
− jω (r −s)
r
− jωs
s
r
t:=r−s
− jωt
t
39
Dirac’s delta function The continuous equivalent of the impulse sequence δn is known as Dirac’s delta function δ (x). It is a gen genera eraliz lized ed function function,, defi defined ned such that
{ }
δ (x) =
∞
0,
∞
x=0 , x=0
1
δ (x) dx = 1
−∞
0
x
and can be thought of as the limit of function sequences such as δ (x) = lim n
or
→∞
0, n/22, n/
δ (x) = lim n
→∞
n
√
π
|x| ≥ 1/n |x| < 1/n − n2 x2 e
The delta fun functi ction on is mat mathem hemati atical cally ly spea speakin king g not a fun functi ction, on, but a distribution, that is an expression that is only defined when integrated. 40
Some properties of Dirac’s delta function:
∞
f ((x)δ (x f
−∞ ∞
− a) dx
e±2
jf t
π
= f f ((a)
df = δ (t)
−∞ ∞ 1 e± jωt dω = δ (t) 2π −∞ Fourier transform:
F {δ(t)}(ω)
∞
=
δ (t) e− jωt dt = e0 = 1
−∞ ∞ 1 F −1{1}(t) = 2π 1 −∞
·
jωt
·e
dω
= δ (t)
http://mathworld.wolfram.com/DeltaFunction.html
41
Sine and cosine in the frequency domain cos(2π f f 0 t) =
1 2 e 2
F {cos(2 F {sin(2
jf 0 t
π
jf 0 t
f 0 t) (f f f )) =
} f t)}(f f f ))
π
0
=
1 2 jf 0 t 1 −2 e e 2j 2j 1 f 0 ) + δ (f + f 0 ) 2 j f 0 ) + δ (f + f 0 ) 2
sin(2π f f0 t) =
π
π
−
1 δ (f 2 j δ (f 2
− −
ℜ
1 2
1 2
ℑ
j
f 0
1 2
f
−f 0
−
π
ℜ 1 2
−f 0
1 + e−2 2
jf 0 t
π
ℑ
j
f 0
f
As any real-valued signal x(t) can be represented as a combination of sine and cosine functions, the spectrum of any real-valued signal will show the symmetry X (e jω ) = [X (e− jω )]∗ , where ∗ denotes the complex conjugate (i.e., negated imaginary part). 42
Fourier transform symmetries We call a function x(t) odd if x( t) = even if x( t) =
− −
−x(t) x(t)
and ∗ is the complex conjugate, such that (a (a + jb)∗ = (a Then
·
x(t) is real x(t) is imaginary x(t) is even x(t) is odd x(t) is real and even x(t) is real and odd x(t) is imaginary and even x(t) is imaginary and odd
⇔ ⇔ ⇔ ⇔ ⇔ ⇔ ⇔ ⇔
− jb).
X ( f f )) = [X (f f )] )]∗ X ( f f )) = [X (f f )] )]∗ X (f f )) is even X (f f )) is odd X (f f )) is real and even X (f f )) is imaginary and odd X (f f )) is imaginary and even X (f f )) is real and odd
− −
−
43
Exampl Exa mple: e: amp amplit litude ude modul modulati ation on Communication channels usually permit only the use of a given frequency interval, such as 300–3400 Hz for the analog phone network or 590–598 MHz for TV channel 36. Modulation with a carrier frequency f c shifts the spectrum of a signal x(t) into the desired band. Amplitude modulation (AM): y (t) = A cos(2π tf tf c ) x(t)
·
·
X (f )
Y (f )
=
∗ −f l
0 f l
f
−f c
f c
f
−f c
0
f c
f
The spectrum of the baseband signal in the interval f l < f < f l is shifted by the modulation to the intervals f c f l < f < f c + f l . How can such a signal be demodulated?
± −
−
±
44
Sampling using a Dirac comb The loss of information in the sampling process that converts a continuouss functi tinuou function on x(t) into a discrete sequence xn defined by
{ }
xn = x(ts n) = x(n/f s )
·
can be modelled through multiplying x(t) by a comb of Dirac impulses s(t) = ts
∞
·
δ (t
n=
−∞
− t · n) s
to obtain the sampled function x ˆ(t) = x(t) s(t)
·
The function x ˆ(t) now contains exactly the same information as the discrete sequence xn , but is still in a form that can be analysed using the Fourier transform on continuous functions.
{ }
45
The Fourier transform of a Dirac comb s(t) = ts
∞
·
δ (t
n=
−∞
− t · n) s
=
∞
e2
jnt/ts
π
n=
−∞
is another Dirac comb
F · · −
S (f f )) =
∞
ts
n=
ts
∞ ∞
−∞
n=
−∞
δ (t
δ (t
− t n) s
ts n) e
2π jf t
(f f )) =
∞
−
dt =
δ f
n=
−∞
−∞
s(t)
−2ts −ts
n ts
.
S (f )
0
ts
2ts
t
−2f s −f s
0
f s
2f s f 46
Sampling and aliasing sample cos(2π tf) cos(2π t(k⋅ f ± f)) s
0
Sampled at frequency f s , the function cos(2πtf cos(2 πtf )) cannot be distinguished from cos[2πt cos[2 πt((kf s f f )] )] for any k Z.
±
∈
47
Frequency-domain view of sampling x(t)
s(t)
· 0
x ˆ(t)
= ...
t
...
−1/f s
X (f )
0 1/f s
...
t
0
f
−1/f s
0 1/f s
t
ˆ (f ) X
S (f )
∗
...
= ...
... ...
−f s
f s
f
...
−f s 0
f s
f
Sampling a signal in the time domain corresponds in the frequency domain dom ain to con convol volving ving its spectrum spectrum wit with h a Dir Dirac ac comb. The resulti resulting ng copies of the original signal spectrum in the spectrum of the sampled signal are called “images”. 48
Nyquist limit and anti-aliasing filters If the (double-sided) bandwidth of a signal to be sampled is larger than the sampling frequency f s , the images of the signal that emerge during sampling may overlap with the original spectrum. Such an overlap will hinder reconstruction of the original continuous signal by removing the aliasing frequencies with a reconstruction filter . Therefore, it is advisable to limit the bandwidth of the input signal to the sampling frequency f s befo before re sampli sampling, ng, using an anti-aliasing filter . In the common case of a real-valued base-band signal (with frequency content down to 0 Hz), all frequencies f that occur in the signal with non-zero power should be limited to the interval f s /2 < f < f s /2. The upper limit f s /2 for the single-sided bandwidth of a baseband signal is known as the “Nyquist limit”.
−
49
Nyquist limit and anti-aliasing filters Without anti-aliasing filter
With anti-aliasing filter
single-sided bandwidth
X (f )
Nyquist limit = f s /2
X (f ) anti-aliasing anti-aliasing filter
0
f
double-sided double-sided bandwidth bandwidth
ˆ (f ) X
−2f s −f s
−f s ˆ (f ) X
0
f s
2f s
f
−2f s −f s
0
f s
f
reconstruction reconstruction filter
0
f s
2f s
f
Anti-aliasing and reconstruction filters both suppress frequencies outside f < f s /2.
||
50
Reconstruction of a continuous band-limited waveform The ideal anti-aliasing filter for eliminating any frequency content above f s /2 before sampling with a frequency of f s has the Fourier transform H (f f )) =
1 if f < 0 if f >
|| ||
f s 2 f s 2
= rect(t rect(ts f f )).
This leads, after an inverse Fourier transform, to the impulse response sin tf tf 1 · h(t) = f = · sinc tf tf t π
s
π
s
s
s
t ts
.
The original band-limited signal can be reconstructed by convolving this with the sampled signal x ˆ (t), which eliminates the periodicity of the frequency domain introduced by the sampling process: x(t) = h(t) xˆ(t)
∗
Note that sampling h(t) gives the impulse function: h(t) s(t) = δ (t).
·
51
Impulse response of ideal low-pass filter with cut-off frequency f s /2:
0
−3 −2.5 −2 −1.5 −1 −0.5
0 t⋅ f
0.5
1
1.5
2
2 .5
3
s
52
Reconstruction filter example sampled signal interpolation result scaled/shifted sin(x)/x pulses
1
2
3
4
5 53
Reconstruction filters The mathematically ideal form of a reconstruction filter for suppressing aliasing aliasi ng frequ frequencie enciess interpola interpolates tes the sample sampled d signa signall xn = x(ts n) back into the continuous waveform
·
x(t) =
∞
n=
−∞
xn
·
sin π (t ts n) . ts n) π (t
− · − ·
Choice of sampling frequency Due to causality and economic constraints, practical analog filters can only approximate such an ideal low-pass filter. Instead of a sharp transition between the “pass band” (< f s /2) and the “stop band” (> f s /2), they feature a “transition band” in which their signal attenuation gradually increases. The sampling frequency is therefore usually chosen somewhat higher than twice the highest frequency of interest in the continuous signal (e.g., 4 ). On the other hand, the higher the sampling frequency, the higher are CPU, power and memory requireme requi rements. nts. There Therefor fore, e, the choice of sampl sampling ing frequency frequency is a tradeoff between between signal quality, analog filter cost and digital subsystem expenses.
×
54
Exercise 6 Digit Digital-to al-to-analo -analogg converters converters cannot output Dirac pulses. Instead, for each sample, they hold the output voltage (approximately) constant, sta nt, until the next sample sample ar arriv rives. es. Ho How w can thi thiss behav behaviou iourr be mode modeled led mathematically as a linear time-invariant system, and how does it affect the spectrum of the output signal? Exercise 7 Many DSP systems use “oversampling” to lessen the requirements men ts on the design design of an anal analog og reconstr reconstruct uction ion filter. filter. They use (a fini finite te approximation of) the sinc-interpolation formula to multiply the sampling frequency f s of the initial sampled signal by a factor N before passing it to the dig digita ital-t l-to-a o-analo nalogg con conver verter ter.. Whi While le thi thiss req require uiress mo more re CP CPU U opera operatio tions ns and a faster D/A converter, the requirements on the subsequently applied analog reconstruction filter are much less stringent. Explain why, and draw schematic representations of the signal spectrum before and after all the relevant signal-processing steps. Exercise 8 Similarly, explain how oversampling can be applied to lessen the requirements on the design of an analog anti-aliasing filter. 55
Band-pass signal sampling Sampled signals can also be reconstructed if their spectral components remain entirely within the interval n f s /2 < f < (n + 1) f s /2 for some n N. (The baseband case discussed so far is just n = 0.)
·
∈
||
·
In this case, the aliasing copies of the positive and the negative frequencies will interleave instead of ove overla rlap, p, and can the theref refor oree be rem remove oved d aga again in wit with h a rec recons onstruc tructio tion n filt filter er wit with h the impulse impulse response sin π tf 2n + 1 tf s /2 cos 2π tf h(t) = f s tf s 4 tf s /2 π tf
·
X (f )
− 54 f s
n=2
„
«
= (n + 1) f s
sin π t(n + 1) f s π t(n + 1) f s
ˆ (f ) X
anti-aliasing filter
0
5 4 f s
f
−f s
−f s 2
tnf s tnf − nf s sin tnf . tnf s π
π
reconstruction filter
0
f s
2
f s
f
56
Exercise 9 Reconstructing a sampled baseband signal:
•
Generate a one second long Gaussian noise sequence rn MATLAB function randn) with a sampling rate of 300 Hz.
{ } (using
•
fir1(50, (50, 45/ 45/150) 150) function to design a finite impulse reUse the fir1 spons spo nsee lo loww-pa pass ss filter filter wi with th a cu cutt-off off freque frequency ncy of 45 Hz. Us Usee th thee filtfilt function in order to apply that filter to the generated noise signal, resulting in the filtered noise signal xn .
{ }
•
Then sample xn at 100 Hz by setting all but every third sample value to zero, resulting in sequence yn .
•
Generate another low-pass filter with a cut-off frequency of 50 Hz and apply it to yn , in order to interpola interpolate te the recons reconstruct tructed ed filtered noise signal zn . Mu Mult ltipl iplyy th thee re resu sult lt by three three,, to compens compensat atee th thee energy lost during sampling.
•
Plot xn , yn , and zn . Finally compare xn and zn .
{ }
{ }
{ } { }
{ }{ }
{ }
{ }
{ }
Why should the first filter have a lower cut-off frequency than the second? 57
Exercise 10 Reconstructing a sampled band-pass signal:
•
Generate a 1 s noise sequence rn , as in exercise 9, but this time use a sampling frequency of 3 kHz.
•
Apply to that a band-pass filter that attenuates frequencies outside the interval 31–44 Hz, which the MATLAB Signal Processing Toolbox function ch cheby eby2( 2(3, 3, 30 30, , [31 44 44]/1 ]/150 500) 0) will design for you.
•
Then sample the resulting signal at 30 Hz by setting all but every 100-th sample value to zero.
•
cheb eby2( y2(3, 3, 20, [3 [30 0 45] 45]/1 /1500 500) ) another band-pass Generate with ch filter for the interval 30–45 Hz and apply it to the above 30-Hzsampled sam pled signal, signal, to reconstr reconstruct uct the orig riginal inal signal. signal. (Y (You’ ou’llll hav havee to multiply it by 100, to compensate the energy lost during sampling.)
•
Plot all the produced sequences and compare the original band-pass signal and that reconstructed after being sampled at 30 Hz.
{ }
Why does the reconstructed waveform differ much more from the original if you reduce the cut-off frequencies of both band-pass filters by 5 Hz? 58
Spectrum of a periodic signal A signal x(t) that is periodic with frequency f p can be factored into a sing si ngle le period period x˙ (t) convolved with an impulse comb p(t). Th This is corr correesponds in the frequency domain to the multiplication of the spectrum of the single period with a comb of impulses spaced f p apart. x(t)
x˙ (t)
p(t)
= ...
∗
...
−1/f p 0
1/f p
t
t
t
1/f p
P (f )
= f p
...
−1/f p 0
˙ (f ) X
X (f )
−f p 0
...
·
f
...
f
...
−f p 0
f p
f 59
Spectrum of a sampled signal A signal x(t) that is sampled with frequency f s has a spectrum that is periodic with a period of f s . x(t)
s(t)
· 0
ˆ ( t) x
= ...
t
...
−1/f s
X (f )
...
t
0 1/f s
−1/f s
0
f
0 1/f s
t
ˆ (f ) X
S (f )
∗
...
= ...
...
−f s
f s
f
...
...
−f s 0
f s
f 60
Continuous vs discrete Fourier transform
• Sampling a continuous signal makes its spectrum periodic • A periodic signal has a sampled spectrum We sample a signal x(t) with f s , getting x ˆ(t). We take n consecutive samples of x ˆ(t) and repeat these periodically, getting a new signal ¨x(t) ¨ (f with period n/f s . It Itss spect spectru rum m X f )) is sampled (i.e., has non-zero value) at frequency intervals f s /n and repeats itself with a period f s . ¨ (f Now both x ¨(t) and its spectrum X f )) are finite vectors of length n. ¨ (f ) X
x ¨(t)
...
... ... f s−1 0 f s−1
−n/f s
t
n/f s
...
−f s/n 0
−f s
f s /n
f s
f 61
Discrete Fourier Transform (DFT) n 1
X k =
−
2π j ik n
xi e−
·
i=0
1 xk = n
n 1
−
·
F n =
F n
1 1 1 1 .. .
1
1
1 e−2 j 2 e−2 j 3 e−2 j .. . −1 1 e−2 j
π
n
n
π
π
n
n
π
π
n
π
n
n
x0 x1 x2 .. .
xn−1
=
n
π
n
·
1
2 e−2 j 4 e−2 j 6 e−2 j .. . 2( −1) e−2 j
π
n
X 0 X 1 X 2 .. .
X n−1
,
3
n
π
n
π
n
π
..
n
·
1
·
X 0 X 1 X 2 .. .
−1
e−2 j 2( −1) e−2 j 3( −1) e−2 j .. . ( −1)( −1) e−2 j π
n
n
π
n
n
π
n
n
.
···
n
1 F n∗ n
× n matrix
··· ··· ··· ···
e−2 j 6 e−2 j 9 e−2 j .. . 3( −1) e−2 j π
·
i=0
The n-point DFT multiplies a vector with an n
ik
X i e2π j n
n
π
n
n
X n−1
=
x0 x1 x2 .. .
xn−1
62
Discrete Fourier Transform visualized
·
x0 x1 x2 x3 x4 x5 x6 x7
=
X 0 X 1 X 2 X 3 X 4 X 5 X 6 X 7
The n-point DFT of a signal xi sampled at frequency f s contains in the elements X 0 to X n/ n/2 of the resulting frequency-domain vector the frequency components 0, f s /n /n,, 2f s /n /n,, 3f s /n /n,, . . . , f s /2, and contains in X n−1 downto X n/ corre rrespond sponding ing neg negati ative ve fre freque quenci ncies. es. Not Notee n/2 the co that for a real-valued input vector, both X 0 and X n/ n/2 will be real, too.
{ }
Why is there no phase information recovered at f s /2? 63
Inverse DFT visualized
1 8
·
·
X 0 X 1 X 2 X 3 X 4 X 5 X 6 X 7
=
x0 x1 x2 x3 x4 x5 x6 x7
64
Fast Fourier Transform (FFT) n 1
F { } − − · − − − · − − ·− F { } − − · F { F { } − − − · F { n 1 n xi i=0 k
xi e
=
2π j ik n
i=0
n 2
1
=
x2i e
ik 2π j n/ 2
+ e
k 2π j n
i=0
=
n 2
1
x2i+1 e
ik 2π j n/ 2
i=0
n 2
n 2
n 2
1 x2i i=0
k
n 2
1 x2i i=0
k
n 2
+ e
k 2π j n
+ e
k 2π j n
n 2
n 2
x2i+1 x2i+1
n 2
−1
} − } i=0
k
,
n 2
1 i=0
k
− n2
k< , k
n 2
≥ n2
The DFT over n-element vectors can be reduced to two DFTs over 2-element vectors plus n multiplications and n additions, leading to n/2-element n/ log2 n rounds and n log2 n additions and multiplications overall, compared to n2 for the equivalent matrix multiplication. A high-performance FFT implementation in C with many processor-specific optimizations and support for non-power-of-2 sizes is available at http://www.fftw.org/ . 65
Efficient real-valued FFT The symmetry properties of the Fourier transform applied to the discrete −1 = n xi n−1 have the form Fourier transform X i ni=0 i=0
{ }
F { } ∀i : xi = ℜ(xi) ⇐⇒ ∀i : X n−i = X i∗ ∀i : xi = j · ℑ(xi) ⇐⇒ ∀i : X n−i = −X i∗
These two symmetries, combined with the linearity of the DFT, allows us to calculate two real-valued n-point DFTs
{X i′}ni=0−1 = F n{x′i}ni=0−1
{X i′′}ni=0−1 = F n{x′′i }ni=0−1
simultaneously in a single complex-valued n-point DFT, by composing its input as xi = x′i + j x′′i and decomposing its output as 1 X i′ = (X i + X n∗−i ) 2
·
1 X i′′ = (X i 2
− X n∗−i)
To optimize the calculation of a single real-valued FFT, use this trick to calculate the two half-size real-value FFTs that occur in the first round. 66
Fast complex multiplication Calculating the product of two complex numbers as (a + jb) (c + jd) = (ac
·
− bd bd)) + j(ad j(ad + bc bc))
involves four (real-valued) multiplications and two additions. The alternative calculation (a + jb) (c + jd) = (α
·
− β ) + j(α j(α + γ )
with
α = a(c + d) β = d(a + b) γ = c(b a)
−
provides the same result with three multiplications and five additions. The latter may perform faster on CPUs where multiplications take three or more times longer than additions. This trick is most helpful on simple simplerr microco microcontroll ntrollers. ers. Speciali Specialized zed signal-processin signal-processing g CPUs (DSPs) feature 1-clock-cycle multipliers. High-end desktop processors use pipelined multipliers that stall where operations depend on each other. 67
FFT-based convolution Calculating the convolution of two finite sequences xi of lengths m and n via
{ }
m 1 i=0
− and {yi }n−1 i=0
min m 1,i
zi =
{ − }
j =max 0,i (n 1)
{ − − }
x j yi− j ,
0
·
≤i
takes mn multiplications. Can we apply the FFT and the convolution theorem to calculate the convolution faster, in just O(m log m + n log n) multiplications? 1
{z } = F − (F {x } · F {y }) i
i
i
There is obviously no problem if this condition is fulfilled:
{x } and {y } are periodic, with equal period lengths i
i
In this case, the fact that the DFT interprets its input as a single period of a periodic signal will do exactly what is needed, and the FFT and inverse FFT can be applied directly as above. 68
In the general case, measures have to be taken to prevent a wrap-over: −1
A
B
F [F(A)⋅F(B)]
A’
B’
F [F(A’)⋅F(B’)]
−1
Both sequences sequences are padded with zero values to a length of at least m + n 1. This ensures that the start and end of the resulting sequence do not overlap.
−
69
Zero padding is usually applied to extend both sequence lengths to the next higher power of two (2 ⌈log2 (m+n−1)⌉ ), which facilitates the FFT. With a causal sequence, simply append the padding zeros at the end. With a non-causal sequence, values with a negative index number are wrapped around the DFT block boundaries and appear at the right end. en d. In th this is ca case se,, ze zero ro-p -pad addi ding ng is ap appl plie ied d in th thee ce cente nterr of the bl block ock,, between the last and first element of the sequence. Thanks to the periodic nature of the DFT, zero padding at both ends has the same effect as padding only at one end. If both sequences can be loaded entirely into RAM, the FFT can be applied to them in one step. However, one of the sequences might be too large for that. It could also be a realtime waveform (e.g., a telephone signal) that cannot be delayed until the end of the transmission. In such cases, the sequence has to be split into shorter blocks that are separately convolved and then added together with a suitable overlap. 70
Each block is zero-padded at both ends and then convolved as before:
∗
∗
∗
=
=
=
The regions originally added as zero padding are, after convolution, aligned to overlap with the unpadded ends of their respective neighbour blocks. The overlapping parts of the blocks are then added together. 71
Deconvolution A signal u(t) was distorted by convolution with a known impulse response h(t) (e.g., through a transmission channel or a sensor problem). The “smeared” result s(t) was recorded. Can we undo the damage and restore (or at least estimate) u(t)?
∗
=
∗
=
72
The convolution theorem turns the problem into one of multiplication: s(t) =
u(t
− τ τ )) · h(τ τ )) · dτ
s = u h
F {s} F {u}
= =
u =
∗ F {u} · F {h} F {s}/F {h} F − {F {F{{s}/F {h}} 1
In practice, we also record some noise n(t) (quantization, etc.): c(t) = s(t) + n(t) =
u(t
− τ τ )) · h(τ τ )) · dτ + n(t) Problem – At frequencies f where F {h}(f f )) ap app pro roac ache hess ze zero ro,, the noise will be amplified (potentially enormously) during deconvolution: u˜ =
1
1
F − {F {F{{c}/F {h}} = u + F − {F {F{{n}/F {h}} 73
Typical workarounds: the Fourier transform of the impulse response, such that → Modify |F{{h}(f f ))| > ǫ for some experimentally chosen threshold ǫ. |F → Ifspectrum |F{{s}(f f ))| and the noise estimates of the signal spectrum |F |F{{n}(f f ))| can be obtained, then we can apply the |F “Wiener filter” (“optimal filter”)
2
|F {s}(f f ))| |F{ W ((f W f )) = |F{{s}(f f ))| + |F |F |F{{n}(f f ))| 2
2
before deconvolution: u˜ =
1
F − {W · F {c}/F {h}}
Exercise 11 Use MATLAB to deconvolve the blurred stars from slide 26. The files stars-blurred.png with the blurred-stars image and stars-psf.png with the impulse response (point-spre (point-spread ad function) are available available on the course-materia course-materiall we web b page. You may find the MATLAB functions imread , double, imagesc , circshift, fft2, ifft2 of use. Try different ways to control the noise (see above) and distortions near the margins (windowing). [The MATLAB MATLAB image processing toolbox provides provides ready-made ready-made “pro “professio fessional” nal” functions functions forr suc such h tas tasks. ks. Do not use the these, se, except except perdeconvwnr, deconvreg, deconvlucy , edgetaper, fo haps to compare their outputs with the results of your own attempts.] 74
Spectral estimation Sine wave 4×f /32
Discrete Fourier Transform
s
1
15 10
0
5 −1
0
10
20
30
0
0
Sine wave 4.61×f /32
10
20
30
Discrete Fourier Transform
s
1
15 10
0
5 −1
0
10
20
30
0
0
10
20
30 75
We introduced the DFT as a special case of the continuous Fourier transform, where the input is sampled and periodic . If the input is sampled, but not periodic, the DFT can still be used to calculate an approximation of the Fourier transform of the original contin con tinuou uouss sig signal nal.. Ho Howe wever ver,, the there re are two two effe effects cts to con consid sider. er. The Theyy are particularly visible when analysing pure sine waves. Sine waves whose frequency is a multiple of the base frequency (f ( f s /n /n)) of the DFT are identical to their periodic extension beyond the size of the DFT. They are, therefore, represented exactly by a single sharp peak in the DFT. All their energy falls into one single frequency “bin” in the DFT result. Sine waves with other frequencies, which do not match exactly one of the output frequency bins of the DFT, are still represented by a peak at the output bin that represents the nearest integer multiple of the DFT’s base frequency. However, such a peak is distorted in two ways:
→ Its amplitude is lower (down to 63.7%). → Much signal energy has “leaked” to other frequencies.
76
35 30 25 20 15 10 5 0 0
5
16 10
15
15.5
20
25
30
15
input freq.
DFT index
The leakage of energy to other frequency bins not only blurs the estimated spectrum. tru m. The peak am amplit plitude ude also changes changes significan significantly tly as the frequenc frequencyy of a ton tonee changes from that associated with one output bin to the next, a phenomenon known as scalloping . In the above graphic graphic,, an input sine wave gradua gradually lly changes changes from the frequency of bin 15 to that of bin 16 (only positive frequencies shown). 77
Windowing Sine wave
Discrete Fourier Transform 300
1
200 0 100 −1 0
200
400
0
0
200
400
Sine wave multiplied with window function Discrete Fourier Transform 100 1
0
50
−1
0
0
200
400
0
200
400 78
The reason for the leakage and scalloping losses is easy to visualize with the help of the convolution theorem: The operation of cutting a sequence of the size of the DFT input vector out of a longer original signal (the one whose continuous Fourier spectrum we try to estimate) is equivalent to multiplying this signal with a rectangular function. This destroys all information and continuity outside the “window” that is fed into the DFT. Multiplication with a rectangular window of length T in the time domain is equivalent to convolution with sin(πf T )/(πf T ) in the frequency domain. The subsequent interpretation of this window as a periodic sequence by the DFT leads to sampling of this convolution result (sampling meaning multiplication with a Dirac comb whose impulses are spaced f s /n apart). Where the window length was an exact multiple of the original signal period, sampling of the sin(πf T )/(πf T ) curve leads to a single Dirac pulse, and the windowing causes no distortion. In all other cases, the effects of the convolution become visible in the frequency domain as leakage and scalloping losses. 79
Some better window functions 1 0.8 0.6 0.4 0.2 Rectangular window Triangular window Hann window Hamming window
0
0
0.2
0.4
0.6
0.8
1
All these functions are 0 outside the interval [0,1]. 80
Rectangular window (64−point)
Triangular window
) 20 B d ( 0 e d u t i n −20 g a M−40
) 20 B d ( 0 e d u t i n −20 g a M−40
−60
−60
) 20 B d ( 0 e d u t i n −20 g a M−40
) 20 B d ( 0 e d u t i n −20 g a M−40
−60
−60
0 0.5 1 Normalized Frequency (×π rad/sample) Hann window
0 0.5 1 Normalized Frequency (×π rad/sample) Hamming window
0 0.5 1 Normalized Frequency (×π rad/sample)
0 0.5 1 Normalized Frequency (×π rad/sample)
81
Numerous alternatives to the rectangular window have been proposed thatt red tha reduce uce leakage leakage and sca scallo llopin pingg in spect spectral ral estimatio estimation. n. The These se ar aree vectors multiplied element-wise with the input vector before applying the DFT to it. The Theyy all force force the sig signal nal amplitud amplitudee smoot smoothly hly down down to zero at the edge of the window, thereby avoiding the introduction of sharp jumps in the signal when it is extended periodically by the DFT. −1 are: Three examples of such window vectors wi ni=0 Triangular window (Bartlett window):
{ }
wi = 1
− − 1
i n/22 n/
Hann window (raised-cosine window, Hanning window): wi = 0.5
Hamming window: wi = 0.54
− 0.5 × cos
− 0.46 × cos
2π
i
n
2π
− 1
i
n
−1
82
Zero padding increases DFT resolution The two figures below show two spectra of the 16-element sequence si = cos(2π 3i/16) + cos(2π 4i/16),
·
i
·
∈ {0, . . . , 15}.
The left plot shows the DFT of the windowed sequence xi = si wi ,
i
·
∈ {0, . . . , 15}
and the right plot shows the DFT of the zero-padded windowed sequence x′i =
where wi = 0.54
si wi , i i 0,
·
∈ {0, . . . , 15} ∈ {16, . . . , 63}
− 0.46 × co coss (2πi/ 15) is the Hamming window.
DFT without zero padding
DFT with 48 zeros appended to window 4
4 2
2
0
0
0
5
10
15
0
20
40
60 83
Applyi Appl ying ng th thee di disc scre rete te Fou ouri rier er tr tran ansf sfoorm to an n-el -eleme ement nt lon longg rea reallvalued sequence leads to a spectrum consisting of only n/ n/22 + 1 discre discrete te frequencies. Since the resulting spectrum has already been distorted by multiplying the (hypothetically longer) signal with a windowing function that limits its length to n non-zero values and forces the waveform smoothly down to zero at the window boundaries, appending further zeros outside the window will not distort the signal further. The frequency resolution of the DFT is the sampling frequency divided by the block size of the DFT. Zero padding can therefore be used to increase the frequency resolution of the DFT. Note that zero padding does not add any additional information to the sign si gnal al.. Th Thee spe spect ctru rum m ha hass al alre read adyy bee been n “l “looww-pa pass ss fil filte tere red” d” by bei being ng convolved with the spectrum of the windowing function. Zero padding in the time domain merely samples this spectrum blurred by the windowing step at a higher resolution, thereby making it easier to visually distinguish spectral lines and to locate their peak more precisely. 84
Frequency inversion In order to turn the spectrum X (f ) of a real-valued signal xi sampled at f s into an inverted spectrum X ′ (f ) = X (f s /2 f ), we merely have to shift the periodic spectrum by f s /2: X ′ (f ) X (f )
−
= ...
∗
... ...
−f s
0
f s
...
f
0
−f s
− f 2
f
f s
s
0
f s 2
f
This can be accom accomplished plished by multiplying multiplying the sampl sampled ed sequen sequence ce xi with yi = fs t = cos π i, which is nothing but multiplication with the sequence cos π f . . . , 1, 1, 1, 1, 1, 1, 1, 1, . . .
−
−
−
−
So in order to design a discrete high-pass filter that attenuates all frequencies f outside the range f c < f < f s /2, we merely have to design a low-pass filter that attenuates all frequencies outside the range f c < f < f c , and then multiply every second value of its impulse response with 1.
||
−
−
85
Window-based design of FIR filters Recall that the ideal continuous low-pass filter with cut-off frequency f c has the frequency characteristic H (f f )) =
1 if f < f c = rect 0 if f > f c
|| ||
f 2f c
and the impulse response sin2π tf tf c h(t) = 2f c = 2f c sinc(2 sinc(2f f c t). 2π tf tf c
·
·
Sampling this impulse response with the sampling frequency f s of the signal to be processed will lead to a periodic frequency characteristic, that matches the periodic spectrum of the sampled signal. There are two problems though:
→ the impulse response is infinitely long → this filter is not causal, that is h(t) = 0 for t < 0
86
Solutions:
→ Make the impulse response finite by multiplying the sampled h(t) with a windowing function → Make the impulse response causal by adding a delay of half the window size The impulse response of an n-th order low-pass filter is then chosen as hi = 2f c /f s
(i − n/ n/2) 2)f f /f ] · sin[2 · w 2 (i − n/ n/2) 2)f f /f π
c
π
c
s
s
i
where wi is a windowing sequence, such as the Hamming window
{ }
wi = 0.54
− 0.46 × cos (2 (2πi/n πi/n))
with wi = 0 for i < 0 and i > n. Note that for f c = f s /4, we have hi = 0 for all even values of i. The Theref refor ore, e, this special case case requires requir es only half the number of multip multiplicati lications ons during the convolution. convolution. Such “half-band” “half-band” FIR filters are used, for example, as anti-aliasing filters wherever a sampling rate needs to be halved. 87
FIR low-pass filter design example Impulse Response 0.3
t 1 r a P 0.5 y r a 0 n i g −0.5 a m I −1
e d u t i l p m A
30
−1
0.2 0.1 0
−0.1 0
0 1 Real Part
0
10 20 n (samples)
30
0
) s e e r g −500 e d ( e −1000 s a h P
) B d ( e −20 d u t i n −40 g a M
−60
0 0.5 1 Normalized Frequency (×π rad/sample)
−1500
order: n = 30, cutoff frequency ( 6 dB): f c = 0.25
−
0 0.5 1 Normalized Frequency (×π rad/sample)
× f s/2, window: Hamming
88
We truncate the ideal, infinitely-long impulse response by multiplication with a window sequence. In the frequency domain, this will convolve the rectangular frequency response of the ideal lowpass filter with the freque frequency ncy characteris characteristic tic of the window. window. The width of the main lobe determ determines ines the width of the transition band, and the side lobes cause ripples in the passband and stopband.
Converting a low-pass into a band-pass filter To obtain a band-pass filter that attenuates all frequencies f outside the range f l < f < f h , we first design a low-pass filter with a cut-off frequency (f (f h f l )/2 and multiply its impulse response with a sine wave of frequency (f ( f h + f l )/2, before applying the usual windowing:
−
hi = (f h
n/2)( n/ 2)(f f − f )/f ] − f )/f · sin[ (i(i−−n/ · cos[ n/2)( 2)(f f − f )/f π
l
s
h
π
h
l
l
s
(f h + f l )] wi
·
π
s
H (f )
= −f h −f l
0 f l
f h
f
∗ − f −2 f h
l
f h −f l 2
f
− f +2 f h
l
0
f h +f l 2
f 89
Exercise 12 Explain the difference between the DFT, FFT, and FFTW. Exercise 13 Push-button telephones use a combination of two sine tones to signal, which button is currently being pressed: 697 770 852 941
Hz Hz Hz Hz
1209 Hz 1 4 7 *
1336 Hz 2 5 8 0
1477 Hz 3 6 9 #
1633 Hz A B C D
(a) You receive a digital telephone signal with a sampling frequency of 8 kHz. You cut a 256-sample window out of this sequence, multiply it with a windowing function and apply a 256-point DFT. What are the indices where ) will show the highest amplitude if the resulting vector (X 0 , X 1 , . . . , X255 butto but ton n 9 wa wass pushed pushed at the time time of of the reco recordin rding? g? (b) Use MATLAB to determine, which button sequence was typed in the touch tones recorded in the file touchtone.wav on the course-material web page. 90
Polynomial representation of sequences We can represent sequences xn as polynomials:
{ }
∞
X (v ) =
xn v n
n=
−∞
Example of polynomial multiplication: (1 + 2v + 3v 2 )
·
(2 + 1 v )
2 + 4v + 6v 2 + 1v + 2v 2
+ 3v 3
= 2 + 5v + 8v 2
+ 3v 3
Compare this with the convolution of two sequences (in MATLAB): conv([1 conv([ 1 2 3], [2 1]) equals [2 5 8 3] 3] 91
Convolution of sequences is equivalent to polynomial multiplication:
{h } ∗ { x } n
n
=
H (v ) X (v ) =
·
{y } ⇒ ↓ ↓
k=
−∞
hk xn−k
·
· · · ∞
n=
=
yn =
n
∞
∞
∞
hn v n
−∞
n=
∞
n=
−∞ k=−∞
xn v n
−∞
hk xn−k v n
Note how the Fourier transform of a sequence can be accessed easily from its polynomial form: X (e (e− jω ) =
∞
n=
xn e− jωn
−∞ 92
Example of polynomial division: 1 1
= 1 + av + a2 v 2 + a3 v 3 +
− av 1
−
··· =
∞
an v n
n=0
1 + av + a2 v 2 + av 1 1 av av av a2 v 2 a2 v 2 a2 v 2 a3 v 3
···
−
−
−
···
Rational functions (quotients of two polynomials) can provide a convenient closed-form representations for infinitely-long exponential sequences, in particular the impulse responses of IIR filters. 93
The z -transform The z -transform of a sequence xn is defined as:
{ }
X (z ) =
∞
xn z −n
n=
−∞
Note that is differs only in the sign of the exponent from the polynomial representation discussed on the preceeding slides.
Recall that the above X (z ) is exactly the factor with which an exponentiall seque nentia sequence nce z n is multiplied, if it is convolved with xn :
{ }
n
{ }
{z } ∗ {x } = {y } ⇒
n
yn =
n
∞
k=
−∞
z n−k xk = z n
∞
· k=
−∞
z −k xk = z n X (z )
·
94
The z -transform defines for each sequence a continuous complex-valued surface over the complex plane C. For finite sequences, its value is always defined across the entire complex plane. For infinite sequences, it can be shown that the z -transform converges only for the region
xn+1 xn+1 lim < z < lim n→∞ n→−∞ xn xn
||
The z -transform identifies a sequence unambiguously only in conjunction with a given region of other er wo words, rds, there exist different different seq sequen uences ces,, tha thatt hav havee the same exp expres ressio sion n as convergence . In oth their z -transform, but that converge for different amplitudes of z .
The z -transform is a generalization of the Fourier transform, which it contains on the complex unit circle ( z = 1):
||
jω
F {x }(ω) = X (e(e n
)=
∞
xn e− jωn
n=
−∞ 95
The z -transform of the impulse response hn of the causal LTI system defined by
xn
{ }
k
m
l=0
al yn−l =
·
with yn = hn rational function
l=0
z −1
b0 b1
xn−1
bl xn−l
·
z −1
···
{ } { } ∗ {x } is the
z −1
n
··· bm
xn−m
1 a− 0
−a1 ··· −a
k
yn z −1 yn−1 z −1
··· z −1 yn−k
+ bm z −m b0 + b1 z −1 + b2 z −2 + H (z ) = a0 + a1 z −1 + a2 z −2 + + ak z − k (bm = 0, ak = 0) which can also be written as
··· ···
H (z ) =
zk zm
m m l l=0 bl z . k k l l=0 al z
− −
H (z ) has m zeros and k poles at non-zero locations in the z plane, plus k m zeros (if k > m) or m k poles (if m > k ) at z = 0.
−
−
96
This function can be converted into the form m
b0 H (z ) = a0
m
·
(1
l=1 k
(1
l=1
1
− c · z− ) l
1
− d · z− ) l
b0 k−m = z a0
·
·
(z
−c)
(z
−d)
l
l=1 k
l=1
l
where the cl are the non-zero positions of zeros (H ( H (cl ) = 0) and the dl are the non-zero positions of the poles (i.e., z dl H (z ) ) of H (z ). Exc Except ept for for a con consta stant nt factor, factor, H (z ) is entirely characterized by the position of these zeros and poles.
→ ⇒|
| →∞
As with the Fourier transform, convolution in the time domain corresponds to complex multiplication in the z -domain:
{x } •−◦ X (z), {y } •−◦ Y Y ((z) ⇒ {x } ∗ {y } •−◦ X (z) · Y Y ((z) n
n
n
n
Delaying a sequence by one corresponds in the z -domain to multiplication with z −1 : xn−∆n X (z ) z −∆n
{
} •−◦
·
97
2 1.75 1.5 1.25
|
z
) (
H |
1
0.75 0.5 0.25 0 1 0.5 0 −0.5 imaginary
−1
−1
0
−0.5
0.5
1
real
This example is an amplitude plot of xn 0.8 H (z ) =
1
0.8 0.8z = 0.2 z −1 z 0.2
− ·
−
which features a zero at 0 and a pole at 0. 0 .2.
yn
0.2
z −1 yn−1 98
H (z ) = z−z0.7 = 1−0.17·z−1 t r a 1 P y r a 0 n i g a m−1 I
z Plane
Impulse Response 1
e d u t i l 0.5 p m A
−1
0 1 Real Part
0
0
10 20 n (samples)
30
H (z ) = z−z0.9 = 1−0.19·z−1 t r a 1 P y r a 0 n i g a m−1 I
z Plane
Impulse Response 1
e d u t i l 0.5 p m A
−1
0 1 Real Part
0
0
10 20 n (samples)
30 99
H (z ) = z−z 1 = 1−1z−1 t r a 1 P y r a 0 n i g a m−1 I
z Plane
Impulse Response 1
e d u t i l 0.5 p m A
−1
0 1 Real Part
0
0
10 20 n (samples)
30
H (z ) = z−z1.1 = 1−1.11·z−1 t r a 1 P y r a 0 n i g a m−1 I
z Plane
Impulse Response 20
e d u t i l 10 p m A
−1
0 1 Real Part
0
0
10 20 n (samples)
30 100
H (z ) = t r a 1 P y r a 0 n i g a m−1 I
− ·
· − ·
1 = 1−1.8 co cos( s(π/6)z −1 +0.92 ·z −2
z Plane
Impulse Response e d u t i l p m A
2
−1
H (z ) = t r a 1 P y r a 0 n i g a m−1 I
z2 (z 0.9 e jπ/6 ) (z 0.9 e− jπ/6 )
−
· −
0
−2
0 1 Real Part
z2 (z e jπ/6 ) (z e− jπ/6 )
2
0
30
1 = 1−2 co cos( s(π/6)z −1 +z −2
z Plane
Impulse Response e d u t i l p m A
2
−1
10 20 n (samples)
5 0
−5
0 1 Real Part
0
10 20 n (samples)
30 101
H (z ) = t r a 1 P y r a 0 n i g a m−1 I
− ·
· − ·
1 = 1−1.8 co = cos( s(π/2)z −1 +0.92 ·z −2
z Plane
−1
0 1 Real Part
z z +1
=
·
1 0
−1
0
10 20 n (samples)
30
1 1+z −1
z Plane
Impulse Response e d u t i l p m A
−1
1 1+0.92 z −2
Impulse Response e d u t i l p m A
2
H (z ) = t r a 1 P y r a 0 n i g a m−1 I
z2 (z 0.9 e jπ/2 ) (z 0.9 e− jπ/2 )
0 1 Real Part
1 0
−1
0
10 20 n (samples)
30 102
IIR Filter design techniques The design of a filter starts with specifying the desired parameters:
→ The is the frequency range where we want to approximate a gain of one. is the frequency range where we want to approx→ The imate a gain of zero. of a filter is the number of poles it uses in the → The z -domain, and equivalently the number of delay elements necpassband
stopband
order
essary to implement it.
→ Both passband and stopband will in practice not have gains of ex exac actl tlyy on onee an and d ze zero ro,, re respe spect ctive ively ly,, bu butt ma mayy sh shoow se seve vera rall deviations from these ideal values, and these ripples may have a specified maximum quotient between the highest and lowest gain. 103
will in practice not be an abrupt change of gain between → There passband and stopband, but a where the fretransition band
quency response will gradually change from its passband to its stopband value. The designer can then trade off conflicting goals such as a small transition band, a low order, a low ripple amplitude, or even an absence of ripples. Design techniques for making these tradeoffs for analog filters (involving capacitors, resistors, coils) can also be used to design digital IIR filters:
Butterworth filters Have no ripples, gain falls monotonically across the pass and transition band. Within the passband, the gain drops slowly down to 1 1/2 ( 3 dB). Outside the passband, it drops asymptotically by a factor 2 N per octave (N (N 20 dB/decade).
−
·
−
104
Chebyshev type I filters Distribute the gain error uniformly throughout the passband (equiripples) and drop off monotonically outside.
Chebyshev type II filters Distribute the gain error uniformly throughout the stopband (equiripples) and drop off monotonically in the passband.
Elliptic filters (Cauer filters) Distribute the gain error as equiripples both in the passband and stopband. This type of filter is optimal in terms of the combination of the passband-gain tolerance, stopband-gain tolerance, and transition-band width that can be achieved at a given filter order. All these filter design techniques are implemented in the MATLAB Signal Processing Toolbox in the functions butter , cheby1, cheby2 , and ellip, which output the coefficients an and bn of the difference equation that describes the filter. These can be applied with filter to a sequence, or can be visualized with zplane as poles/zeros in the z -domain, with impz as an impulse response, and with freqz as an amp amplitu litude de and phase spectrum. spectrum. The comman commands ds sptool and fdatool provide interactive GUIs to design digital filters. 105
Butterworth filter design example Impulse Response 0.8
1
t r a P 0.5 y r a 0 n i g a −0.5 m I
e 0.6 d u t i l 0.4 p m A0.2
−1
−1
0 Real Part
0
1
0
) s e e r g e d ( e s a h P
) B d ( e −20 d u t i n −40 g a M
−60
0 0.5 1 Normalized Frequency (×π rad/sample)
order: 1, cutoff frequency ( 3 dB): 0.25
−
× f s /2
0
10 20 n (samples)
30
0
−50
−100
0 0.5 1 Normalized Frequency (×π rad/sample)
106
Butterworth filter design example Impulse Response 0.3
1
t r a P 0.5 y r a 0 n i g a −0.5 m I
e d u t i l p m A
−1
−1
0 Real Part
0.1 0
−0.1 0
1
0
0.2
10 20 n (samples)
30
0
) s e e r g −200 e d ( e −400 s a h P
) B d ( e −20 d u t i n −40 g a M
−60
0 0.5 1 Normalized Frequency (×π rad/sample)
order: 5, cutoff frequency ( 3 dB): 0.25
−
−600
0 0.5 1 Normalized Frequency (×π rad/sample)
× f s /2
107
Chebyshev type I filter design example Impulse Response 0.6
1
t r a P 0.5 y r a 0 n i g a −0.5 m I
e d u t i l p m A
−1
−1
0 Real Part
0
) B d ( e −20 d u t i n −40 g a M
−60
0.2 0
−0.2 0
1
10 20 n (samples)
30
0
) s e e r g −200 e d ( e −400 s a h P
0 0.5 1 Normalized Frequency (×π rad/sample)
order: 5, cutoff frequency: 0 .5
0.4
−600
0 0.5 1 Normalized Frequency (×π rad/sample)
× f s /2, pass-band ripple: −3 dB
108
Chebyshev type II filter design example Impulse Response 0.6
1
t r a P 0.5 y r a 0 n i g a −0.5 m I
e d u t i l p m A
−1
−1
0 Real Part
0.2 0
−0.2 0
1
0
0.4
10 20 n (samples)
30
100
) s e e 0 r g e d ( −100 e s a −200 h P
) B d ( e −20 d u t i n −40 g a M
−60
0 0.5 1 Normalized Frequency (×π rad/sample)
order: 5, cutoff frequency: 0 .5
−300
0 0.5 1 Normalized Frequency (×π rad/sample)
× f s /2, stop-band ripple: −20 dB
109
Elliptic filter design example Impulse Response 0.6
1
t r a P 0.5 y r a 0 n i g a −0.5 m I
e d u t i l p m A
−1
−1
0 Real Part
0
) B d ( e −20 d u t i n −40 g a M
−60
0.2 0
−0.2 0
1
10 20 n (samples)
30
0
) s e e −100 r g e d ( −200 e s a −300 h P
0 0.5 1 Normalized Frequency (×π rad/sample)
order: 5, cutoff frequency: 0 .5
0.4
−400
0 0.5 1 Normalized Frequency (×π rad/sample)
× f s /2, pass-band ripple: −3 dB, stop-band ripple: −20 dB
110
Exercise 14 Draw the direct form II block diagrams of the causal infiniteimpulse response filters described by the following z -transforms and write down a formula describing their time-domain impulse responses: 1 (a) H (z ) = 1 12 z −1
− 1 − 41 z −4 ′ (b) H (z ) = 1 − 14 z −1 4
(c) H ′′ (z ) =
1 1 −1 1 −2 + z + z 2 4 2
Exercise 15 (a) Perform the polynomial division of the rational function given in exercise 14 (a) until you have found the coefficient of z −5 in the result. (b) Perform the polynomial division of the rational function given in exercise 14 (b) until you have found the coefficient of z −10 in the result. (c) Use its z -transform to show that the filter in exercise 14 (b) has actually a finite impulse response and draw the corresponding block diagram. 111
Exercise 16 Consider the system h : xn yn with yn + yn−1 = xn xn−4 . (a) Draw the direct form I block diagram of a digital filter that realises h. (b) What is the impulse response of h? (c) What is the step response of h (i.e., h u)? (d) Apply the z -transform to (the impulse response of) h to express it as a rational ratio nal functi function on H (z ). (e) Can you eliminate a common factor from numerator and denominator? What does this mean? (f) For what values z C is H (z ) = 0? (g) How many poles does H have in the complex plane? (h) Write H as a fraction using the position of its poles and zeros and draw their location in relation to the complex unit circle. (i) If h is applied to a sound file with a sampling frequency of 8000 Hz, sine waves of what frequency will be eliminated and sine waves of what frequency will be quadrupled in their amplitude?
{ }→{ }
−
∗
∈
112
Random sequences and noise A discrete random sequence xn is a sequence of numbers
{ }
. . . , x−2 , x−1 , x0 , x1 , x2 , . . . where each value xn is the outcome of a random variable xn in a corresponding sequence of random variables . . . , x−2 , x−1 , x0 , x1 , x2 , . . . Such a collection of random variables is called a random process . Each individual random variable xn is characterized by its probability distribution function P xn (a) = Prob(xn a)
≤
and the entire random process is characterized completely by all joint probability distribution functions P xn1 ,...,xnk (a1 , . . . , ak ) = Prob(xn1
≤ a ∧... ∧x ≤ a ) 1
nk
k
for all possible sets xn1 , . . . , xnk .
{
}
113
Two random variables xn and xm are called independent if P xn ,xm (a, b) = P xn (a) P xm (b)
·
and a random process is called stationary if P xn1 +l ,...,xnk +l (a1 , . . . , ak ) = P xn1 ,...,xnk (a1 , . . . , ak ) for all l, that is, if the probability distributions are time invariant. The derivative pxn (a) = P x′ n (a) is called the probability density function, and helps us to define quantities such as the
→ → → →
expected value (xn ) =
E
apxn (a) da 2
| |
2
E | | a p (a) da variance Var(x ) = E [|x − E (x )| ] = E (|x | ) − |E (x )| correlation Cor(x , x ) = E (x · x∗ ) mean-square mean-squa re value (average power) ( xn ) = n
n
n
m
n
n
2
n
xn
2
n
2
m
Remember that ( ) is linear, that is (ax) = a (x) and (x + y) = (x) + (y). Al Also so,, 2 Var(ax) = a Var(x) and, if x and y are independent, Var(x + y) = Var(x) + Var(y).
E ·
E
E
E
E
E
114
A stationary random process xn can be characterized by its mean value mx = (xn ),
{ } E
its variance
σx2 = ( xn
2
E | − m | ) = γ x
xx x x (0)
(σx is also called standard deviation ), its autocorrelation sequence φxx (k ) = (xn+k x∗n ) and its autocovariance sequence
E
·
2
− m ) · (x − m )∗ ] = φ (k ) − | m | A pair of stationary random processes {x } and {y } can, in addition, [(xn+k γ xx xx (k ) = [(
E
x
n
x
n
xx
x
n
be characterized by its crosscorrelation sequence φxy (k ) = (xn+k yn∗ ) and its crosscovariance sequence γ xy [(xn+k xy (k ) = [(
E
E
·
− m ) · (y − m )∗] = φ x
n
y
xy (k )
− m m∗ x
y
115
Deterministic crosscorrelation sequence For deterministic sequences xn and yn , the crosscorrelation sequence is ∞
{ }
cxy (k) =
{ }
xi+k yi .
i=
−∞
After dividing through the overlapping length of the finite sequences involved, cxy (k) can be used to estimate, from a finite sample of a stationary random sequence, the underlying φxy (k ). MATLAB’s xcorr function does that with option unbiased .
yn−l ), then If xn is similar to yn , but lags l elements behind (xn cxy (l) will be a peak in the crosscorrelation sequence. It is therefore widely calculated to locate shifted versions of a known sequence in another one. The deterministic crosscorrelation sequence is a close cousin of the convolution, with just the second input sequence mirrored:
{ }
{ }
≈
{cxy (n)} = {xn} ∗ {y−n} It can therefore be calculated equally easily via the Fourier transform:
∗ C xy xy (f ) = X (f ) Y (f )
·
Swapping the input sequences mirrors the output sequence: cxy (k ) = cyx ( k ).
−
116
Equivalently, we define the deterministic autocorrelation sequence in the time domain as
∞
cxx (k ) =
xi+k xi .
i=
−∞
which corresponds in the frequency domain to C xx f )) = X (f f )) X ∗ (f f )) = X (f f )) 2 . xx (f
·
|
|
In other words, the Fourier transform C xx f )) of the autocorrelation xx (f sequence cxx (n) of a sequence xn is identical to the squared amplitudes of the Fourier transform, or power spectrum, of xn . This suggests, that the Fourier transform of the autocorrelation sequence of a random process might be a suitable way for defining the power spectrum of that random process.
{
}
{ }
{ }
What can we say about the phase in the Fourier spectrum of a time-invariant random process? 117
Filtered random sequences Let xn be a random sequence from a stationary random process. The output
{ }
yn =
∞
k=
−∞
hk xn−k =
·
∞
k=
−∞
hn−k xk
·
of an LTI applied to it will then be another random sequence, characterized by my = mx
∞
hk
k=
−∞
and φyy (k ) =
∞
i=
−∞
φxx (k i)chh (i),
−
where
φxx (k ) = chh (k ) =
E (x · x∗ ) n+k
∞
i=
n
−∞ hi+k hi . 118
In other words:
{ φ (n)} {y } = {h } ∗ {x } ⇒ Φ (f f )) yy
n
n
n
yy
{c (n)} ∗ {φ (n)} |H (f f ))| · Φ (f f ))
=
hh
xx
2
=
xx
Similarly:
{ φ (n)} {y } = {h } ∗ {x } ⇒ Φ (f f )) yx
n
n
n
yx
= =
{h } ∗ {φ (n)} H (f f )) · Φ (f f )) n
xx
xx
White noise A random sequence xn is a white noise signal, if mx = 0 and
{ }
φxx (k ) = σx2 δk . The power spectrum of a white noise signal is flat: Φxx (f f )) = σx2 . 119
Application example: Where an LTI yn = hn xn can be observed to operate on white noise xn with φxx (k ) = σx2 δk , the crosscorrelation between input and output will reveal the impulse response of the system:
{ } { }∗{ } { }
φyx (k ) = σx2 hk
·
where φyx (k ) = φxy ( k ) = (yn+k x∗n ).
−
E
·
120
DFT averaging
The above diagrams show different types of spectral estimates of a sequence xi = sin(2π j 8/64) + sin(2π j 14.32/64) + ni with φnn (i) = 4δi . Left is a single 64-element DFT of xi (wi (with th rec rectan tangul gular ar window). window). The flat spectrum of white noise is only an expected value. In a single discrete Fourier transform of such a sequence, the significant variance of the noise spectrum becomes visible. It almost drowns the two peaks from sine waves. After cutting xi into 1000 windows of 64 elements each, calculating their DFT, and plotting the average of their absolute values, the centre figure shows an approximation of the expected value of the amplitude spectrum, with a flat noise floor. Taking the absolute value before spectral averaging incoherent rent averag averaging ing , as the phase information is thrown away. is called incohe
×
×
{ }
{ }
121
The rightmost figure was generated from the same set of 1000 windows, but this time the complex values of the DFTs were averaged before the absolu abs olute te value was taken. taken. This is cal called led coherent averaging and, because of the linearity of the DFT, identical to first averaging the 1000 windows and then applying a single DFT and taking its absolute value. The windows start 64 samples apart. Only periodic waveforms with a period that divides 64 are not averaged away. This periodic averaging step suppresses both the noise and the second sine wave.
Periodic averaging If a zero-mean signal xi has a periodic component with period p, the periodic component can be isolated by periodic averaging :
{ }
1 x ¯ i = l im k →∞ 2k + 1
k
xi+ pn
n= k
−
Periodic averaging corresponds in the time domain to convolution with a Dirac comb n δi− pn . In the frequency domain, this means multiplication with a Dirac comb that eliminates all frequencies but multiples of 1/p.
122
Image, video and audio compression Structure of modern audiovisual communication systems:
signal
-
sensor + sampling
-
perceptual coding
-
entropy coding
channel coding
-
?
noise
human senses
display
perceptual decoding
entropy decoding
channel
-
?
channel decoding
123
Audio-visual lossy coding today typically consists of these steps:
→ → → → →
A transducer converts the original stimulus into a voltage.
→
An entropy coder turns that into an apparently-random bit string, whose length approximates the remaining entropy.
This analog signal is then sampled and quantized . The digiti digitizatio zation n par paramete ameters rs (samp (sampling ling freque frequency ncy,, quanti quantizatio zation n level levels) s) are pre preferabl ferablyy chosen generously beyond the ability of human senses or output devices.
The digitized sensor-domain signal is then transformed into a perceptual domain. This step often mimics some of the first neural processing steps in humans.
perceptuall model of what This signal is quantized again, based on a perceptua level of quantization-noise humans can still sense.
The resulting quantized levels may still be highly statistically dependen pend ent. t. A prediction or decorrelation transform exploits this and produces a less dependent symbol sequence of lower entropy.
The first neural processing steps in humans are in effect often a kind of decorrelation transform; our eyes and ears were optimized optimized like any other AV AV communications communications system. system. This allows allows us to use the same transform for decorrelating and transforming into a perceptually relevant domain. 124
Outline of the remaining lectures
→ →
Quick review of entropy coding Transfo ransform rm coding: techni techniques ques for converting sequences sequences of highly highly-dependent symbols into less-dependent lower-entropy sequences.
• • •
→ →
run-length coding decorrelation, decorrela tion, KarhunenKarhunen-Lo` Lo`eve eve transf transform orm (PCA) other orthogonal transforms (especially DCT)
Introduction to some characteristics and limits of human senses
• • •
perceptual scales and sensitivity limits colour vision human hearing limits, critical bands, audio masking
Quantization techniques to remove information that is irrelevant to human senses 125
→
Image and audio coding standards
• • • •
A/µ-law coding (digital telephone network) JPEG MPEG video MPEG audio
Literature
→ → →
D. Salomon: A guide to data compression methods. ISBN 0387952608, 2002. L. Gul Gulick ick,, G. Gescheide Gescheider, r, R. Fris risina: ina: Hea Hearing ring.. ISB ISBN N 01 01950 950430 43073 73,, 1989. H. Schiffman: Sensation and perception. ISBN 0471082082, 1982.
126
