Mathematische Methoden des Visual Computing …...Phase (1) [ ] [ ] Imaginary Part Real Part 1 1 1 1...

DSP-1-Introduction 1

Mathematische Methoden des Visual Computing (Signaltheorie)

VU (1+1h) 186.849

http://ti.tuwien.ac.at/rts/teaching/courses/mmvc

Herbert Grünbacher Institut für Technische Informatik (E182)

[email protected]


Exercises & Exam

• (Home) exercises based on Matlab – Have exercises ready for (oral) exam

• Registration for exam:

Matlab Campus (Win/Linux/Mac)

• INTU (€ 18,00): MATLAB /Simulink including 14 selected tool boxes (DVD)

• Download (€13,90): www.zid.tuwien.ac.at/studentensoftware


Alternative solution public domain http://www.scilab.org

Books • James H. McClellan, Ronald Schafer, Mark Yoder:

Signal Processing First Prentice Hall

• John G. Proakis, Dimitris G. Manolakis Digital Signal Processing Prentice Hall

• Vinay K. Ingle und John G. Proakis Digital Signal Processing Using MATLAB Thomson Learning



(Digital) Signal (Processing)

• Signal A physical quantity varying with time, space, or any other independent variable. E.g.: Temperature, pressure, audio, …


Representation of Signals • Signals appear in different forms and

representations, e.g. speech / acoustic signal Microphone / electric signal Magnetization / magnetic signal Sequence of numbers / digital signal (CD)

cont

inuo

us

disc

rete


Composition of Signals

• Impulse Function Representation in the Time Domain

• Fourier Series, Fourier Transformation Representation in the Frequency Domain a.k.a. spectral representation


Processing / System

• A system is a facility which influences, changes, records, plays back, transmits… signals, e.g. CD player

• Systems operate on signals and generate new signals or signal representations


Digital • Representation by numbers

– Sampling: Quantization in the time domain – Digital: Quantization in the amplitude domain

Time Amplitude

continuous discrete

Analog (T&A cont.) digital signal (T&A discrete)


Discretization

• In Time Sampling (Shannon-) Theorem

• In Amplitude no similar theorem but quality requirements


Representation

• We represent Signals and Systems in mathematical form to – describe and understand Signals und

Systems and – design and implement systems with

requested properties

Content of lecture Tuesday 16:00 – 17:30


16.10. Introduction, Complex numbers (13+30) 23.10. Systems, Time functions (16+20) 30.10. Signals & Spectrum (75) 06.11. Sampling (50) 13.11. FIR filters (90) 20.11. z-Transform (30) 27.11. IIR filters (60)

Complex Numbers

Quadratic Equations

Complex Numbers 2

( ) ( )

( ) ( ) ( )

( )

21,2 1 2

21,2 1 2

21

2

,2 1,2

2 3 0 1 1 3 3; 1

2 1 0 1 1 1 1; 1

2 5 0 1 1 5 1 2

3 1 0

1 1 1 0

1 2

x x x x x

x x x x x

x x

x x

x x x x j

x

x j x

− + =

− − = − =

− − −

− − = = ± + = = −

− + = = ± − = =

− + = = ± − =

±

two real roots

double root

two conjugate complex roots

( )1 2 0j+ =

Complex Numbers 3

Real and Imaginary Part

( , ) Re Im z x y x jy z j z= + = +

Complex Numbers 4

Phasor (Zeiger) ≠ Vector (Vektor)

• Vector: Physical Quantity with Direction e.g. Force, Acceleration, Impulse

• Phasor: Representation of a Complex Number

• Algebra only partially the same (e.g. Addition), but not Multiplication (inner and outer product)

Complex Numbers 5

Magnitude and Angle (Phase) z r ϕ= ∠

compass(z)

ϕ

r

Complex Numbers 6

Angle: Calculation Imagination

• We calculate using radian measure • We visualize using degrees

360 2

Representation: [rad] [°]=45 0.785180

ππα

° =

⋅ =

Complex Numbers 7

Cartesian Polar Representation

( )2 2

cos sin cos sinImaginary Part tan arctan

Real Part

x r y r z r j ry yr x yx x

ϕ ϕ ϕ ϕ

ϕ ϕ

= = = +

= + = = =

Complex Numbers 8

Cartesian Polar

• Addition • Subtraction • Conjugate

• Multiplication • Division • Power • Root • Conjugate

Complex Numbers 9

Attention Phase (1)

[ ][ ]

Imaginary PartReal Part

11

11

sin( /180)

arctan Division by Zero

defined only for 90 90

De

a

gree/Radian:

Real Part 0:

arctarctan 0.7854 45

arctan( ) 0.7854 45

n

rad

rad

α π

ϕ

ϕ

ϕ

ϕ

−−

= =

° ⋅

= ⇒

− ° ≤

= °

°

= =

≤

= ° = −135°

Complex Numbers 10

Attention Phase (2) Phase defined only from 0 2 but physical quantities often have: 2n

ϕ πϕ π ϕ

≤ ≤+ ⋅ =

unwrap(phase)

Complex Numbers 11

Attention Phase (3)

.Runding errors may lead tophase jum

and represent the same phase

Magnit

ps .

If (or very small)ude is Zero Phase has no me !g

,anin

π π

π π→

−

−

Complex Numbers 12

Magnitude is positive!

Complex Numbers 13

Matlab (1) MATLAB knows (only) complex numbers:

3 + 4i oder 3 + 4j

Avoid using i or j as variable: 4+3*i 4+3i i=3; i = 4+3*i 13 but 4+3i 4.00 + 3.00i Restore i as imaginary unit: i = sqrt(-1)

Use notation 4 + 3*1i to avoid troubles.

Complex Numbers 14

Matlab (2)

real(z) Real Part of z real(3-4i) 3 imag(z) Imaginary Part of z imag(3-4i) -4 abs(z) Magnitude of z abs(3-4i) 5 angle(z) Phase of z angle(3-4i) -0.9273 conj(z) Conjugate of z conj(3-4i) 3+4i

angle(z) defined from -180° to 180°

Attention transpose z'

Complex Numbers 15

z' is conjugate transpose z.' is non-conjugate transpose

[3+4i, −1 − 1i, 10 − 10i]' [3+4i, − 1 − 1i, 10 − 10i].'

Complex Numbers 16

Euler (1)

z r ϕ= ∠ cos sinje jϕ ϕ ϕ= +

Complex Numbers 17

Euler (2)

c

2 3 4 5

2 3 4 2 4 3

2

os s n

4 3

i

12! 3! 4! 5!

5 51 1 ( )2! 3! 4! 5! 2! 4! 3! 5!

5cos 1 sin 2! 4! 3! 5!

j

x

j

x x x xe x

e j j j j

ϕ ϕ

ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕϕ ϕ

ϕ ϕ ϕ ϕϕ ϕ ϕ

+

= + + + + + +

= + − − + + + = − + − + − + −

= − + − = − + −

Proof algebraic!

Algebra

Geometry

Complex Numbers 18

Euler (3) cos sincos sin2cos

j

j

j j

e je j

e e

ϕ

ϕ

ϕ ϕ

ϕ ϕ

ϕ ϕ

ϕ

−

−

= +

= −

+ =

cos2

sin2

j j

j j

e e

e ej

ϕ ϕ

ϕ ϕ

ϕ

ϕ

−

−

+=

−=

Complex Numbers 19

1je π = −

Complex Numbers 20

2*exp(1i*45*pi/180) 1.4142 + 1.4142i 3.6056*exp(2.1588i)*10.8167*exp(-0.9828i) 15.00 +36.00i compass(exp((i*30*pi/180)*(0:11)))

0.2

0.4

0.6

0.8

1

30

210

60

240

90

270

120

300

150

330

180 0

MATLAB results always in Cartesian form, i is imaginary unit; input in Euler form possible, i or j as imaginary unit allowed.

abs(3*exp(1i*45*pi/180)) 3 (180/pi)*angle (3*exp(1i*45*pi/180)) 45.00

Complex Numbers 21

Complex Algebra (1)

Addition Subtraction

1 2 1 1 2 2

1 2 1 2

( ) ( ) ( ) ( )z z x jy x jy

x x j y y+ = + + +

= + + +1 2 1 1 2 2

1 2 1 2

( ) ( ) ( ) ( )z z x jy x jy

x x j y y− = + − +

= − + −

Complex Numbers 22

Complex Algebra (2)

Multiplication Division

1 2

1 2

1 2 1 2( )

1 2

j j

j

z z r e r e

r r e

ϕ ϕ

ϕ ϕ+

× = ×

=

1

2

1 2

11 2

2

( )1

2

j

j

j

r ez zr er er

ϕ

ϕ

ϕ ϕ−

÷ =

=

Complex Numbers 23

Complex Algebra (3)

1

2

3

4

5

30

210

60

240

90

270

120

300

150

330

180 0

Conjugate Complex, transpose

1 1 1 1

1 1

( ) ( )z z x jy

x jy

′∗ ∗= = += −

( )1

1

1 1

1

j

j

z r e

re

ϕ

ϕ

∗∗

−

=

=

Complex Numbers 24

Complex Algebra (4)

( )NN j N jNz re r eϕ ϕ= =

Power

Complex Numbers 25

Complex Algebra (4)

2

1 0,1, 2,..., 1

j nNN e

n N

π

== −

Root

2( )

0,1, 2,..., 1

nN NjN j Nz re re

n N

ϕ πϕ += == −

1 0Nz − =

Complex Numbers 26

je ϕ

je ϕ−

cosϕ cosϕ

sinj ϕ

cos sin2cos

j

j j

e je e

ϕ

ϕ ϕ

ϕ ϕ

ϕ −

= +

= +

Complex Numbers 27

0j te ω

0cos tω

0sinj tω

0 0( ) 2t t f tϕ ω π= = ⋅

ϕ

t

90° 180° 270°

Complex Function

Complex Numbers 28

Rotating Phasor

-1 0 1-1

-0.5

0

0.5

1

Realteil

Imag

inär

teil

0 0.5 1-1

-0.5

0

0.5

1

t

Imag

inär

teil

-1 0 10

0.2

0.4

0.6

0.8

1

Realteil

t

0

0

0 0cos sin

cos Re

j t

j t

e t j t

t e

ω

ω

ω ω

ω

= +

=

Complex Numbers 29

2cos j t j tt e eω ωω −= +

DSP-3-Systems 1

Systems

DSP-3-Systems 2

Block diagram

1 1

2 2

m n

y xy x

y x

= ⋅

A

Excitation Response

DSP-3-Systems 3

Tasks : Excitation & System

? : Response System Analysis : Excitation & Response

? : System System Synthesis

: System & Response ? : Excitation Measurement

Cause Effect

DSP-3-Systems 4

Linearity (1) Cause EffectCause Effectk k

→⋅ → ⋅

homogeneous

Cause1 Effect1Cause2 Effect2

Cause1+Cause2 Effect1+Effect2

→→→

Superposition principle

1 2 1 2Cause1 + Cause2 Effect1 + Effect2k k k k⋅ ⋅ → ⋅ ⋅

No interaction between causes.

additive

DSP-3-Systems 5

Mental arithmetic

x 4 1024 ?

1000 x 4 4000 20 x 4 80 4 x 4 16

4096

DSP-3-Systems 6

Linearity (2) 1. »Complicated« signals composed of simple

(basic) signals

2. Compute system responses for all basic signals

3. Sum of all system responses of basic signals is system response of complicated signal

Important basic signals: Sinusoidal Function und Impulse Function

Nonlinear Systems (1)

DSP-3-Systems 7

0 0.5 1 1.5 2 2.5 3 3.5 40

2

4

6

8

10

12

14

16

x y

y = x2

2

2 2

27 20 720 400

27 7297 49

400 49 449

= +

= =

≠ =+ =

Superposition principle violated Nonlinear System

DSP-3-Systems 8

Nonlinear systems (2)

0 0.5 1 1.5 2 2.5 3 3.5 40

2

4

6

8

10

12

14

16

x y

y = x2

[ ]

( )( ) ( ) ( )

( )12cos cos( ) cos co

2 2 2

2

2

2

s

2

cos ( ) 2 cos( ) cos( ) cos ( )1 11 cos 2 2 cos( )cos( ) 1 cos 2

cos( ) cos( )

2 2t t t t

A t A t B t B t

A t B t

A t AB t t B tω ω ω

ω ω

ω

ω ωΩ = Ω− + Ω+

= + Ω + Ω =

= + + ⋅ Ω + + Ω

+ Ω =

DSP-3-Systems 9

0 5 100

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Spektrum Eingangssignal0 5 10

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Spektrum Ausgangssignal

0

ω

Ω

2Ω

2ω

Ω-ω Ω+ω

2 2 2[ cos( )] [ cos( )] [ cos( ) cos( )]A t B t A t B tω ω+ Ω + Ω≠

Spectral representation

DSP-3-Systems 10

( ) ( )

( ) ( ) ( )

( )

( )

222 12

2 2212

2 2 212

2cos2

12

cos

2

2

1 cos2

j t j t

A B

j t j t j t j t

j t j t

t

t e e

e e e e

e e

t

ω ω

ω ω ω ω

ω ω

ω

ω

ω

−

+

− −

−

= +

= + +

= + +

= +

Complex Algebra

DSP-3-Systems 11

( ) ( )( ) ( )2 2

2

2 3 31 12 2 2 2

2 23 3 3 31 1 1 12 2 2 2 2 2 2 2

2

(cos 3cos )

2

j t j t j t j t

A B

j t j t j t j t j t j t j t j t

ABA B

t t e e e e

e e e e e e e e

ω ω

ω ω ω ω

ω − Ω − Ω

− − Ω − Ω Ω − Ω

+ Ω = + + + =

+ + + + + +

( ) ( ) ( ) ( )

[ ]( ) ( )

3 31 12 2 2 2

2 2 22 21 1 1 12 2 2 2

cos2

( ) ( ) ( ) ( )3 3 3 31 1 1 12 2 2 2 2 2 2 2

cos(

212

)

2

cos( )

2

2

3 32 2

2

2 cos 2

2

3 cos( ) cos( )

2

2j t j t

t

j t j t j t j t

t t

e e t

e e e e

t t

ω ω

ω

ω ω ω ω

ω ω

ω

ω ω

−

Ω+ − Ω+ Ω− − Ω−

Ω+ Ω−

+ + = +

+ + + =

= Ω+ + Ω−

+

( ) ( )( ) ( )

2 22 23 232

22 312

2

2

4.5 cos 2

2 2 5

2j t j te e tΩ − Ω+ = + Ω

=+

DSP-3-Systems 12

Linear Systems (3)

• Electric circuits consisting of R, C and L • Amplifiers and Filters • Mass-Spring-Damping-Systems • Series and parallel resonance • Differentiation und Integration • Propagation of electromagnetic and acoustic

waves in isotropic media

Linear systems are described by linear differential or difference equations

DSP-3-Systems 13

Linearity (4)

Linearity only within certain physical limits

Beyond those limits systems become nonlinear (e.g. clipping) or get destroyed

DSP-3-Systems 14

Nonlinear Systems (3)

• Electrical power P = RI2 • Electronic circuits like comparators,

discriminators, … • Clipping, nonlinear amplification • (Magnetic) hysteresis • Multiplication of signals • Digital logic gates • Propagation of electromagnetic and acoustic

waves anisotropic media

Nonlinear systems are described by nonlinear differential or difference equations

DSP-3-Systems 15

LinearTimeInvariant-Systems

Delay n0 System

System Delay n0

x[n] x[n-n0] w[n]

y[n] y[n-n0]

TI if: w[n] = y[n-n0]

(L)TI: input signal delayed by n0 causes output signal delayed by n0

DSP-3-Systems 16

• Linear : Nonlinear

• Time independent : Time dependent

• Dynamic : Non-dynamic

• Causal : Non-causal

• Lumped : Distributed

• Time continuous : Time discrete

• Analog : Digital

• Stable : Non-stable

DSP-4-Time functions 1

Time functions


• It is not possible to represent signals (voice, music, physical values, ...) as (simple) time functions.

• Therefore we use simple »test signals«, which can be represented mathematically and investigate the response of the system to the test signal.

• We compose real signals of those test signals (remember the superposition principle in linear systems) and add (superpose) the responses of the test signals to gain the system response of the real signal.


Elementary Functions

1

0 0( )

1 0t

tt

δ−

<= ≥

0

0

( ) 0 0

( ) 1

t t

t dt

δ

δ∞

−∞

= ≠

=∫

u(t)

o


Impulse “Function“ • The Impuls Function is no true function and

doesn‘t describe a distinct function either. • The Impulse Function is zero, but at t = 0, but

it is not defined at t = 0 either. • The impact of the Impulse Function is the

action on other functions: We are not interested in how the Impulse Function looks like or what it is, we are only interested how the Impulse Function acts on other functions – the signal functions s(t).


1 1"form( ) ( ) ( ) ( )al"t

i i i idt d t tdt

δ δ τ τ δ δ+ +−∞

= =∫

-2 -1 1 2

0.2

0.4

0.6

0.8

1Signal Plot

-2 -1 1 2

0.2

0.4

0.6

0.8

1Signal Plot

-2 -1 1 2

0.5

1

1.5

2Signal Plot

-2 -1 1 2

0.5

1

1.5

2Signal Plot

δ0 δ-1

δ-3 δ-2


Impulse Step

1 0

0 0( ) ( )

1 0

t tt dt

tδ δ τ−

−∞

<= = ≥∫


Convolution (1) 1 2 1 2 1 2( ) ( ) ( ) ( ) ( ) ( )g t g t g g t d g t g dτ τ τ τ τ τ∞ ∞

−∞ −∞

∗ = − = −∫ ∫


Convolution (Folding) (2)

• g1(t), g2(t) g1(τ), g2(τ )

• t is independent Variable hence τ

• g1(-τ) »Folding«

• g1(-τ) delayed - by t g1(t - τ) - multiplied - integrated new value of g1∗ g2

1 2 1 2 1 2( ) ( ) ( ) ( ) ( ) ( )g t g t g g t d g t g dτ τ τ τ τ τ∞ ∞

−∞ −∞

∗ = − = −∫ ∫


Convolution (3) 1 2

1 2 1

1 2

2

( ), ( ) [ ], [ ] kDiscrete:

[ ] [ ] [ ] [ ]k

g t g t

g n g n g

g

k

n

n g

n

k

gτ

∞

=−∞

∗

⇒⇒

= −∑


Convolution with elementary function

0

10

Step:

Impulse: 0 ( ) ( )

1 ( ) ( ) ( )

( )t

i g t

i g t t g t

t g d

δ

δ τ τ+

−−= − ∗ =

= ∗ =

∫

The Impulse Function δ0(t) has the property of a sampling function: Convolving with δ0(t) delivers a sample of g(t) at t.


Mathematical representation of the sampling process

0 00

( ) ( ) ( ) ( ) ( )t

g t g t t g t dδ τ δ τ τ+

−

= ∗ = −∫


Sinusoidal Function

0

00

( ) sin( )12 2

s t A t

fT

ω ϕ

ω π π

= +

= =

0

2 3( ) sin4

s t A tTπ π

=

−


Complex Exponential

0

0

( )0 0

( )0

( ) cos sin

( R) ce os

j t

j t

x t Ae A t jA t

x t Ae A t


Euler ( )( ) cos( ) sin( )oj t

o os t Ae A t jA tω ϕ ω ϕ ω ϕ+= = + + +

-1 0 1-1

-0.5

0

0.5

1

Realteil

Imag

inär

teil

0 0.5 1-1

-0.5

0

0.5

1

tIm

agin

ärte

il

-1 0 10

0.2

0.4

0.6

0.8

1

Realteil

t

( )( ) Re cos( )oj tos t Ae A tω ϕ ω ϕ+= = +

DSP-4-Time functions 15 15

12cos( ) j t j tt e e

Signale

Left and right turning phasors generate (real) cosine function. Positive and negative frequencies!

Inverse Euler

Why Cosine?

• Cosine allows representation of DC component (mean value)


cos(2 ) cos( )A f t A t Aπ = =0


Complex Exponential ( )

( ) ;

( )

o

st j

j t

s t Ke K Ae

s t Ae

s j

ϕ

ϕ

ω

σ ω

+=

= = = +

( )oe sR cst tKe K e tσ ω ϕ= +


Time delay Phase delay

1

0 3

1

0 3

1

0 3

2 5

-1 2

( )x t

( 2)x t

( 1)x t

t

t

t

12

0( ) 3 1 3

0 else

t tx t t t

Delay towards »right« (positive time delay)

1( )x t t

negative sign


0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

Phase = -72°

0( ) cosx t A t

Frequency (Period)

0 00

10.025 40T fT

Maximum of Cosine if argument zero

00 11 0 ; 2 40 0.005 1.2566 72

t t

DSP_5-Signals&Spectrum 1

Signals & Spectrum


• A signal is a physical quantity varying with time, carrying information. E.g. voice, audio, video, temperature, …

• Signals can be changed, stored, transmitted Energy consumption

• Signals appear in different forms and representations, e.g. voice/acoustical signal microphone/electrical signal magnetization/magnetic signal series of numbers/digital signal (CD)


• Signals can be processed (analog and/or digital) by signal processing systems, e.g. cassette recorder, CD-player,…

• Systems have signals as input(s) and generate new signals or signal representations as output(s)


• Signals and systems are represented in mathematical form to describe and analyze them to design and implement systems with specified behavior.

• There is an infinite number of signals. To simplify analysis and design, signals are decomposed in simple elementary (basic) signals.


Composition of »sounds«


ASP

ASP

DSP


Synthesis of periodic signals

20

1 1

( ) cos(2 ) k

k

N Nj f t

k k k o

jk

kk k

k

x t A A f

X A

t X e X e

e

Periodic signals can be composed of sinusoids consisting of the fundamental component f0 = 1/T0 and harmonic components. Harmonic components are integer multiples of the fundamental component.


»Complex« Magnitude Phase in real representation

2 2

2

cos(2 )

kk k k

k

k k

j f t j f tk

j f

k

jk

kt

A f t

A e A e

X

e

e

Complex magnitude in complex exponential


Two-sided spectrum Based on inverse Euler we get:

*2 2

1

( )2 2

k k

Nj f t j f tk k

ok

X Xx t X e e

x(t) generated by right and left turning phasors Two-sided spectrum

1, 21, 2, 21

21

,1 12 2

( ,0), , , , ,1 12 2

...o X f X ffX X X f

Spectral representation

• Signals can be composed of sinusoidal components

• The waveform of sinusoids is known, therefore it is sufficient to represent only magnitude(s) and phase(s) as function of frequency spectral lines



Time domain and frequency domain (spectral) representation

t

s(t)

S(f)


Analysis of periodic signals


1. Odd periodic signals can be represented by sine-functions

01 0 2 0 0

00

1

( ) sin sin 2 ... sin2

sin( )2

n

N

kk

as t b t b t b n t

a b k t


2. We multiply signal s(t) by sinω0t and average over period


We remember from mathematics

2

20

0 für sin sin

für

1sin 1 cos2

2 2

T

o o

T

T T

m nn t m t dt

m n

t dt t dtT

0

sin sin 3 0T

t tdt


2 1sin 1 cos2

21 1

2o o

T T

tT

t dt t d

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2-1

-0.5

0

0.5

1

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2-1

-0.5

0

0.5

1

T


3. All integrals but sin2ω0 t are zero multiplication by sinω0 t and integration „sucks“ spectral component ω0 from s(t)

4. For further spectral components we multiply by sin(kω0 t ), k = 2, 3,…

0

2( )sin

T

k ob s t k t dt

T

Real representation of Fourier series


• We repeat for even signals and get in real-valued representation

0

2( )cos

T

k oa s t k t dt

T

• Signals which are neither even nor odd must be composed of sine and cosine components and we need to compute ak and bk .

19

Complex-valued Fourier series

0 0

00

2

2

00

( )

1 ( )

jk tjk t T

k kk k

j k tT Tk

x t X e X e

X x t e dtT

DSP_5-Signals&Spectrum

20

Σ X0 e j k ω0t

k = −∞

∞


21

Σ Xk e j k ω0t

k = −∞

∞


22

Σ Xk e j k ω0t

k = −∞

∞

022 fTπω π= =


23

0

00

00

( )

1 ( )

j kk

k

t

kj tkT

x t X e

X x t e dtT


Phasor interpretation


X1

X2

X3

Xn

ω0

2ω0

3ω0

nω0


01 0 2 0 0

1 0 2 0

00

0

0

0

( ) cos cos2 ... cos2

+ sin sin 2 ... sin

2 ( )cos = mean valu

2 (

e

)

of

s

( )2

in

n

nT

k

k

as t a t a t a n t

b t b

a s

t b n t

a f t k tdtT

b f t k tdT

t

0

0

0

0 1 0 1 2 0 2 0

2 2

12

( ) cos cos 2 ... cos

, tan

( )

1 ( )

T

n n

kk k k k

k

kjk t

kk

jk tk

t

s t c c t c t c n tbc a b

a

s t X e

aX f t e dt

T

012

0 12

for 0( ) for 0

( ) for 0

T

k k

k k

ka jb k

a jb k

real: sine + cosine

real: cosine + phase

complex


Fourier coefficients Fourier coefficients are often based on period 2π

0

13

2 22

4( ) sin sin

1

...1 3

Ta

y t t t


Scaling of time axes

1

0

13

22

2

4( ) sin sin ...3

f

t t

ay t t t


Spectral analysis

Both signals have the same sound.

Hearing doesn‘t recognize phase!

http://www.jhu.edu/~signals/listen-new/listen-newindex.htm









Hearing • Audio band: 16 Hz − 20 kHz • 16 Hz: vibration

3.500 Hz: fundamental frequency Piccolo Flute • Threshold of hearing 20 µPascal (frequency dependent)

Threshold of pain 130 dB • Hearing properties: resolution1/60 of a whole step

= 3 Hz (medium frequency range)

• Logarithmic magnitude characteristic, ear differentiates approx. 325 sound levels

29 DSP_5-Signals&Spectrum


Spectrum numerical

01 ( ) , 0, 1, 2,...jk t

kT

X s t e dt kT

ω−= = ± ±∫

Ts

T

/ stripssN T T


0

0

1

00

1 ( ) , 0, 1, 2,...

1lim ( ) s

s

jkk

TN

jks sT

n

t

nT

X s t e dt kT

s nT e TT

st nT


0

1

00

1 (m )lis

Nn

nT

kk

jsX s nT e

N

Discrete Fourier Transform (DFT)

Neglecting error of numerical integration leads to


% Fourier coefficients |sinus| with FFT N = 256; t = linspace (0,1,(N+1)); t = t([1:N]); s = abs(sin(pi*t)); subplot (311), plot(t,s) d = fft(s)/N; % complex spectrum dMag = abs(d); dPhase=angle(d); % magnitude, phase M = 10; % Fourier coefficients of M spectral lines d0 = dMag(1); dM=2*dMag(2:M); % single-sided spectrum cMag = [d0,dM]; cPhase=dPhase(1:M); subplot (312), stem((0:(M-1)),cMag), xlabel 'Betragsspektrum' subplot (313), stem((0:(M-1)),cPhase), xlabel 'Phasenspektrum'


2 / 2)T T


2 4 cos cos cos( ) ...1.3 3.5 5.7

0.6366 0.4244cos2 0.0849cos4 0.0364cos6 ..

2 4 6

.

t t ts t

t t t

2 / 2T T

sin is function Cosine tereven ms onlyt


Magnitude Phase Magnitude Phase (FFT) (FFT) (Series) (Series)

0.6366 0 0.6366 0

0.4244 3.1416 0.4244 π

0.0849 3.1416 0.0849 π

0.0364 3.1416 0.0364 π

0.0202 3.1416 0.0202 π

0.0129 3.1416 0.0129 π

0.0089 3.1416 0.0089 π

0.0065 3.1416 0.0065 π

0.0050 -3.1416 0.0050 π ≡ - π

0.0040 3.1416 0.0039 π


Dirichlet conditions • For the series to exist, the coefficients must be finite, which is

the case if f(t) is absolutely integrable over one period (weak Diriclet condition):

( )T

s t dt

• s(t) have only a finite number of maxima and minima in one period, and only a finite number of discontinuities in one period (strong Diriclet condition)

• Any periodic waveform generated in a lab satisfies strong Diriclet conditions.

Discrete periodic signals


Signal periodic if [ ] [ ] for all (Shifting property)

A discrete complex exponential can be periodic with period ,but the frequency must be an integer multiple of 2 /

x n N x n n

NNπ

− =

( )

0 0 0 0 0

0

0

1

20

2

1 2

only N different frequencies 0,1,2, , 1

j n N j n j n j N j n

j N j k

kj nj n N

e e e e e

e e N k

e ek N

ω ω ω ω ω

ω π

πω

ω π

− −

=

− −

= → =

= = → =

=→ = −


Non-periodic signals


Periodic extension

Non-periodic signal converted to periodic signal (to use Fourier series) and let tPeriod ∞.


Fourier Transform

( ) ( )

( ) ( )

1( ) ( )2

j t

j t

f t F

F f t e dt

f t F e d


• The Fourier Transform (non-periodic signals) delivers continuous spectra (spectral density).

• The Fourier Series (periodic signals) delivers

discrete spectra.


/ 2

/ 2 / 2

/ 2

2sin1 2

sin2

(

)

s2

in

2

c

j t j je dt eF ej


Spectrum Impulse function

0 )( 1() ( ) j tF tF t e dt

The impulse contains all frequencies with equal magnitude and is therefore an excellent (theoretical) test signal for systems.


Understanding Fourier

• Dk or rather F(ω) magnitude

• Sum if discrete spectral lines kω0

• Integral if spectral density ω

00 0

0

00

1

0

1

1 ( )

1 ( )

( )

(

( )

( ) ( ) 2

)

T

k

N

k k

j n

jk t jk t

sn

j t

kk

D D

D

dtT

s nT eN

f t f t

f dt

e e

e f tF F

j tde t

FS:

FT:


Discrete Fourier transform (DFT) periodic

discrete [n], discrete (ω or k)

1(2 / )

0

1 1

0 0

1(2 / )

0

0,1,

1[ ] [ ]

22 ;

1[ ] [

]

2,...,[

0,1,2,...,

] [ ]

X

[ ] [

1

1

]k n k n

k

n

k s n

N Nj t j t

k

Nj n

N

n

j

k

N

n

k Nx X

k kf t nTN NT

x n X k e k x

X x n

k eN

n

k e k N

N

n

n

eN

π

ω

π

ω

πω π

− −

−−

−

= =

=

−

=

= = =

= ⇔

=

= −

=

=

=

−

∑ ∑

∑

∑


(2 0 / )

(2 / )

(2 2 / )

(2 ( 1) / )

0,1,2,..., 1

1[ ] [0]

1 + [1]

1 [2]

1 + [ 1]

------------------------

---------------------------[ ]

---

-- -0

-[

n

n

n

j N

j N

j N

nj N N

x X eN

X eN

X eN

X

n n N

N

xkX

N e

π

π

π

π −

= = −

=

+

−

(2 0 / )

(2 1/ )

(2 2 / )

(2 ( 1) / )

] [1]

[2]

[ 1

0,1,2,..

]

1

.,

j N

j

k

k

k

k

N

j N

j N N

ex ex e

x

k

N e

Nπ

π

π

π

−

−

−

− −

= −

+

+

+ −

1

0

N

k

−

=∑

1

0

N

n

−

=∑

1(2 / )

0

1[ ] [ ] nN

j k N

knx X k e

Nπ

−

=

= ∑


2

[ ] , 0,1,2, , 1kj n

Nx n e k Nπ

= = −

Re

0 5-1

0

1

k = 0

0 5-1

0

1

k = 10 5

-1

0

1

k = 20 5

-1

0

1

k = 30 5

-1

0

1

k = 4

0 5-1

0

1

k = 70 5

-1

0

1

k = 60 5

-1

0

1

k = 5

0 5-1

0

1

k = 80 5

-1

0

1

k = 90 5

-1

0

1

k = 100 5

-1

0

1

k = 11

N = 8 Sample of next period

1 ↔ 7 2 ↔ 6 3 ↔ 5

8 ↔ 0 9 ↔ 1 10 ↔ 2 11 ↔ 3


( )

( ) ( )

12 /

01

0

cos si[ ] [ ]

[ ] [ ] cos 2 / sin

n

2 /

jakN

j n N

nN

n

k

k

eX x n e

X k

j

x n n N j n Nk

π α α

π π

−−

=

−

=

=

= −

= +∑

∑

( ) ( )3

0

[ ] [ ] cos 2 / 4 sin 2 / 4

[ ] [0]cos(2 0 / 4) [0]sin(2 0 / 4)

0 0 00 [1]cos(2 1 / 4) [1]sin(2 1 / 4)

[2]cos(2 2 / 4) [2]sin(2 2 / 4) [3]cos(2 3

000 /

0

nX x n n j n

X x jxx jxx jxx

k k kπ π

π ππ ππ ππ

=

= −

= ⋅ ⋅ − ⋅ ⋅ ++ ⋅ ⋅ − ⋅ ⋅ ++ ⋅ ⋅ − ⋅ ⋅ ++ ⋅ ⋅

∑

4) [3]sin(2 3 / 4)[ ] [0]cos(2 0 / 4) [0]sin(2 0 / 4)

[1]cos(2 1 / 4) [1]sin(1 1 1

1 11

2 1 / 4) [2]cos(2 2 / 4) [2]sin(2 2 / 4) [3]cos(2 3 / 4) [3]si

1n(1 3

0

2

jxX x jx

x jxx jxx jx

ππ ππ ππ ππ π

− ⋅ ⋅= ⋅ ⋅ − ⋅ ⋅ ++ ⋅ ⋅ − ⋅ ⋅ ++ ⋅ ⋅ − ⋅ ⋅ ++ ⋅ ⋅ − ⋅ ⋅ / 4)

[ ] [0]cos(2 0 / 4) [0]sin(2 0 / 4) [1]cos(2 1 / 4) [1]sin(2 1 / 4) [2]cos(2 2 / 4) [2]sin(2 2 / 4) [3]cos(2 3 / 4) [3]sin

2

(2 3 / 4)[ ] [0]

2 22

1

c

22 22

s2

3 o

X x jxx jxx jxx jx

X x

π ππ ππ ππ π

= ⋅ ⋅ − ⋅ ⋅ ++ ⋅ ⋅ − ⋅ ⋅ ++ ⋅ ⋅ − ⋅ ⋅ ++ ⋅ ⋅ − ⋅ ⋅= (2 0 / 4) [0]sin(2 0 / 4)

[1]cos(2 1 / 4) [1]sin(2 1 / 4)

3 33 33[2]cos(2 2 / 4) [2]sin(2 2 / 4)

[3]cos(2 3 / 4) [3]sin(2 3 / )3 43

3

jxx jxx jxx jx

π ππ ππ ππ π

⋅ ⋅ − ⋅ ⋅ ++ ⋅ ⋅ − ⋅ ⋅ ++ ⋅ ⋅ − ⋅ ⋅ ++ ⋅ ⋅ − ⋅ ⋅

1

0

N

n

−

=∑

1

0

N

n

−

=∑

1

0

N

n

−

=∑

1

0

N

n

−

=∑


Each X[k] DFT-term is the sum of the point-to-point product of the input sequence x[n] and the complex exponential sequence in the form cos(ϕ) – j sin(ϕ).

The “frequency” [k] depends on

• fs the sampling frequency of the original signal and • N the number of samples.

For a sampling frequency of 480 Hz and a 16 point-DFT the frequency spacing is fs/N = 480/16 = 30 Hz and the frequency component are:

X[0]= 0 Hz, X[1]= 30 Hz, X[2]= 60 Hz, …, X[15]= 450 Hz


Matrix representation of DFT

1(2 / )

1(2 / )

0

1

0

2 /0

usi

0,1,2,..., 1

0,1,2,

[ ] [ ] 0,1,2,..., 1

0,1,2,...

[ ] [ ]

1[ ] [ ]

1[ ] [ ]

ng leads to:

..., 1

j

Nj k N

k

Nk

Nk

Nj n N

nN

N

n

Nn

n

k

k

n n N

n n N

x

X

X k eN

x X k WN

x

W e

n e

X n

k k N

k kWx

π

π

π

−

=

−

−−

=

−

−

=

=

=

=

= −

=

=

= −

= −

=

∑

∑

∑1

0

, 1N

nN

−

=

−∑


( )

1 2 1

2 12 4

1 1

[0] [0] 1 [1] 1 [2] 1 [ 1] 1[1] [0] 1 [1] [2] [ 1]

[2] [0] 1 [1] [2] [ 1]

[ 1] [0] 1 [1] [2] [ 1]__________________________

NN N N

NN N N

N NN N

X x x x x NX x x W x W x N W

X x x W x W x N W

X N x x W x W x N W

−

−

− −

= ⋅ + ⋅ + ⋅ + − ⋅

= ⋅ + ⋅ + ⋅ + − ⋅

= ⋅ + ⋅ + ⋅ + − ⋅

− = ⋅ + ⋅ + ⋅ + − ⋅

0

_____________________1NW =

1

0

[ ] [ ] 0,1,2,..., 1N

Nkn

nk kn NX x W

−

=

= −= ∑


2 /j NNW e π−=

1 2 1

2 4 2( 1)

1 2( 1) 1) ( 1)

1 1 1 1[0] [0]1[1] [1]1[2] [2]

1[ 1] [ 1]

NN N N

NN N N

N N N x NN N N

X xW W WX xW W WX x

W W WX N x N

−

−

− − − −

= ⋅ − −

NW

= NX W x

1 1= N

− −= 1 *N NNx W WW X

0.2

0.4

0.6

0.8

1

30

210

60

240

90

270

120

300

150

330

180 0

N=8


Matlab support for WN function DFTmatrixDarst sig=[1 1 1 1 0 0 0 0]; N=length(sig);

SIG=sig*dftmtx(N); sigrueck=(1/N)*(SIG*conj(dftmtx(N))); figure subplot(1, 2 ,1) stem(sig) title 'Signal' subplot(1, 2, 2) stem(abs(sigrueck)) title 'Signal transf & rücktransf'

0 2 4 60

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Signal0 2 4 6

0

0.5

1

1.5

2

2.5

3

3.5

4

Betragsspektrum0 2 4 6

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Signal => transf => rücktransf


34( ) sin(2 . . ) 0.5cos(2 . . )

2 2( ) [ ] cos

1

sin

2x t t t

X k x n kn j knN N

π π π

π π

= + +

= − ∑

Example:

See Excel-sheet


0 1 2 3 4 5 6 7 8-1.5

-1

-0.5

0

0.5

1

1.5

[ ] [0.3536 0.3536 0.6464 1.0607 0.3536 -1.0607 -1.3536 -0.3536 0.3536]x n =


x[n] "2pi/8" n k=1 cos sin x[n].cos x[n].sin

0,3535 0,785398 0 0,000 1,000 0,000 0,354 0,000

0,3535 0,785398 1 0,785 0,707 0,707 0,250 0,250

0,6464 0,785398 2 1,571 0,000 1,000 0,000 0,646

1,0607 0,785398 3 2,356 -0,707 0,707 -0,750 0,750

0,3535 0,785398 4 3,142 -1,000 0,000 -0,354 0,000

-1,0607 0,785398 5 3,927 -0,707 -0,707 0,750 0,750

-1,3535 0,785398 6 4,712 0,000 -1,000 0,000 1,354

-0,3535 0,785398 7 5,498 0,707 -0,707 -0,250 0,250

0,000 -4,000 4 ∠ -90 Grad


x = 0.3535 0.3535 0.6464 1.0607 0.3535 -1.0607 -1.3535 -0.3535

Matlab: X = fft(x)

Magnitude spectrum

34( ) sin(2 . . ) 0.5cos(21 2. . )x t t tπ π π= + +

1 2

[1] 4 [2] 2 1 0.5X X

A A= =

= =

real2 complex

rA

cA

NX A

X AN

=

=

Attention!

Real- & Imaginary part of DFT


0 1 2 3 4 5 6 7-4

-2

0

2

4

Realteil

0 1 2 3 4 5 6 7-4

-2

0

2

4

Imaginärteil

Matlab


FFT Discrete Fourier transform. FFT(X) is the discrete Fourier transform (DFT) of vector X. For matrices, the FFT operation is applied to each column. For N-D arrays, the FFT operation operates on the first non-singleton dimension.

IFFT Inverse discrete Fourier transform. IFFT(X) is the inverse discrete Fourier transform of X.

• fft(x) represents spectrum from 0 to 2π

• Use Matlab function fftshift(X) to shift zero-frequency component to center of spectrum



• All types of FT are based on the complex exponential function (sequence).

• The complex exponential reaches from -∞ to ∞.

• The synthesis of a (continuous) aperiodic signals requires ∞ frequency components (which cancel out by superposition if signal is zero)

• It is not possible to store ∞ components in a digital computer

• The DFT represents discrete, periodic signals.


(Complex) DFT

0 N-1 0 N/2 N-1

0 N-1 0 N/2

N-1

Time domain Frequency domain

Realteil Realteil

Imaginärteil Imaginärteil

Symmetry of DFT


1 bis 1 2 2N Nk k k= = − >

Frequency components

[ ] *[ ]X k X N k= −

Real part even, Imaginary part odd

Frequency spacing of DFT


After sampling the frequency axes is „lost“

ˆ[ ] ( ) cos( ) cos( )s sx n x nT A nT A nω ϕ ω ϕ= = + = +ˆ snTω =

scontinous

k ffN

= Spacing of spectral lines sfN

No frequency components outside spacing raster!

(Circular) Time-shifting


[ ] [ ]x n X k⇔2

shifted[ ] [ ] [ ]kj m

Nx n m X k X k eπ

−− ⇔ =

Each spectral component is shifted by

Phase shift is proportional to frequency (linear phase shift)

2 kj mNeπ

−

68

[ ] [1.0607 0.3536 -1.0607 -1.3536 -0.3536 0.3536 0.3536 0.3536 0.6464 ]x n =

1 90 , 2 45k k= → − ° = → + °

-3 -2 -1 0 1 2 3 4 5 6 7-1.5

-1

-0.5

0

0.5

1

1.5

0 1 2 3 4 5 6 7-4

-2

0

2

4

Betrag

0 1 2 3 4 5 6 7

-50

0

50

Phase

2 138

1352 23

8

270

90 45

+45 315 45

j

j

e

e

π

π

°

°

⇒ − ° + = + °

⇒ ° + = ° = − °

[ ] [0.3536 0.3536 0.6464 1.0607 0.3536 -1.0607 -1.3536 -0.3536 0.3536]x n


DFT


21

0

[ ] [ ]nN j k

N

nX k x n e

π−

=

= ∑−

21

0

1[ ] [ ]kN j n

N

kx n X e

Nk

π−

=

= ∑

iDFT

Computation of iDFT with DFT


2 2

2

*1 1

* *

0 0*

1*

0

1 1[ ] [ ] [ ]

1[ ] [ ]

N N

N

N Nj kn j kn

k k

Nj kn

k

x n X k e X k eN N

x n X k eN

-1FFT

-1

1/N

1/N

Re [ ]X k

Im [ ]X k

Re [ ]x n

Im [ ]x n

( )*

*

*

1.

2.

3. Scaling by 1 /(4. not required for real sequences)

FFT

N

→X X

X

x


FFT

1/N

1/N

Re [ ]X k

Im [ ]X k

Re [ ]x n

Im [ ]x n

Non-periodic signals (Zero Padding)


DFT based on periodic sequences signal representation periodic in time and frequency domain

0 5 10 150

0.2

0.4

0.6

0.8

1Zeitbereich L= 16 N= 16

0 5 10 150

2

4

6

8

Frequenzbereich

20.5 [1 cos 0 1[ ]

0 sonst

n n Lx n L

π − ≤ ≤ − =


0 5 10 15 20 25 30 350

0.2

0.4

0.6

0.8


0 5 10 15 20 25 30 350

2

4

6

8

Frequenzbereich

L = 16 N = 32


0 20 40 60 80 100 120 1400

0.2

0.4

0.6

0.8


0 20 40 60 80 100 120 1400

2

4

6

8

Frequenzbereich

L = 16 N = 128


iscrete ime

(

[ ]) j n

n

N

x n eX ωω∞

−

=−∞

→ ∞

= ∑

(D T FT)


1DTF

1[ ] 0.5 [ ] ( )

1T

0.5

n

jx n n X

e

-10 -8 -6 -4 -2 0 2 4 6 8 10

0.8

1

1.2

1.4

1.6

1.8

2

Betrag

-10 -8 -6 -4 -2 0 2 4 6 8 10-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

Phase

0 0.5 1 1.5 20

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

DSP_6-FIR 1

Digital Analog

1. Digital signals less sensitive to noise than analog signals (signal state only 0 or1)

2. Coding (compression or security) and processing easier as for analog signals

DSP_6-FIR 2

3. Duplication of digital data is easy, no loss of quality, copying analog data (audio tape, film) leads to reduced quality due to noise

4. Storing digital signals easier and more cost effective than analog signals

5. Digital hardware (chip technology) easier and more cost effective to build than analog HW

DSP_6-FIR 3

7. Analog signal processing still required (aliasing and interpolation)

8. Most signals are analog 9. Sampling theorem and A/D- and conversion

limit digital applications (chip technology) 10. Today: audio signal processing digital, video

signal processing a good portion DSP

DSP_6-FIR 4

A/D conversion

S&H circuit A/D converter

DSP_6-FIR 5

Finite Impulse Response Filter

Discrete-Time System

Input Output

[ ] [ ]y n x n[ ]x n

DSP_6-FIR 6

3 point running average

13

13

13

[0] ( [0] [1] [2])[1] ( [1] [2] [3])

...[ ] ( [ ] [ 1] [ 2])

y x x xy x x x

y n x n x n x n

Filter equation

Output before input: non causal!

n n < -2 -2 -1 0 1 2 3 4 5 6 7 n > 7

x[n] 0 0 0 2 4 6 4 2 0 0 0 0

y[n] 0 2/3 2 4 14/3 4 2 2/3 0 0 0 0

DSP_6-FIR 7

Causality

0 5 10 15 20 250

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Now

Past | Future

13[ ] ( [ ] [ ] [ ])

backward vs foreward avera

1 2

ge

y n x n x n x n

DSP_6-FIR 8

General FIR Filter

0

1 1 13 3 3

[ ] [ ]

Number of filter coefficients: Length of filter 1

e.g.: [ ] [ ] [ 1] [ 2

Order of f

]

ilter M,

M

kk

y n b x n k

L M

y n x n x n x n

DSP_6-FIR 9

Example

DSP_6-FIR 10

• We observe a M-point run-into phase (and run-off phase (not shown)).

• Size of higher frequency component has been reduced Low pass behavior

• Output signal phase-shifted to the right

DSP_6-FIR 11

• FIR-Application (1): Suppress high frequency components, ==> low pass

• FIR-Application (2): Suppress low frequency components, ==> high pass

DSP_6-FIR 12

13[ ] ( [ ] [ 1] [ 2])y n x n x n x n

1 1 13 3 3( [ 1] [ 2] [ 3[ 2] ] ][ )x n x n xx ny nn

Av subtracted from aerage ctual va vlue alue

bk = -1/3 2/3 -1/3

bk = 1/3 1/3 1/3

DSP_6-FIR 13

Unit (Im)Pulse Impulse function, Dirac function

1 0[ ]

0 0n

nn

n … -2 -1 0 1 2 3 δ[n] 0 0 0 1 0 0 0 δ[n-2] 0 0 0 0 0 1 0

-4 -3 -2 -1 0 1 2 3 40

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Einheitsimpuls

DSP_6-FIR 14

Representing signals by impulse functions

[ ] [ ] [ 1] [ 2]2 4 6 [ 3 2 44 ] [ ]x n n n n n n

[ ][ ] [ ]k

x kx n n k

DSP_6-FIR 15

Unit Impulse Response Input: Unit Impulse Output: Impulse Response

0

0

[ ] [ ]

[ ]0,1,2,...,

[ ] [ ]0 el e

[ ]

[ ]s

M

kk

Mn

kk

y n b x n k

x nb n M

y n n kbh n

n

bk ][[

]x nn ]

[[

]y nh n

DSP_6-FIR 16

Example 1

0 1 2 3

1 2

Running Avarage Filter [ ] [ ] [ 1] [ 2]

o

b b by n b n b n b n

Impulse response is sequence of filter coefficients

[ ] 0 for 0 and for h n n n M

Length of impulse response ist finite finite impulse response

0 0.5 1 1.5 20

0.05

0.1

0.15

0.2

0.25

0.3

0.35

n (samples)

Am

plitu

de

Impulse Response

B = [1 1 1]/3 A = 1 impz(B,A)

DSP_6-FIR 17

Delay Filter

0

1 2

Filter coefficients 0,0,1

[ ] [ ] [ 1] [ 2]

[ ] [ ]

k

o

by n b x n b x n b x n

y n x n n

0 0 1 [ 2]x n

[ ] [ 2]y n x n Delay by two time steps

0 0.5 1 1.5 20

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

n (samples)

Am

plitu

de

Impulse Response

impz([0 0 1],1)

DSP_6-FIR 18

Convolution (Faltung)

0

[ ] [ ]k

M

k

y n x n kb

Impulse response Filter coefficients

0

replaced by [ ]

[ ][ ] [ ]M

k

k

y n xh k

b h

k

k

n

DSP_6-FIR 19

x = [1 2 3 -2 2 1] h = [3 -1 2 1]

0 1 2 3-1

0

1

2

3

h[k]0 2 4 6

-2

-1

0

1

2

3

x[k]

-6 -4 -2 0-2

-1

0

1

2

3

x[-k]

folded

DSP_6-FIR 20

-5 0 5 10-202

n = 00 5 10

05

1015

x[n]*h[n] n = 0

-5 0 5 10-202

n = 10 5 10

05

1015

x[n]*h[n] n = 1

-5 0 5 10-202

n = 30 5 10

05

1015

x[n]*h[n] n = 3

-5 0 5 10-202

n = 80 5 10

05

1015

x[n]*h[n] n = 0

h = [3 -1 2 1] x = [1 2 3 -2 2 1]

h = [3 ] * x = [1] 3

h = [3 -1 ] * x = [2 1] 5

h = [3 -1 2 1] * x = [ -2 3 2 1 ] - 4

3 5 9 -4 16 0 1 4 1

DSP_6-FIR 21

Superposition of input pulses Weight Weight x (delayed Impulse Response)

x[0] = 1 y[0] = [1]h[n-0] 3 -1 2 1 0 0 0 0 0 0

x[1] = 2 y[1] = [2]h[n-1] 0 6 - 2 4 2 0 0 0 0 0

x[3] = 3 y[3] = [3]h[n-3] 0 0 9 - 3 6 3 0 0 0 0

x[4] = -2 y[4] = [-2]h[n-4] 0 0 0 - 6 2 - 4 - 2 0 0 0

x[5] = 2 y[5] = [2]h[n-5] 0 0 0 0 6 - 2 4 2 0 0

x[6] = 1 x[6] = [1]h[n-6] 0 0 0 0 0 3 - 1 2 1 0

y[n] = ∑x[n]h[n-k] 3 5 9 - 4 16 0 1 4 1 0

Each input sample generates is own impulse response. The system response is the superposition of weighted und time-delayed impulse responses.

DSP_6-FIR 22

Composition of output signal

h[0] = 3 y[0] = [3]x[n-0] 3 6 9 - 6 6 3 0 0 0 0

h[1] =-1 y[1] = [-1]x[n-1] 0 - 1 - 2 - 3 2 - 2 - 1 0 0 0

h[3] = 2 y[3] = [2]x[n-3] 0 0 2 4 6 - 4 4 2 0 0

h[4] = 1 y[4] = [3]x[n-4] 0 0 0 1 2 3 - 2 2 1 0

y[n] = ∑h[k]x[n-k] 3 5 9 - 4 16 0 1 4 1 0

Output signal is sum of weighted and time-delayed input signals. Input signal weighted by filter coefficients.

DSP_6-FIR 23

Convolution and Matlab x = [1 2 3 -2 2 1] h = [3 -1 2 1] conv(h,x) 3 5 9 -4 16 0 1 4 1

1 2 3 4 5 6 7 8 9-4

-2

0

2

4

6

8

10

12

14

16

stem(conv(h,x))

DSP_6-FIR 24

Willkommen!

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5x 104

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0 0.5 1 1.5 2 2.5x 104

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

DSP_6-FIR 25

DSP_6-FIR 26

Remember: Multiplication of polynomials

2 21 2 1 2

21 0 2 0

2 31 1 1 2 1

2 3 42 1 2 2 2

( )( )

o o

o o

o

o

a a x a x b b x b x

a b a b x a b x

a b x a b x a b x

a b x a b x a b x

A = [a0,a1,a2] B = [bo,b1,b2] conv(A,B)

DSP_6-FIR 27

Implementation of FIR-Filters

0

[ ] [ ]M

kk

y n b x n k

Addition Multiplication Delay (Storage)

x[n] ×

y[n]

b

+

x1[n]

y[n] x2[n]

x[n] y[n] Unit Delay

DSP_6-FIR 28

Block diagram (direct form) 3

0

0 1 2 3

[ ] [ ]

[ ] [ ] [ 1] [ 2] [ 3]

kk

y n b x n k

y n b x n b x n b x n b x n

×

Unit Delay

Unit Delay

× × ×

+ + +

0b 1b 2b3b

[ ]x n [ 1]x n [ 2]x n [ 3]x n

[ ]y n

Unit Delay

DSP_6-FIR 29

Block diagram (transposed form)

×

Unit Delay

Unit Delay

Unit Delay

× × ×

+ + +

3b 2b 1b0b

[ ]x n

[ ]y n

1[ ]v n2[ ]v n3[ ]v n

0 1

1 1 2

2 2 3

3 3

[ ] [ ] [ 1][ ] [ ] [ 1][ ] [ ] [ 1][ ] [ ]

y n b x n v nv n b x n v nv n b x n v nv n b x n

2 2 3

1 1 2 3

0 1 2 3

[ ] [ ] [ 1][ ] [ ] [ 1] [ 2][ ] [ ] [ 1] [ 2] [ 3]

v n b x n b x nv n b x n b x n b x ny n b x n b x n b x n b x n

DSP_6-FIR 30

Convolution and LTI systems [ ][ ]

x nn

[ ][ ]

y nh nLTI System

[ ] [ ] [ ] general case l

x n x l n l l

Input signal is a sequence of weighted and time-delayed impulses. The system response is:

[0] [0] [0] [ ][1] [ 1] [1] [ 1][2] [ 2] [2] [ 2][ ] [ ] [ ] [ ]

x x h nx n x h nx n x h nx l n l x l h n l

DSP_6-FIR 31

Convolution and LTI Systems (2)

Application of superposition principle leads to:

[ ] [ ] [ ] [ ] [ ] [ ]l l

x n x l n l y n x l h n l Most general case: Neither x[n] nor h[n] are of finite duration:

[ ] [ ] [ ]l

y n x l h n l

A LTI-System can be represented by its convolution sum.

DSP_6-FIR 32

Properties of LTI systems

0

1 2

[ ] [ ] [ ] [ ] [ ]

Convolution operator

Convolution with an impulse [ ] [ ] [

Commutativ:

]

[ ] [ ] [ ] [ ]

[ ] [ ] [ ] [ ] [ ]

[ ] [ ]Associative:

l

o

l

y n x n h n x l h n l

x n n n x n n

x n h n h n x n

y n h n x n h l x n l

x n x n x

3 1 2 3[ ] [ ] [ ] [ ]n x n x n x n

DSP_6-FIR 33

Cascaded LTI Systems LTI 1 LTI 2

LTI

LTI 2 LTI 1

1 2 1 2

2 1 2 1

[ ] ( [ ] [ ]) [ ] [ ] [ ] [ ]

[ ] [ ] [ ] ( [ ] [ ]) [ ]

y n x n h n h n x n h n h n

x n h n h n x n h n h n

1[ ] [ ]x n h n[ ]x n 1 2( [ ] [ ]) [ ]x n h n h n

2[ ] [ ]x n h n[ ]x n 2 1( [ ] [ ]) [ ]x n h n h n

1 2[ ] [ ]h n h n

1[ ]h n 2[ ]h n

1[ ]h n2[ ]h n

Frequency Response of FIR Filters (Frequenzgang)

DSP_6-FIR 35

Sinusoidal input signals

ˆ is a dimension-less

ˆ[ ] ( ) cos( ) cos( )ˆ ˆ is the normalized frequenc

uniy

t.

s s

s

x n x nT A nT A nT

The index in x[n] is also dimensionless. A sampled (continous) signal x(t)carries no time information : A time - discrete signal is only a sequence ofnumbers and doesn't contain information

To reconstruct the signal the s

about the sampl

ampling period m

i

u

ng

st

int

be

erva

kn

ll.

own!

DSP_6-FIR 36

ˆ

ˆ (

0

)

0

0

( ) sampling continous discrete

ˆˆ /2

[ ]

[ ]

[ ] []

]

[

M M

j j t

s s

j j n

jM

jk k

k

kn

kk k

j

b k

x t Ae e

f fT

x n Ae e

x n x n

Ae e

f f

y n

b e

h k k

b

ˆˆ

0

( )ˆ

j jM

k

k

nAe e

H

DSP_6-FIR 37

The complex exponential is the only signal passing the linear system without changing the wave form. Only amplitude and phase (complex amplitude) change, but the sinusoidal property remains. The expression is only valid for the input signal and is meaningless for any other signals.

ˆ[ ] [ ]( )y n x n= Hˆj j nAe e

DSP_6-FIR 38

Frequency response ˆ

ˆ ˆ

0 0

ˆ( )

ˆ( )

ˆ ˆ( ) describes response of the LTI-

[ ]

Sys

[

tems for ev

]

ery .

j j n

M Mj k j k

k kk k

y n Ae e n

b e h k e

= H

H

Hˆ ˆj n j ne e A B

same frequency

but different (complex) amplitudes

DSP_6-FIR 39

Example

ˆ ˆ 2

ˆ

ˆ

ˆ ˆ

1 2 1

ˆ( ) 1 2 ( 2 ) ˆ2 2cos

ˆ

(ˆ

) ˆ) (

j

k

j j

j j j

b

e ee

ee e

H

H

ˆ

0

ˆ( )M

j kk

kb e

H

DSP_6-FIR 40

Frequency response (2)

ˆ( )

ˆ( ) ... complex funktionˆ ˆ ˆ( ) ( ) ( )

ˆ ( ) j

e j m

e

H

H

H H

H

H

Magnitude, Betrag Gain, Verstärkung

Phase

DSP_6-FIR 41

Bode diagram

0.5 ½

fsampling π

DSP_6-FIR 42

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9-180

-160

-140

-120

-100

-80

-60

-40

-20

0

Normalized Frequency (×π rad/sample)

Phas

e (d

egre

es)

Phase Response

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.90

0.5

1

1.5

2

2.5

3

3.5

4


Mag

nitu

de

Magnitude Response

Low pass bk = [1 2 1]

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-200

-150

-100

-50

0


Pha

se (d

egre

es)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-100

-80

-60

-40

-20

0

20


Mag

nitu

de (d

B)

Magnitude in logarithmic representation

freqz([1 2 1],1)

fvtool([1 2 1],1)

1

Nomalized Frequency (∙ π rad/sample) 1

1

1

DSP_6-FIR 43

Logarithmic representation

log [Bel]out

in

P aP

2

2log 2log [Bel] 20log [dB]out out out

in in in

U U UU U U

To describe properties of filters one is interested in the ratio of input power to output power. Because this power ratio spans several powers of ten logarithmic representation is appropriate. The ratio is measured in Bel.

Comparing (voltages) amplitudes leads to:

10 20 dB

1/1000 -60 dB 1/100 -40 dB 1/10 -20 dB

-3 dB 1 0 dB

3 dB 2

1/ 2

DSP_6-FIR 44

Example FIR Filter

ˆ

ˆ ˆ

ˆ ˆ ˆ2 3 4

ˆ ˆ2 2

ˆ2

ˆ2

ˆ2

2 ˆˆ ˆ2 ˆ

1, 2,4, 2,1

ˆ( ) 1 2 4 2after separating we get

ˆ(

ˆ2c

2 2

ˆ2 2 4cosˆ4cos

) 4

ˆ( )ˆ( ) 4

os2 ˆ2cos

4

42

j j

j j

j j

j j

k

j j j j

j

j

j

ee e

b

H e e e ee

H e

H e

e

e ee

H

e

ˆ ˆ ˆ ˆcos 2cos2 ( ) 2

DSP_6-FIR 45

bk = [1 -2 4 -2 1]

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.451

2.8

4.6

6.4

8.2

10

Frequency (Hz)

Mag

nitu

deMagnitude and Phase Responses

-7

-5.6

-4.2

-2.8

-1.4

0

Phas

e (r

adia

ns)

DSP_6-FIR 46

In example [1 2 1] and [1 -2 4 -2 1] the phase ˆ ˆ ˆ ˆ( ) ( ) 2

is a linear function of

FIR-Filters have linear phase if the filter coefficients are:

[1 -2 4 -2 1]

für 0,1,...,k M kb b k M

0

*0, 0 *

0, 0

If is a zero of H( ) is, than1 1 , ,

are also zeros of ( ).

z z

z zz z

H z

DSP_6-FIR 47

Meaning of linear phase

0 0.5 1 1.5 2 2.5 3 3.5-1

-0.5

0

0.5

1

0 0.5 1 1.5 2 2.5 3 3.5-1

-0.5

0

0.5

1

0 0.5 1 1.5 2 2.5 3 3.5-1

-0.5

0

0.5

1

0 0.5 1 1.5 2 2.5 3 3.5-1.5

-1

-0.5

0

0.5

1

1.5

0 0.5 1 1.5 2 2.5 3 3.5-1

-0.5

0

0.5

1

0 0.5 1 1.5 2 2.5 3 3.5-1

-0.5

0

0.5

1

Phase shift constant ϕ1 = ϕ 2 = ϕ3 = 60°

Phase linear (proportional to frequency) ϕ1 = 30 ° ϕ2 = 60° ϕ3 = 90°

Linear phase maintains wave form of signal

DSP_6-FIR 48

Sounds are phase independent: All 4 wave forms sound the same

0 20 40 60 80 100 120 140 160 180 200-3

-2

-1

0

1

2

3

0 20 40 60 80 100 120 140 160 180 200-3

-2

-1

0

1

2

3

0 20 40 60 80 100 120 140 160 180 200-3

-2

-1

0

1

2

3

0 20 40 60 80 100 120 140 160 180 200-3

-2

-1

0

1

2

3

Test signal for linear phase

DSP_6-FIR 49

DSP_6-FIR 50

0 200 400 600 8000

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Tscheb 5. Ordnung, 1 dB

200 400 600 800

0.5

1

1.5

0 200 400 600 8000

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Eingangssignal

200 400 600 800

0.5

1

1.5

0 200 400 600 8000

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

LinPhase 40.Ordnung, max. flach

200 400 600 800

0.5

1

1.5

Ringing

DSP_6-FIR 51

Superposition (Überlagerung)

*

01

01

*

0

0

ˆ ˆ

ˆ ˆ

1

[ ]2 2

cos

... (for real sequences [ ]ˆ ˆ( ) ( )

ˆ ˆ(0

)

[ ]2 2

) ( ) ( )

( 0)

ˆ

k k

k k

Nk k

k

N

k kk

Nk k

k

j n j n

k

j n j nk k

X Xx n X

X X X

x n

X Xy n X

e

n

X

e

e e

*H H

H H H

H 1

ˆˆ ˆ( ) (cos )N

k kk kk

k X n X

H H

DSP_6-FIR 52

Example

4

910

1,2,1

[ ]

(0)

4

3.414

94 3cos 2cos4 2 10

9

4

0.097

4 3.414 0.03cos 2cos

49

4

10

[ ]2

94

7

k

j

j

b

x n

e

e

y n

n n

n

H

H

H

00 119n

DSP_6-FIR 53

Using Matlab

[H,w] = freqz(B,A,w) [H,w] = freqz( [1 2 1],1, [0,(pi/4),(pi*9/10)]) [4,(2.4142-2.4142i),(-0.0931-0.0302i] abs(H) und angle(H) [4ej0,3.41420e-j0.7854,0.0979e-j2.8274] [4,3.41420e-jπ/4,0.0979e-j9π /10]

DSP_6-FIR 54

FREQZ Digital filter frequency response. H = FREQZ(B,A,W) returns the frequency response at frequencies designated in vector W, in radians/sample (normally between 0 and pi) of the filter:

ˆ ˆ0 1

ˆ ˆ0 1

....ˆ( )....

Numerator and denominator coefficients given in vectors B and A.

j j mm

j j nn

b b e b eHa a e a e

A = [1] for FIR-Filters

DSP_6-FIR 55

Frequency Domain

ˆ( )HSpectrum

of signals is modified by the

frequency response of the system

Single sampli NOTng values considered.

DSP_6-FIR 56

Properties of frequency response

ˆ

0 0

ˆ 2

0 1

ˆ

ˆ ˆ( ) (

ˆ[ ] [ ] [ ] ( ) [ ]

ˆ ( 2 )

ˆ ˆ( ) ( )

) wenn reel

l

M Mj k

k kM

j k j kk

k

k

h n h k n k h k e

b e e

b

*

H H

H

H

H H

H(ω) Spectrum periodic

ˆ( )

ˆe ( ) ... ˆ ˆ ˆ( ) ( ) ( )

ˆe ( )ˆm ( )

ˆFrequency respose only needed for 0 .

ˆm ( ) . .

..

H

H

H

H

H

H

Magnitude ... evenPhase ... odd

Re evenIm odd

. Symmetry reasons

DSP_6-FIR 57

Cascaded LTI Systeme

1 2

ˆ ˆ ˆ ˆ ˆ( ) ( ) ( ) ( ) ( )ˆ ˆ[ ] [ ] ( ) ( )h n h n

1 2 2 1

1 2

H H H H H

H H

Convolution in time domain Multiplication in frequency domain

DSP_6-FIR 58

ˆ ˆ ˆ ˆ2 3 4

ˆ ˆ ˆ2 3

Cascaded systemsˆ( ) 2 4 6 4 2ˆ( ) 1 2 2

ˆ ˆ ˆ( ) ( ) ( )

j j j j

j j j

e e e ee e e

1

2

1 2

H

H

H H H

w=-pi:(pi/500):pi; b1=[2 4 6 4 2]; b2=[1 -2 2 -1]; b=conv(b1, b2); H1=freqz(b1,1,w); H2=freqz(b2,1,w); H=freqz(b,1,w);

-4 -3 -2 -1 0 1 2 3 40

5

10

15

20

-4 -3 -2 -1 0 1 2 3 40

2

4

6

-4 -3 -2 -1 0 1 2 3 40

5

10

15

H1

H2 H3

DSP_6-FIR 59

First Difference System bk =[1 -1]

ˆ

[ ] [ ] [ 1]ˆ ˆ ˆ[ ] 1 1 cos sinˆ ˆ[ ] 2 sin( / 2)

ˆsinˆ[ ] arctanˆ1 cos

j

y n x n x ne j

H

H

H

DSP_6-FIR 60

High pass bk =[1 -1]

-15 -10 -5 0 5 10 150

0.5

1

1.5

2

ω

Bet

rag

-15 -10 -5 0 5 10 15-2

-1

0

1

2

ω

Pha

se

2π 4ππ

DSP_6-FIR 61

Running-Average Filter

1

01

ˆ

0

1

0

ˆ ˆ ˆ/2 /2 /2ˆ1ˆ

ˆ ˆ ˆ ˆ/2 /2 /20

1[ ] [ ]

( )

L-point running averager

1( )

1Using 1

1 1 1 11

L

kL

j k

kLL

k

k

j L j L j Lj LLj k

j j j jk

y n x n kL

eL

e e eeeL L e L e e e

H

H

Phase

Magnitude

ˆ 1 /2

ˆsi

n / 2

ˆs n

/ 2=

i

LjLL

e

DSP_6-FIR 62

ˆ ( 1)/2ˆ( )ˆsin / 2

ˆsin /

ˆ( )

ˆ(

ˆ(

2)

)

j

L

L

LLD

Le

L

H

D

D

diric

is called Diriclet - or periodic sinc -Function,the corresponding Matlab function is .

DSP_6-FIR 63

x=-pi:(pi/200):pi; plot (x,diric(x,11));

Periodic sinc-Function from -3π to 3π -4 -3 -2 -1 0 1 2 3 4

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

-10 -8 -6 -4 -2 0 2 4 6 8 10-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

π -π

DSP_6-FIR 64

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.90

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1


Mag

nitu

de

Magnitude Response

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9-200

-150

-100

-50

0

50

100


Phas

e (d

egre

es)

Phase Response

X: 0.422 Y: -20.2

X: 0.418 Y: -16.6

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9-900

-800

-700

-600

-500

-400

-300

-200

-100

0


Con

tinuo

us P

hase

(deg

rees

)

Continuous Phase Response

Phase ±180° Phase continuous

Run

ning

Avg

.-Filt

er b

k =

ones

(1,1

1)/1

1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9-90

-80

-70

-60

-50

-40

-30

-20

-10

0


Mag

nitu

de (d

B)

Magnitude Response (dB)

Amplitude linear Amplitude in dB

DSP_6-FIR 65

Transient behavior (1) • Frequency domain analysis only for

sinusoidal signals, i.e. input (output) signal reaches from -∞ to ∞.

• But real signals are switched on/off. Steady state of FIR filter only after all delay elements are »filled« with values from input signal.

• Until this state is reached analysis of the system behavior must be done in the time domain (transient analysis).

DSP_6-FIR 66

Transient behavior (2) bk = [1,2,4,2,1]

DSP_6-FIR 67

Low pass – high pass

DSP_6-FIR 68

Band pass – band stop

DSP_6-FIR 69

Multiband

DSP_6-FIR 70

From ideal to real low pass filter

ˆ1ˆ( )

ˆ0c

iLPc

Hω ω

ωω ω π

≤= < ≤

π ωc

ˆ1 1 sinˆ ˆ( ) , 2 2

c cj n j nj n c

iLP iTPe e nh H e d n

jn jn n

π ω ωω

π

ωω ωπ π π

−

−

= = − = −∞ < < ∞

∫

Impulse response not causal!

Fourier transform

DSP_6-FIR 71

-10 -5 0 5 10-0.2

0

0.2

0.4

0.6

0.8

1

1.2

Time (sec)

Am

plitu

deSinc Function

-M M

DSP_6-FIR 72

( )( )

sin 0 1, 2 1

0 else

c n Mn M

LP

n N N Mh

ωπ

− −

≤ ≤ − = +=

Truncating at – M and M and shifting to the right by M produces a (different but) causal low pass filter

Computing frequency response of filter hLP (again by using the Fourier transformation)

DSP_6-FIR 73

Above expression corresponds to multiplication with rectangular window

LP idealLP Rh h w= ×

( )( )

sin 0 1, 2 1

0 else

c n Mn M

LP

n N N Mh

ωπ

− −

≤ ≤ − = +=

DSP_6-FIR 74

LP idealLP Rh h w= ×

Frequency domain

idealLPh ⇒ Ideal filter

ˆ(2 1)2ˆˆ2

1 0[ ]

0 else

sinˆ( )

sin

R

MMj n

Rn M

n Mw n

W eω

ωω

ω+

−

=−

≤ ≤=

= =∑

DSP_6-FIR 75

ˆ( )1ˆ ˆ ˆ( ) ( ) ( ) ( ) ( )2

j jLP idealLP R idealLP RH H W H e W e d

πφ ω φ

π

ω ω ω φπ

−

−

= ∗ = ∫

LP idealLP Rh h w= ×Multiplication in the time domain ↓ convolution in the frequency domain

DSP_6-FIR 76

DSP_6-FIR 77

Gibbs‘ phenomenon

DSP_6-FIR 78

Windowing (1)

Time

Frequency

DSP_6-FIR 79

Windowing (2)

Hamming2[ ] 0.54 0.46cos

2 1nw n M n M

Mπ = + − ≤ ≤ +

Time

Frequency

DSP_6-FIR 80

( )( )

20

Kaiser0

1 /,

I n Mw M n M

I

β

β

− = − ≤ ≤

0 50 1000

0.2

0.4

0.6

0.8

1

β = 00 50 100

0

0.2

0.4

0.6

0.8

1

β = 2

0 50 1000

0.2

0.4

0.6

0.8

1

β = 40 50 100

0

0.2

0.4

0.6

0.8

1

β = 6

DSP_6-FIR 81

Properties of windows

Window B A Rectangular 4π/(2M+1) -13 dB Hann 8π/(2M+1) -32 dB Hamming 8π/(2M+1) -43 dB Blackmann 12π/(2M+1) -58 dB

B

A

Sharper edge from pass to stop band Lower attenuation A in stop band

DSP_6-FIR 82

10 20 30 40 50 600

0.2

0.4

0.6

0.8

1

Samples

Ampl

itude

Time domain

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9-140

-120

-100

-80

-60

-40

-20

0

20

40


Mag

nitu

de (d

B)

Frequency domain

Rechteck

von Hann

Dreieck

Hamming

Comparing windows in Matlab with wintool

DSP_6-FIR 83

Matlab function FIR1

B = FIR1(N,Wn) designs an N'th order lowpass FIR digital filter and returns the filter coefficients in length N+1 vector B. The cut-off frequency Wn must be between 0 < Wn < 1.0, with 1.0 corresponding to half the sample rate. The filter B is real and has linear phase. The normalized gain of the filter at Wn is -6 dB.

DSP_6-FIR 84

B = FIR1(N,Wn,'high') designs an N'th order highpass filter. You can also use B = FIR1(N,Wn,'low') to design a lowpass filter. If Wn is a two-element vector, Wn = [W1 W2], FIR1 returns an order N bandpass filter with passband W1 < W < W2. You can also specify B = FIR1(N,Wn,'bandpass'). If Wn = [W1 W2], B = FIR1(N,Wn,'stop') will design a bandstop filter.

DSP_6-FIR 85

B = FIR1(N,Wn,WIN) designs an N-th order FIR filter using the N+1 length vector WIN to window the impulse response. If empty or omitted, FIR1 uses a Hamming window of length N+1. For a complete list of available windows, see the help for the WINDOW function.

DSP_6-FIR 86

% Lowpass filter

% Ordnung 20

%

b=fir1(20,0.2);

disp (b);

[h,w]=freqz(b,1,1024);

mag=20*log10(abs(h));

plot(w/pi,mag);grid 0 0.2 0.4 0.6 0.8 1-120

-100

-80

-60

-40

-20

0

fir (20, 0.2)

DSP_6-FIR 87

Practical filter design

Design software allows specification of pass and stop bands (tolerance diagrams). Filter characteristic are very sensitive to location of poles and zeros.

DSP_6-FIR 88

Matlab fdatool

DSP_6-FIR 89

Summary(1) • FIR filters use only values from input sequence x[n]. • Filter coefficients b[k] determine behavior of filter (low pass, high pass, band pass, band stop). • Impulse response and frequency response to describe properties of (FIR-) filters. • Impulse response (and transient response) has finite duration.

DSP_6-FIR 90

Summary (2)

• FIR filters may have (and usually have) linear phase, a property important for maintaining the signal wave form. • FIR filters are always stable and easy to realize in hardware and software. • FIR filters of the same order as IIR filters are less selective. In other words: To realize a given selectivity using FIR filters a higher filter order (higher cost) is required.

z-Transform

z-Transform 2

Why z-Transform? The z-Transform introduces polynomials and rational functions in the analysis of linear time-discrete systems and has a similar importance as the Laplace transform for continuous systems Convolution becomes a multiplication of polynomials Algebraic operations like division, multiplication and factoring correspond to composing and decomposing of LTI systems Carrying out the z-transform in general leads to functions consisting of nominator and denominator polynomials The location of the roots of these polynomial determins important properties of digital filters

z-Transform 3

ˆ, , Domainn z

ˆ Frequency dom

Impulse response, difference equations, "real" signal domain

Frequency response, spectral representation, analysis of sou

a

nd

Opera Doma

Tim

in

e do

ton

ma

ir

in

z z

n

s, poles and zeros,

mathematical analysis and synthesis

z-Transform 4

Input signals

Impulse Complex exponential function

• x[n] = zn

ˆj nAe

z-Transform 5

z-Transform and FIR- filter

0

0

[ ] [ ] (using filter coefficients)

or notation using convolut

for a

ion sum

[ ] [ ] [ ] [ ] [ ]

We now use the "signal" [ ] ll n ,

M

kk

M

kk

n

y n b x n k

y n x n h n h n b n k

x n z z

0

0

0 0 0

0

arbitrary (complex) numberas input signal

[ ] [ ]

...[ ] system fu( ) nction

M M Mn k n k

Mk

kk

M

nk k k

k

Mk

k

kk

k

k

k

y n b x n k b z b z z b z

b z h k z

z

H z

ˆ

Remember:jz e

z-Transform 6

0 0

[ ] [ ] (

T

)

For the input signal we get[ ] [ ] ( )( ) is the Transform

he system function ( ) is the Transformof the impulse response:

of [ ].

M Mk

k kk k

n

n n

h n b n k H z b z

zy n h n z H z zH z z h n

H z z

z-Transform 7

0 1 2 1 20 1 2

2

2 2

[ ] ( )1,2,1 ( ) 1 2 1

2 1 2 1 =1+z

k

h n H zb H z b z b z b z z z

z zz z

The system function is a rational function numerator polynomial

denominator polynomial

2 2 1z z

2z

z-Transform 8

0 1 2

21 2 0 1 2

0 1 2 2

Impulse response: [ ] [ ] [ 1] [ 2]

System function: ( )

h n b n b n b n

b z b z bH z b b z b zz

z - Transform

z-Transform 9

Representation of signals

0

A signal of finite lenght can be represented as

[ ] [ ] [ ]

the Transform of this signal is

is an arbitrary complex number, e.g., is the independent(complex) variable of

0( ) [ ]

N

k

N

x n x k n k

z

k

kz z

X z x k z

1

01

the -Transform.Alternatively we may write

( ) [ ]( )

which means, that X( ) is polynomial of order N of the variable .

Nk

k

z

X z x k z

z z

0 0

... System function( ) [ ] M M

k kk

k kH h kz b z z

z-Transform 10

Transform The transition from n z is called z-Transform of x[n].

X [z] is a polynomial in z-1 , the coefficients are the values of the sequence x[n]. e.g.:

n n < -1 -1 0 1 2 3 4 5 n > 5x[n] 0 0 2 4 6 4 2 0 0

1 2 3 4( ) 2 4 6 4 2X z z z z z

z-Transform 11

0 0

[ ] [ ] [ ] ( ) [ ]

[ ] ( )

N Nk

k k

n Domain z Domain

x x k kn n X

zXn

x kz

x

z

0

Example:

[ ] [ ]

[ ]o

n

x n n n

X z z

0

[ ]

of [ ]

j

Nj k j k

kk

z e

z re x k r e

x k r

ω

ω ω− −

=

−

= →

= →∑

DFT

DFT

z-Transform 12

The Transform of a FIR-Filters is apo

The system function ( ) is a function

lynom of degree and ha

of the comp

s zeros(fundame

lex

nta

variable

l theore

m of algebra).

. z

M M

H zz

1 11 2 1 1 3 2

2

1 13 2

Example:[ ] 6 [ ] 5 [ 1] [ 2]

( ) 6 5 (3 )(2 ) 6

Zeros at and

y n x n x n x nz z

H z z z z zz

z-Transform 13

Properties of the z-Transform

1 2 1 2

1

[ ] [ ] ( ) ( )

in the domain correponds to a time de

Linea

layof 1 in the domain (time do

rity

Time

mai

-delay

n)

ax n bx n aX z bX z

z zn

n n < 0 0 1 2 3 4 5 n > 5 x[n] 0 3 1 4 1 5 9 0

1 2 3 4 5

1

0 1 2 3 4 5 6

( ) 3 4 5 9( ) ( )

= (0 ) 3 4 5 9

X z z z z z zY z z X z

z z z z z z z

z-Transform 14

z-Transform as Operator

1

1

1

1 1

The ex

Unit-Delay

pression

Operator [ ] [ ] [ 1][ ] for all [ ] [ ] [ ]

Nervertheless it is common to use as nota[ ] is only valid for [ ]

ti !

o

n

n n

n

n

z

y n x n x nx n z ny n x n z z z z z

x

x n

zn x n z

DD

D D

1

n for theunit-delay operator sym

bol .

[ ] [ ] [ 1]y n z x n x nD

z-Transform 15

Unit Delay

x[n] x[n-1] z -1

X(z) z-1X(z)

×

z-1

× × ×

+ + +

0b 1b 2b3b

[ ]x n [ 1]x n [ 2]x n [ 3]x n

[ ]y n

z-1 z-1

z-Transform 16

Convolution and z-Transform

Domain Domain[ ] [ ] [ ] ( ) ( ) ( )

n zy n h n x n Y z H z X z

Cascading of systems

LTI 1 LTI 2 1[ ] [ ]x n h n[ ]x n 1 2( [ ] [ ]) [ ]x n h n h n

1[ ]h n 2[ ]h n

1 2 1 2[ ] [ ] [ ] ( ) ( ) ( )h n h n h n H z H z H z

z-Transform 17

1 2 3

12 1

1 2 31 2

2 11

1

1

1 2

Example:( ) 1 2 2

One root of this polynomoial is at , and we can write( ) ( )(1 ) ( ) ( )

( ) 1 2 2 1( ) 1

( ) 1 1

11

H z z z z

H z H z z H zH z z z zH z zH z z

H z

z

z z

z

z

11( ) (1 )H z z 1 2

2 ( ) 1H z z z

Partitioning

z-Transform 18

ˆ

ˆ

0 0

ˆ( )

ˆ

(

..

)

.

M Mk j k

k kk k

jz

H z b

e

z H b e

ω

ωω− −

= =

=

= =

⇔

⇔∑ ∑

z - dom

Unit

ain

circle in

the

co

ω-

mplex

domain

plane

z-Transform 19

Poles and Zeros

3 3

3 3

1 2 3

1 1 1

3 2

3 3

e.g.: ( ) 1 2 2

(1 )(1 )(1 )

2 2 1 ( 1)( )( )

j j

j j

H z z z z

z e z e z

z z z z z e z ez z

ο

ο

ο

x3

Triple pole at zero

z-Transform 20

What does a zero mean?

3

1 2 33 3 3 3

23

3

[ ] [ ] 2 [ 1] 2 [ 2] [ 3]

[ ]

[ ] 2 2

(1 2 2 )

1 1 3 1 3 1 0

n

n n n n

n

n

j

j j j j

j j j j

j

y n x n x n x n x n

x n e

y n e e e e

e e e e

e j j

3 31 2 3

The complex input signals

( ) 1, (z ) , (z )get supressed.

n nj jn n nz e e

z-Transform 21

Nulling Filter

0 0

0ˆ ˆ1 1

0 2 2

Zeors in the plan "remove" signalsof the form [ ] .

ˆTo remove a cosine signal cos( )both complex signals must be removed.The zeros are conjugate complex.

n

j n j n

zx n z

n e e

z-Transform 22

PN-Video

z-Transform 23

21 11

10

10 109

1 90

( ) (1 )

10-point running average 101 1( )1 ( 1)

j kL Lk L

kk

k

k

H z z e z

Lz zH z zz z z

z-Transform 24

3 21 2 3

3

2 2 1( ) 1 2 2 z z zH z z z zz

− − − − + −= − + − =

-1 -0.5 0 0.5 1

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

3

Real Part

Imag

inar

y Pa

rt

3 2 /3 /3

/30

2 2 1 ( 1)( ),

01

( )

j j

j

z z z z z e z ez e

π π

π

−− + −

= ±

= − − − =

B = [1 –2 2 -1] A = [1] zplane(B,A)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0

1

2

3

4

5

6

Normalized Frequency ( ×π rad/sample)

Mag

nitu

de

Magnitude Response

z-Transform 25

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0

1

2

3

4

5

6


Mag

nitu

de

Magnitude Response

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

-20

-15

-10

-5

0

5

10

15


Mag

nitu

de (

dB)

Magnitude Response (dB)

lin log

z-Transform 26

Inverse Filtering (1)

Channel

Output signal known

Channel known

Input signal wanted

[ ] [ ] [ ][ ] ?

y n h n x nx n

= ∗= Inverse convolution

Deconvolution

z-Transform 27

Inverse Filtering (2) Solution in the time domain difficult or impossible

Solution in the z-domain

2Correction

1 2 2

1Ch

1

annel

( )( ) ( ) ( ) ( )

1( ) ( ) 1 )

)

(( )

(Y z X z H z X z

H z H z H

z

zH z

zH H= =

= ⇒ =

z-Transform 28

Inverse Filtering (3) 1 2

1

2 1 21

( ) (1 0.5 )1 1H (z) = =

H ( ) (1 0.5 )

H z z z

z z z

− −

− −

= − +

− +

Zeros in denominator (outside of zero) Polstellen

-1 -0.5 0 0.5 1

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

2

Real Part

Imag

inar

y Pa

rt

z-Transform 29

-1

-0.5

0

0.5

1-1

-0.50

0.51

0

2

4

6

8

10

12

14

16

ImRe

H(z

)

11

1 0.9( )

zH z −−

=

-1 -0.5 0 0.5 1

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

Real Part

Imag

inar

y Pa

rtPole/Zero Plot

0 10 20 30 40 50 60 70 80 900

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Samples

Impulse Response

Am

plitu

de

PN-diagram

Impulse response

z-Transform 30

11

1 0.9( )

zH z −−

=

-1 -0.5 0 0.5 1

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

Real Part

Imag

inar

y Pa

rt

0

π

DSP-8-IIR 1

Infinite Impulse Response-Filter

DSP-8-IIR 2

First order IIR-Filter 1 1[ ] [ 1] [ ] [ 1]oy n a y n b x n b x n

FIR block

z-1

x

x + +

x

z-1

b0

b1 a1

v[n] x[n] y[n]

y[n-1] Feed-forward block (FIR-Filter)

Feed-back block

DSP-8-IIR 3

Example 1[ ] 0.8 [ 1] 5 [ ] ( 0)

[ ] 2 [ ] 3 [ 1] 2 [ 3]y n y n x n bx n n n n

Input assumed to be zero prior to starting time n0, x[n] = 0 for n < n0 Output assumed to be zero before starting time of the signal y[n] = 0 for n < n0 System initially at rest

DSP-8-IIR 4

0.8 [ 1] 5(2) 0.8 (0) 5 ( ) 0.8 [0] 5 [1] 0.8 (10) 5 ( ) 0.8 [1] 5 [2] 0.8 ( 7) 5 ( ) 0.8 [2] 5

[0] 10[1] 7[2

[3] ] 5.6

[3]

23

0

yyy

yy x

xxy

yy

5.52[4 4.4416

0.8 ( 5.6) 5 ( ) ] 0.8 [3] 5 [4] 0.8 (5.52) 5 ( )

0.8 [4] 5 [5] 0.8 (4.4416) 5 ( ) 0.8 [5] 5 [6] 0.8

[5] 3.5328

2

(3

0

50

[6] .

yy

y xy xyy x

328) 5 ( ) 2. 6 20 82

0 5 10 15-8

-6

-4

-2

0

2

4

6

8

10

Eingang0 5 10 15

-8

-6

-4

-2

0

2

4

6

8

10

Ausgang

3

3 (Input zero)[ ] 0.8[ ]

[ 1][3](0.8)ny n y

ny n y n

B = [5]; A = [1 -0.8] filter(B,A,x)

[ ] 0.8 [ 1] 5 [ ]y n y n x n

DSP-8-IIR 5

01

[[ ] [ ]]N M

klkl

b xa y n ly kn n

FIR: Output is f(Input)

Feed-forward

Output is f(Output)

Feed-back

M is order of filter for FIR filters, N is order of filter for IIR filters.

DSP-8-IIR 6

Linearity, Time-invariance

IIR-Filter 01

[[ ] [ ]]N M

klkl

b xa y n ly kn n

are linear und time-invariant

DSP-8-IIR 7

Impulse response of a 1st order system

The response to a unit pulse characterizes a LTI system completely.

Input signal represented as superposition of weighted and delayed impulses, output signal constructed from weighted and delayed impulse responses:

[ ] [ ] [ ]k

y n x k h n k

DSP-8-IIR 8

1 0 1 0

1

1

1

1

0

Proof by evaluating:[0] [ 1] [0]

[ ] [ 1] [ ]of impulse response

The

( )(0) 1for 1[

[ ] [ 1] [ ]

( ) für 0[ ]

0 für 0

soluti

] [

on is:o

no

o

h n a h n b n

b

y n a y n

a nh

b

nn

x n

h a h b a b bn

h n a h n

Difference equation

0

11 1 1 0 1

1] [ ]

( )n n no o

b n

b a a b a b a

0

1.

2.

DSP-8-IIR 9

Notation using step response

1 1

1

1 for 0[ ]

0 for 0hen[ ] ( ) [ ]

cen

oh a

nn

n

n b n

DSP-8-IIR 10

Step response 1[ ] [ 1] [ ]oy n a y n b x n

Iterate difference equation to calculate output sequence one sample at a time:

0

0 0 12

0 0 1 0 1

2 30 1 1 1

[ ] [ ] 1 1 ( )

01

2

3

1 ( ) ( )

1 1

x n y nbb b ab b a b a

b a

n

a a

DSP-8-IIR 11

21 1 1 0

1

11

0 1

k

1

10

L

=0

[ ] (1 ... )

We recall

1 11

1 1

for 0, if 1[ ]1

1

no

nk

k

n

L

ky n b a a a b

r rr

L

ay

r

nb

a

ana

r

DSP-8-IIR 12

We identify three cases:

11

11

0

1

1

1If 1 , than decays

If 1 , than dominate

to zero as

stab

s and [ ] getslarger without bound

l

1 [ ] ( 1) output [

e system lim [

==> unstable syst

]1

emn

n

n

o

a a y n

ay n n

a a nby n

a

b y

] grows as 1

[ ] if even [ ] 0 if odd

o

n na

y n b ny n n

11

01

1[ ]1

nay n ba

DSP-8-IIR 13

0 5 100

0.5

1

1.5

2

a = 0.5 ... stabil

Am

plitu

deStep Response

0 10 200

50

100

a = 1.1 ... instabil

Am

plitu

de

Step Response

0 5 100

0.5

1

a = -1 ... Grenzfall

Am

plitu

de

Step Response

0 5 100

5

10

a = 1 ... instabil

Am

plitu

de

Step Response

DSP-8-IIR 14

System function of an IIR filter

Domain Domain[ ] [ ] [ ] ( ) ( ) ( )

n zy n h n x n Y z H z X z

The system function of an FIR system is always a polynomial in z

-1. If the difference equation has feedback as in IIR systems, the system function is the ratio of two polynomials (rational function).

DSP-8-IIR 15

1 11 1

1

11

1

1

11

11 1

[ ] [ 1] [ ] [ 1]

( ) ( ) ( ) ( )

( ) ( ) ( ) ( )

( ) ( )( )( 1) ( )

o

o

o

o

y n a y n b x n b x n

Y z a z Y z b X z b z X z

Y z a z Y z b X z b z X z

Y z B zH zX

b b zz z Aa z

Numerator polynomial: Feed-forward coefficients Denominator polynomial: 1 + negative feed-back coefficients

DSP-8-IIR 16

Block diagram 1st Direct form

1 11 1

111

11( ) ( )1

1 1 ( )o

oa z a zb b zH z B z

Ab b z

z

DSP-8-IIR 17

Block diagram 2nd Direct form

1 1( ) ( )( ) ( )

B zz

BA A

zz

DSP-8-IIR 18

Delay elements combined

DSP-8-IIR 19

Poles and Zeros 1

1 11

1 1

1

1

( )1

Zero

Pole

o o

o

b b z b z bH za z z a

bzb

z a

DSP-8-IIR 20

Poles and Stabilität

1

0

0

11 1 1

11

11 1

11 1 11 1

0 00

1

The system function

( )( )

1leads to the impulse response

[ ] ( ) ( )

[ ] [ 1]

no

bo bo

n no

nb n

b b a a n

b zb b zH za z z a

h n b a bn na

DSP-8-IIR 21

1

1

1

The impulse response is proportional to for 1.

The response if and 1 .

The impulsThe location of pole(s) indicat

e responsees decaying

or growing impulse respo

if 1

se.S

.

n

na n

n a

a

decayes

grows

ystems with decaying impuls response are stable systems.

The location of poles of stable system functions is strictly inside the unit circle of the z-plane.

DSP-8-IIR 22

Frequency response of an IIR-filter

11( )

1 0.8H z

z

-1

-0.5

0

0.5

1-1

-0.50

0.51

0

2

4

6

8

10

12

14

16

ImRe

H(z

)

ˆ

ˆ

ˆ

ˆ[ ] ( )ˆ( ) ( ) ( ) j

j n

jz e

y n eH e H z

H

H

If sinus sequence meets pole at the unit circle: bounded input unbounded output (resonance)

DSP-8-IIR 23

z-domain

Frequency response (lin)

Time-domain

DSP-8-IIR 24

ˆ[ ] ( ) ( )jh n H z H e

Poles, Zeros H(z)

Frequency domain

Time-domain Input, Output

h(n) ˆ( )jH e

, k ka b

?

DSP-8-IIR 25

1 2 0 1 2

1 20 1 2

1 21 2

ˆ ˆ2ˆ 0 1 2

ˆ ˆ1 2

[ ] [ 1] [ 2] [ ] [ 1] [ 2]

( )1

( )1

j jj

j j

y n a y n a y n b x n b x n b x n

b b z b zH za z a z

b b e b eH ea e a e

Solution difference equation?

DSP-8-IIR 26

Inverse z-Transform

10 1

11

We consider a first order system

( ) ( ) ( ) ( )1

1. Determine transform X( ) of input [ ]2. Muliply ( ) ( ) to get Y(z)3. Determine inverse transform of ( )

b b zH z Y z H z X za zz z x n

H z X zz Y z

DSP-8-IIR 27

1

1

0

1

0

1

1

1

Step response of a first order system[ ] [ ]

( )

this sum is finite for 1

1( ) ... für 1

U

1[ ] 1

nit step for 1

n

nn n

n n

na

h n a n

H z a z az

az

H z a zaz

a

naz

0

1 für 11

k

k

x xx

DSP-8-IIR 28

1

11 1

11

1 21 1

11

11

1

1

1

1

1

( ) ( ) ( )1 (1 )

(determine and by comparing coefficients)1

or faster using residue meth

1

1

(1(1 )(1

1

ode

)

1oo

o

b b zY z H z X za z a

b b za zz

A B A Bz

b b

a z

z BY z Az

a z

z

11 1

1

1

11

1 1

1 11 1 1

11

11 1

1

1 11

1)

(1 )( )1 1

(1 )

1 1

oz a

z a z a

o o

z a

zb b z BY z A A

z z

a

b b z b b aA

aa

a

z

zz

z

Partial fraction expansion

DSP-8-IIR 29

1 0 11

1

11 1 0 1

1 1 111 1

( ) (1 )1

[ ] [ ] [ ]1 1

z

no

b bB Y z za

b b a b by n a n na a

11 [ ] 1

nA Aa naz

DSP-8-IIR 30

Inverse transform (M<N)

1

11

1

11

1. Factoring numerator polynomial of ( ) (1 ) für 1,2,...,2. Partial fraction expansion

( )1

( )(1 )

3. Inverse transform

[ ] [ ]

k

k

Nk

k k

k k z p

Nn

k kk

H zp z k N

AH zp z

A H z p z

h n A p n

DSP-8-IIR 31

Important transform pairs

0

1 2 1 2

0

0

[ ] [ ] ( ) ( )

[ ] z ( )[ ] [ ] [ ] ( ) ( ) ( )[ ] 1[ - ]

n

ax n bx n aX z bX z

x n n X zy n x n h n Y z X z H z

nn n

0

1 1

1[ ]

1

n

n

z

a naz

DSP-8-IIR 32

1 1

1 2 1 1

1 1

1 2.1 1 2.1( )1 0.3 0.4 (1 0.5 )(1 0.8 )

( )1 0.5 1 0.8

z zX zz z z z

A BX zz z

1 1 11

1 1 1

0.50.5

1 2.1 1 0.5 (1 0.5 )( )(1 0.

0

5 ) 21 0.8 1 0.5 1

10.8

zz

z z B zX z z Az z z

11

10.8

1 2.1( )(1 0.5 ) 11 0.8 z

zB X z zz

1 1[ ] 2( 0.5) [ ] (0.8) [ ]n nx n n n

DSP-8-IIR 33

A = [1 -0.3 -0.4]; B = [1 -2.1]; [R,P,K] = residuez(B,A) R' = -1 2 P' = 0.8 -0.5 K = [] impz(B,A)

0 5 10 15 20 25 30 35 40-2

-1.5

-1

-0.5

0

0.5

1

n (samples)

Am

plitu

de

Impulse Response

DSP-8-IIR 34

ˆ

ˆ ˆ

0 0ˆ

ˆ

0

Frequency response of FIR-filters

[ ] [ ]

ˆ = ( ) Frequency response with z-transform

ˆ( ) ( ) ( )

e.g.: ( )1

j

M Mj k j n

k kk k

j n

jz e

y n b x n k b e Ae

Ae

H e H z

bH z

H

H

0 0 0

0

ˆ 0ˆ1

1 1

ˆ ˆ ˆ0ˆ

1

( )1

[ ] ( ) 1

jj

j j n j nj

bH ea z a e

by n H e e e na e

Frequency response of IIR-Filter

DSP-8-IIR 35

0

0

0

0 0

ˆ1 ˆ 1

0ˆ1 1

1

0 1 0ˆ ˆ

1 1ˆ1

1

Suddenly applied complex exponential input sequence:1[ ] [ ] ( )

11( ) ( ) ( )

1 1After partial fraction expansion:

1( )1 1

j nj

j

j j

j

x n e n X ze z

bY z H z X za z e z

b a ba e a eY z

a z e

0 1z

Steady state - transient response

DSP-8-IIR 36

0 0

0 0

0

0

1 1

1 1

decays if stable steady st

0 1 0ˆ ˆ

1 1ˆ

1

ˆ1

ate

0 1 0ˆ ˆ

1 1

( )1 1

[ ] [ ] [ ]

1[ ]

1

1-

1j

nn

n

j

j

j jj

Y zz z

y n n n

u nz

a e

b a ba e a e

b a ba

a

aa

e ee

a

1

DSP-8-IIR 37

( )

( )

1 10.2

0.3

0.

3131

0.2 0.29712 1

1

1

0.22

1

0

0.9[ ] [ ]0.9

0.4980 [ ]

1[ ] [ ] 0.5533 [ ]1 0.9

0.5533 ( 0.2971) [ ]0.55

0.

cos 0.2sin 033 ( 0.2

1;

9

0.

71)

9

.2

9

0.9

j

j

j n j j

n

jn

n

y n ne

e n

y n e n e ne

nj

b

e

nn

a

π

π ππ

δ

δ

δ

ππ

δ

δ

−

−−

− −−

−

− = = − −

= = = +

−

−

= ++ +

= = −

1

1 2

[ ][ ] [ ] [ ]

ny n y n y n

δ−

= +

DSP-8-IIR 38

0 5 10 15 20-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

b1 = 1 a1 = -0.9

Am

plitu

de

Impulse Response

DSP-8-IIR 39

1Instable system: Pole at 1.1a =

DSP-8-IIR 40

IIR-system second order

1 2 1 2

1 2 1 21 2 1 2

1 21 2

1 21 2

[ ] [ 1] [ 2] [ ] [ 1] [ 2]

( ) ( ) ( ) ( ) ( ) ( )

( )1

o

o

o

y n a y n a y n b x n b x n b x n

Y z a z Y z a z Y z b X z b z X z b z X z

b b z b zH za z a z

DSP-8-IIR 41

Poles and Zeros 1 2 2

0 1 2 0 1 21 2 2

1 2 1 2

( )( )( ) 1

Y z b b z b z b z b z bH zX z a z a z z a z a

A polynomial of degree N has N roots. For real coefficients of the polynomial the roots are either real or complex conjugate

DSP-8-IIR 42

Impulse response 1 2

0 1 21 2

1 2

2 1 21 1

2 1 2

20 1 1 1 2 2 1

2

( )( )( ) 1

( )1 1

[ ] [ ] [ ] [ ]n n

b b z b zY zH zX z a z a z

b A AH za p z p z

bh n n A p n A p na

DSP-8-IIR 43

Real poles

1 2 1 15 1 1 16 6 2 3

1 11 12 3

1 11 12 3

1 1( )1 (1 )(1 )

3 2( )(1 ) (1 )

[ ] 3 [ ] 2 [ ]nn

H zz z z z

H zz z

h n n n

1 2n

k

p pp

For real und the impulse responseconsist of two functions in the form

DSP-8-IIR 44

5 16 6B =[1]; A = 1,- ,

impz(B,A)

0 2 4 6 8 10 120

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

n (samples)

Am

plitu

deImpulse Response

DSP-8-IIR 45

Complex conjugate poles

DSP-8-IIR 46

Complex poles at the unit circle 1 1

/ 4 1 / 4 1 2 2

1 1( )(1 )(1 ) 1 1.4142j j

z zH ze z e z z zπ π

− −

− − − − −

+ += =

− − − +

1.1781 1.1781

/ 4 1 / 4 1

1.3066 1.3066( )(1 ) (1 )

j j

j j

e eH ze z e zπ π

−

− − −= +− −

( ) ( )/ 4 / 41.1781 1.17811 1[ ] 1.3066 [ ] 1.3066 [ ]j n j nj jh n e e n e e nπ πδ δ−−− −= +

1[ ] 2 1.3066cos 1.1781 [ ]4

h n n nπ δ− = ∗ −

DSP-8-IIR 47

Complex conjugate poles at unit circle Second-order oscillator

DSP-8-IIR 48

Complex poles inside the unit circle 1 1

2 / 4 1 / 4 11 12 2

1.2490 1.2490

/ 4 1 / 4 11 12 2

1 1( )1 1 0.5 (1 )(1 )

1.5811 1.5811(1 ) (1 )

j j

j j

j j

z zH zz z e z e z

e ee z e z

π π

π π

− −

− − − −

−

− − −

+ += =

− + − −

= +− −

-1 -0.5 0 0.5 1

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

Real Part

Imag

inar

y Pa

rt B = [1 1]; A = [1 -1 0.5] zplane(B,A)

DSP-8-IIR 49

1[ ] 2 1.5811 cos 1.24942

n

h n nπ = × −

0 5 10 15 20 25-0.5

0

0.5

1

1.5

2

n (samples)

Am

plitu

deImpulse Response

1

2

1( )1 1 0.5

zH zz z

−

−

+=

− +

DSP-8-IIR 50

10

0 1 1 20

10

0 1 1 20

10

0 1 1 2 20

10

0 1

(sin )sin [ ] 1 2(cos )

1 (cos )cos [ ] 1 2(cos )

1 ( cos )cos [ ] 1 (2 cos )

( sin )sin [ ] 1 (2

n

n

zn nz z

zn nz z

r zr n nr z r z

r zr n n

0 0

1 2 20

0 1 1 1

10 1

0 1 1 2 2

cos )

0.5 0.5cos( ) [ ] 1 1

cos( ) [ ] 1 2

j jn

j j

n

r z r zKe KeK r k n

re z re zb b zK r k naz r z

2 2 21 0 0 1 0 1

02 2 2 2

2 arccos arctanb r b ab b ab baKr a r r a

DSP-8-IIR 51

PZ-Video

Running average FIR filter

SigProc-7-IIR 52 Direct form I Direct form II

1 1( ) ( )

( ) ( )B z B z

A z A z

SigProc-7-IIR 53 Transposed direct form II

Attention!

SigProc-7-IIR 54

IIR filters are often represented in the form

In such cases the coefficients a in the block diagrams are negative!

1

01

0

( )( )( ) 1

Lk

kk

Mm

mm

b zY zH zX z a z

DSP-8-IIR 55

Filterdesign

SPtool

DSP-8-IIR 56

• Sine wave 2 kHz + noise (randn) (average = 0 , variance = 0.01)

• Band pass ripple 1.5 dB, stop band 35 dB

SigProc-6-FIR 57

SigProc-6-FIR 58

SigProc-6-FIR 59

No evidence of noise in filtered signal spectrum

SigProc-6-FIR 60

Spectrum filter input

Spectrum filter output

Sampling Theorem

DSP_9-Sampling Theorem 2

Sampling

ˆ[ ] ( ) cos( ) cos( )

ˆ

s s

s

x n x nT A nT A n

T

Normalized angular frequency


Normalized angular frequency

ˆˆ

ω has the dimension rad / sec, ω = ωT has the dimension rad,i.e. !The time axis is lost after sampling,the discrete signal is merely a sequence of numbers anddoesn't c

ω is dimen

ontain inf

sionl

orm

ess

ation about the sampling period.To reconstruct the original signal the sampling frequency must be known.I An infinite number of continuous sine signals is transformed into the n other

same diswords :

crete sine representation.


Number of samples

0 1 2 3 4 5

-2

0

2

Zeit t

Am

plitu

de

0 1 2 3 4 5

-2

0

2

Index

Am

plitu

de

0 1 2 3 4 5

-2

0

2

Index

Am

plitu

de


Signal f(t) multiplied by a sequence of impulses δT(t) with time-step of T seconds (sampling time).

( ) ( ) ( ) ( ) ( )S Tk

f t f t t f kT t kT

The impulse sequence is a periodic function and can be represented as Fourier series.

1( ) [1 2 cos 2 cos2 2 cos 3 ...

2 cos ] ( )

T s s s

s

t t t tT

k t k

Representation of the sampling process


/ 2

/ 2

11 ( )

1 2( )

ss

s

s

s

T jk tk TT

s

jk tT s

ks s

D t e dtT

t eT T

Tω

ω

δ

πδ ω

−

−

∞

=−∞

= =

= =

∫

∑

1( ) [1 2cos 2cos2 2cos3 ... 2 cos ]

( )

T s s s ss

t t t t k tT

k

kδ ω ω ω ω= +

∞

+ + + +

→ (real representation)

-5 0 5 100

0.5

1

x 1/Ts

complex representation


1

( ) ( ) ( ) [ ( ) 2 ( )cos 2 ( )cos2 ...]S T s sf t f t t f t f t t f t tT

Sampled cosine (ω0) signal

0 0 0

1cos( )cos( ) cos( ) cos( )

2s s st t t t

Sampling generates for spectral component

0( )s

0


1( ) ( ) ( ) [ ( )cos 0 2 ( )cos 2 ( )cos2 ...]S T s sf t f t t f t t f t t f t t

T

e.g. :

0 0

0

0

( )

1cos( ) cos( ) cos(

co

)2

s

cos( )s s s

f t

t

t

t t t

ω

ω


Spectrum of sampled signal

( ) ( )1( ) cos [ ( ) ( )]2

( ); ( 2 ), ( 2 ); ( 3 ), ( 3 );...1( ) ( )

s s s

s s s s

s sns

f t F

f t t F F

F F F F F

F F nT

ω

ω ω ω ω ω

ω ω ω ω ω ω ω ω ω

ω ω ω∞

=−∞

⇔

⇔ − + +

− + − +

= −∑

The spectrum of the sampled signal is periodic with period ωs.


The spectrum of the original signal f(t) is contained in the spectrum of the sampled signal fs(t) and can be recovered from Fs(ω) by „cutting out" using an (ideal) low pass filter.


Distance of spectra depends on sampling frequency


Shannon Sampling Theorem

A continuous-time signal with frequencies not higher than fmax can be reconstructed exactly from its samples, x[n] = x(nTs), if the samples are taken at a rate fs = 1/Ts that is greater than 2fmax.


Examples

• Scale 14 kHz [44,1 kHz] • Scale 14 kHz [11 kHz]

• Phasor • Helicopter


Spectrum of sampled sine function Sampling frequency = 90 Hz

100 Hz -100 Hz

f = 10 Hz ! Aliasing


Sampling frequency 110 Hz 100 Hz

-100 Hz

f = – 10 Hz ! Folding


Sampling frequency = 250 Hz 100 Hz

f = 100 Hz

-100 Hz

Selectivity of filter


Sampling frequency = 350 Hz

Filter can be less selective

100 Hz

-100 Hz


Spectrum of a continuous-time signal

Spectrum of the sampled continuous-time signal (continued periodically)

Aliasing (Folding), if sampling theorem violated.


Reconstruction/Interpolation

Conversion discret => continous(Interpolation with pulses)

( ) [ ] ( )

( ) characteristic pulse shape of the converter

sn

y t y n p t nT

p t


Interpolation in the time domain Zero-Order Hold Linear Interpolation


1 x

4 x


Interpolation in the frequency domain


Ideal Filter

• To reconstruct the original signal from the (periodic) spectrum, a filter with rectangular frequency response would be required (ideal filter)

• The input signal of the filter is the sequence of impulses of the sampled signal

• The output signal is the superposition of weighted and time-delayed pulse responses


Impulse response of the ideal low pass

( )2

2

1 | | 2( )

0 | | 21( ) 2 sinc 2

2B j t

B

BH

B

h t e d BT Btπ ω

π

ω πω

ω π

ω ππ −

<= ≥

= =∫

00

H(ω

)

ω

0

0

1

h(t)

Zeit t

-2 πB

1/2B

2πB

T

4/2B3/2B2/2B-4/2B -3/2B -2/2B -1/2B


The output signal is the superposition of time-delayed, weighted impulse responses

( )(( ) ( )) ( )sinc 2k k

kTf kT f kTf t h t Bt kππ= − = −∑ ∑

For clarity only every second sinc-function drawn.


The interpolation with sinc-pulses reconstructs the continuous-time signal exactly.

But the impulse response of the ideal low pass is non-causal, i.e. the filter responds before the pulse is applied.

It is not possible to realize non-causal filters!


Interpolation filters (1) Ideal (analog) low pass filters are non-causal and are therefore not realizable.

A practical solution is to use sampling frequencies higher than the Nyquist-Frequency. In this case the gaps in the periodic spectrum are wider and the selectivity of the interpolation-low-pass-filter can be less than ideal. (The continuous-time signal is not reconstructed exactly.)


Interpolation filters (2) It is possible to design sharp (near ideal) analog filters with high attenuation in the stop band. But it is not possible to design filters which fully repress signals in the whole stop band.

Every practical signal is of finite duration. We know from Fourier-Transform, that a signal with finite duration has spectral components in the frequency range from −∞ to ∞.


No signal can be time-limited and frequency-limited at the same time:

Time-limited (finite duration) means that the spectrum reaches from – ∞ bis ∞ (is not band-limited).

Band-limited (finite bandwidth) signals have infinite signal duration (not time-limited).


Overlapping Spectra


Antialiasing-Filter To avoid overlapping of spectra Antialiasing-Filter limit the bandwidth of signals.


Digitization

The resolution of the Analog-/Digital-Converter depends on quality requirements.

For voice signals a resolution of 8 bit is sufficient, music signals on Audio–CD use 16 bit resolution.

The lower the resolution of the A/D converter the more the digital signal deviates from the analog signal.


Quantization error

0 20 40 60-1

0

1

Zeit t

f(t)

0 20 40 60-1

0

1

Zeit t

f(t)

0 20 40 60-1

0

1

Zeit t

F(t)

0 20 40 60-1

0

1

Zeit t

F(t)

0 20 40 60-0.5

0

0.5Fehler: ∆ = f (t) - F(t)

Zeit t

∆

0 20 40 60-5

0

5x 10

-10 Fehler: ∆ = f (t) - F(t)

Zeit t

∆ Quantization noise


Quantization noise (1) Quantization error

For equal probability of all signal amplitudes we get the average power of the error signal (noise)

/ 22 2

/ 2

112

~ 0.2912 12

eff

eff

e s ds

LSBe LSB

∆

−∆

∆= =

∆∆

= =

∫

8 bit: 1/900 12 bit: 1/14.000 16 bit: 1/226.000

of value range


Quantization noise (2)

[y,fs,nbits]=wavread('file.wav')

sound(y,fs,6) % Play back with 6 bits

wavwrite(a,fs,nbits,'wavefile.wav')

nbits = 8, 16, 32 or 64


Sampling in the frequency domain Sampling in the time domain: band-limited signals Sampling in the frequency domain: time-limited signals


Discrete Fourier-Transform (DFT)

We calculate the DFT of a time-limited signal f(t) with duration (a). The spectrum F(ω) is not band-limited (b).

From the continuous time-limited signal f(t) we get the

discrete signal fs(t) by sampling with the interval T=1/Fs (c). Due to sampling the spectrum becomes periodic with period Fs=1/T and we get spectrum Fs(ω) (d).

Sampling the spectrum with F0=1/T0 (f) leads to the periodic

time-signal width period T0 (e).


00

00

1

001

0

1[ ] [ ]

[ ] [ ]

Njk n

k

Njk n

n

f n F k eN

F k f n e

−Ω

=

−− Ω

=

=

=

∑

∑

DFT

IDFT

Diskrete Fourier-Transformation


0 0

0

0

1 1ˆ

0 0

1( ) [ ] [ ] ( )

1( ) ( ) ( ) ( ) 2

1[ ] [ ] [ ] [ ]

Tj j

j j

N Nj

k t t

t t

n

n

k

k

k

k j

f F e F f e dtT

f F e d F f e dt

t

f F e F

t

t t

n nk e

k

k fN

kω ω

ω

ω ωωπ

ω ω

+∞−

=−∞

∞ ∞−

−∞ −∞

− −−

= =

= =

= =

= =

∑ ∫

∫ ∫

∑ ∑

⇔

⇔

⇔ 0ˆ

0

2ˆ

k n

TN

ω

πω ω= =

Time domain Frequency domain

FR:

DFT:

FT:


Spectral lines of DFT • The spectral lines of the DFT lie on a raster

• There are no spectral lines outside the raster, e.g. 15 Hz. The DFT algorithm must compose a 15 Hz signal from other frequency components!

0

0

1000; e.g.: 10 Samples: 0, 10, 20, ..., 90 Hz100

s

s

T fNT f

= = =


( )

( ) ( )

( )

( )

( )( ) ( )

0

0

0

0

0

ˆ

N-1ˆ 2 /

n=0

ˆ2 /

ˆ2 /

Phase onl

ˆ 02

y

2 / 1 /

[ ] 0,1,2,..., 1 [ ]

X[k]= ...

1 1

ˆsin 2 / / 2

j n

j n j N kn

j k N Nj

j k N

j k N Nj

x n e n Nx n

e e

eee

k N Ne e

ω ϕ

ω ϕ π

π ωϕ

π ω

π ωϕ π ω

+

+ −

− −

− −

− − −

= = −

= =

−= =

−− =

∑

The N- point DFT of is

( )0

Diriclet Function

ˆsin 2 / / 2k Nπ ω−

X[k] contains samples of the Dirichlet Function.

2 2 2

2 sin

1j j jj

j

e e e e

0 1 2 3 4 5 6-0.4

-0.2

0

0.2

0.4

0.6

0.8

1


0 0.5 1 1.5 2-1

-0.5

0

0.5

1Time-Domain (f = 3, T = 1, fs = 8, N = 8)

t

x(t)

0 2 4 6 8 10 12 14 160

1

2

3

4Frequency-Domain

f

|X(f)

|

Sampling of f = 3 delivers a non-zero value only at this frequency. (Cosine function continued periodically.)


The frequency f = 2.2 doesn‘t lie on the raster and must be composed of other frequencies! The Dirichlet Function delivers non-zero values at the raster frequencies. Note the periodic continuation of the cosine function!

0 0.5 1 1.5 2-1

-0.5

0

0.5

1Time-Domain (f = 2.2, T = 1, fs = 8, N = 8)

t

x(t)

0 2 4 6 8 10 12 14 160

1

2

3

4Frequency-Domain

f

|X(f)

|


In the previous representation we calculated the spectrum of ( )0ˆ[ ] j nx n e ω ϕ+=

For the real-valued signal x[n] we must compute the spectrum of

( ) ( )( )0 0ˆ ˆ[ ] j n j nx n e eω ϕ ω ϕ+ − += +

(by adding X [k] for .) 0ω−


0 0.5 1 1.5 2-1

-0.5

0

0.5

1Time-Domain (f = 2.2, T = 1, fs = 8, N = 8)

t

x(t)

0 2 4 6 8 10 12 14 160

1

2

3

4

Frequency-Domain

f

|X(f)

|


We always compute the spectrum of the periodic continuation! The periodic continuation delivers as sine-(cosine-)function, and therefore single spectral line, only if the period fits exactly in the raster, i.e. only for spectral lines on the raster.

Mathematische Methoden des Visual Computing …...Phase (1) [ ] [ ] Imaginary Part Real Part 1 1 1 1...

Documents

Transcript of Mathematische Methoden des Visual Computing …...Phase (1) [ ] [ ] Imaginary Part Real Part 1 1 1 1...