SMpA 2001

Analisi dei Segnali

Spettri di ampiezza;

Densità spettrale di potenza;

Stima della risposta armonica

Funzione di coerenza

Relazioni con le funzioni di correlazione

Cepstrum di ampiezza e di potenza; Coefficienti cepstrali

SMpA 2001 n. 2

Fourier Series

( ) ( ) kj tk

ky t Y e ωω

∞

=−∞= ∑

( ) ( )y t y t nT= +

/ 2

/ 2( ) ( ) k

T j tk T

Y y t e dtωω −−

= ∫

Given a periodic signal:

This is its Fourier series:

The coefficients are computed as:

1

1

2

2 2

h h

Tπ

ω

π π

ω ω ω+

= =

= =− Δ

1 1cos(2 ) cos(6 )3 2

t tπ π− +

1/2 1/2

-1/6

f

13

SMpA 2001 n. 3

Sines and Rotating Vectors

A real sine signal is the sum of two complex rotating vectors

-ωtωt

R

I A/2

( )cos( ) ;2

j jAA e e tθ θθ θ ω ϕ−= + = +

That’s why:

In signal analysis wehave to considerpositive and negative frequencies

Periodic signals havemultiple harmonics(so harmonics’vectors make 2, 3,... rotations as the fundamental rotatesonce)

SMpA 2001 n. 4

What does it Means ? ( ) j ty t e ω−

[ ]1 3( ) cos( ) cos(3 )y t Y t Y tω ω= +

[ ]

( ) ( )

( ) ( )

3 31 3

3 3 331

2 4 631

( ) cos( ) cos(3 )

2 2

12 2

j t j t

j t j t j t j t j t

j t j t j t

y t e Y t Y t e

YY e e e e e

YY e e e

ω ω

ω ω ω ω ω

ω ω ω

ω ω− −

− − −

− − −

= + =

⎡ ⎤= + + + =⎢ ⎥⎣ ⎦⎡ ⎤= + + +⎢ ⎥⎣ ⎦

/ 2 3 3/ 2

( ) 0 0 02

T j tT

Yy t e dtω−−

= + + +∫

We have a signal

Make the product

Integrate itZero, because anywaythey rotate an integernumber of times in T

It’s the same thathaving the R-I framerotating at –3ω

SMpA 2001 n. 5

Power

20

1 ( )T

P y t dtT

= ∫Power of a periodic signal

For a sine2 2 2

sin( )2 2 2k k k

k kY Y YY t Pω −⎛ ⎞ ⎛ ⎞= + =⎜ ⎟ ⎜ ⎟

⎝ ⎠ ⎝ ⎠But for d.c.

as the two lines coincide

20P Y=

So we have 2 20

1 ( ) ( ) / 4T

kk

P y t dt YT

ω∞

=−∞= = ∑∫ It holds because

Harmonics are orthogonal over a period.

It is a first form of theParseval’s Theorem1 1cos(2 ) cos(6 )

3 2t tπ π− +For

1 1 1 1 14 4 36 36 4

P ⎛ ⎞ ⎛ ⎞= + + + +⎜ ⎟ ⎜ ⎟⎝ ⎠ ⎝ ⎠

SMpA 2001 n. 6

Fourier Transform 1

Let: 0T ω→ ∞ ⇒ Δ →

2

1( ) ( ) ( ) ( )2

( )

kj t j tk

k

j f t

y t Y e y t Y e d

Y f e df

ω ω

π

ω ω ωπ

∞ ∞

−∞=−∞

∞

−∞

= → = =

=

∑ ∫

∫

/ 2

/ 2( ) ( ) ( ) ( )k

T j t j tk T

Y y t e dt Y y t e dtω ωω ω∞− −

− −∞= → =∫ ∫

It exists iff y(t) has a finite energy Actually here Y(ω) is a density, the “amplitude” from ω to ω + dωis Y(ω)dω(the transform of a sine is a coupleof pulses)

2 ( )y t dt∞

−∞< ∞∫

Are you really interested in phenomenaextending from t= –inf to t=+inf ?

SMpA 2001 n. 7

Fourier Transform 2

2( ) ( ) j f ty t Y f e dfπ∞

−∞= ∫2( ) ( ) j f tY f y t e dtπ∞ −

−∞= ∫

The transform pair:

shows a perfect symmetry, what applies in one direction, applies also in the reverse one

Fourier series can be read as a Discretized Fourier Trasform, and we have:y(t) periodic Y(ω) discrete

Therefore, reversing direction:y(t) discrete Y(ω) periodic

SMpA 2001 n. 8

Sampled Signals

2( ) ( ) nj f tn

nY f y t e π

∞−

=−∞= ∑

/ 2 2/ 2

1( ) ( )s n

s

f j f tn f

sy t Y f e df

fπ

−= ∫

1/fs=Tsy(t)

t

fs is the sampling frequency

0 tn

Ts = tn+1 – tn is the sampling time

Sampling Theorem (Shannon):

fs > 2 fMAX , the maximum frequency in y(t)

Not necessarily periodic

SMpA 2001 n. 9

Working with Computers...

... all is sampled and limited*: both the signal and the “transform”y(t)

y(tk)

T

Ts

Y(ωk)

fs/2 fs

Δf

1/S sT f=

R sT NT=

sf M f= Δ

1R

NTM f

=Δ

* We are free to choosethe record lenght TR.Here is TR=3T

SMpA 2001 n. 10

FFT

T NTN T TR S

s R= = =; Δω π π2 2

Y j y t e dtj t( ) ( )ω ω= −−∞∞z

Y j m T y nt es sj T nm

n N

Ns( ) ( )

/

/Δ Δω ω= −

=−

−∑

2

2 1

When all is sampled and limited, the true Fourier transform cannot be computed,only its discrete approximation as a series(even if the signal is not periodic).

We can only acquire N samples of the signal with a sampling frequency fsand compute M samples of the integral.

FFT is a very fast algorithm to compute it. FFT requires N=M=2q

SMpA 2001 n. 11

Constraints

1RT

f=

ΔGenerally we use FFT,

so N = M

Other constraint:Shannon theorem fs > 2 fMAX

We need to analyze a signal up to a frequency fMAX

with a resolution ΔfTherefore:

All works fine ifthe signal is periodic AND

record lenght is TR = k T, k integer

(never exactly verified)

So we have N = M = TR fs samples both in time AND in frequency(but half of the frequency samples are useless!)

SMpA 2001 n. 12

Experimental Set-up

anti-aliasfilter FFT

sampling

y(t) y(m Δt) Y(n Δω)

Time (sec)

0 20 400

1

2

3

4

0 2 40

0.5

1

Periodicity

fMAX = fs * 15/16 = 9.375

Δf = FMAX / #bins = 9.375 /(16-1) = 10/16 = 0.625

1 6 11 160

1

2

3

4

0 5 10Hz

*

binsfs = 10

0 0.5 1 1.60

0.5

1

G.U.1

Diapositiva 12

G.U.1 La fig di sn è sbagliata. Analogamente all'altra, dopo l'ultimo (16°) campione andrebbe considerato ancorra un DeltaT, quindi Tr=1.6.Così Deltaf = 1/Tr = 0.625Giovanni Ulivi; 28/10/2003

SMpA 2001 n. 13

Satisfied Constraints

0 1 2

-1

0

1

y=sin(2*pi*t)+0.5*sin(4*pi*t)

t = 0: 0.05: 2-0.05Fs = 20Hz, Δf = 0.5Hz

k

1 3

-20

-10

0 10 20 30 40

10

20

0 2 4 6

1 2 Hz

time samples samples

Length of record TR=2 T

SMpA 2001 n. 14

Unsatisfied Constraints

0 1 2

-1

0

1

y=sin(2*pi*t)+0.5*sin(4*pi*t)

t = 0: 0.05: 2.25-0.05Fs = 20Hz, Δf = 0.444Hz

0 5 15 25 35 45

10

20

1 3 5

0.888Hz

2.22HzLength of record TR=2.25 TIndeed we often have TR ≠ KT

SMpA 2001 n. 15

What’s happening?

The lenght of the record is different from a multiple of the periodTR ≠ KT, K integer

No frequency sample (bin) coincide with one of the signal frequenciesConsidering the interpretation of

no signal component makes exactly K turnsThe series is actually associated to a signal with period TR that is:

( ) j ty t e ω−

0 10 20 30 40 50 60 70 80 90-1.5

-1

-0.5

0

0.5

1

1.5

SMpA 2001 n. 16

How to Study the Problem?

Only a time segment is available for y(t) (from –TR / 2 to TR / 2)

Differently restating the problem in continuous time:

in the frequency domain, a convolution

Should the signal spectrum be an impulse,its estimate would be:

Difficult to single out near components andto detect small components

“Rectangular window”

Y j y t e dt w t y t e dtTj t

TR

TR j t( ) ( ) ( ) ( )//ω ω ω= = ⋅−

−−

−∞∞z z2

20 w t

t T

t T TR

R R0

0 2

1 2 2( )

, /

, / /=

>

= −RST L

Y j Y j W jT ( ) ( ) ( )ω ω ω= ⊗ 0worst casesamples’position

SMpA 2001 n. 17

How to solve it?

… we can try smoothingthe extremities of the record,by scaling it by a “time window”

As the problem ariseat segment stitching, …

Record Replica

SMpA 2001 n. 18

Characteristics

1

0.707

Bandwidth f

Max sidelobeDecay rate

Side-lobe(dB)

decay(dB/ dec)

Band-width

Max |•| error (dB)

Rectangular(boxcar)

-13 -20 1.00 3.9

Hanning -32 -60 1.50 1.4

Hamming -43 -20 1.36 1.8

Kaiser-Bessel -69 -20 1.80 1.0

Flat-top -93 0 3.77 <0.01Time window transform

SMpA 2001 n. 19

Time Windows

-100 -50 0 50 1000

5

10

15

boxcar(15/512)

-100 -50 0 50 1000

10

20

30

boxcar(30/512)

-100 -50 0 50 1000

5

10

15

hamming(30/512)

bins

For a boxcar, the first zero is in f = 1/TR (Hz)

Sidelobes decay is related to continuity: if w(t) is C n-1 decay rate is ω-n

Greatest bandwidth ↔ Reduced sidelobes amplitude(less resolution)

For a complete list see (Harris, IEEE Proc.,1978)

Frequency resolution is mainly affected by TR

SMpA 2001 n. 20

Effect of Noise

We have a measured signal and we want toestimate its spectrumMeasures are always affected by noiseA typical approach is averagingUsually, we cannot syncronize the acquisitions with the signal and we have different segment phases(this is sure when the signal is stochastic, i.e. always) So, averaging must be done in the frequency domain

T’R

T”R

Results ! !single acquisition average

what we expect

SMpA 2001 n. 21

Power Spectrum

A possible explaination: the random position of TR produced random phases of thesingle FFTs, therefore averaging does not work.A solution: cancel the phase information by using the Power Density Spectrum (PDS)

PDS

ADS

2*( ) ( ) ( ) ( )yy X X Xω ω ω ωΦ = =

* *( ) ( ) ( ) ( )xy yxX Yω ω ω ωΦ = = Φ

( ) ( ) ( )

zz xx yy xy yx

z t x t y t= +Φ = Φ + Φ + Φ + Φ

2( ) ( ) ( )

( )yy xx

xy xx

Y H X

H

H

ω ω ω

ω

=

Φ = Φ

Φ = Φ

definition and some properties

SMpA 2001 n. 22

Estimation of Harmonic Responce W(jω)

$ ( ) ( ) ( )( )

W j W j R jU j

ω ω ωω

= +

We need W(jω)butwe cannot gain access to z

The ratio of transformsgives us a wrong result:

Suppose to have L recordsand to use averaging, the error would be

Useful relations

Φ Φ Φ

Φ Φ Φyy uu rr

uy uu ur

W

W

= +

= +

2

u yW

r

z

eR j

L

ii

L

2

2

1( )( )

ωω

= =∑

( ) ( ) ( ) ( )R j Y j W j U jω ω ω ω= −R can be seen as a prediction error:

Let’s try to minimize it

Ipoth: r(t) uncorrelated with u(t)

SMpA 2001 n. 23

W(jω) Estimation - 2

$ ( ) ( )( ) * *

W j Me eR j

L

Y U Me Y U Me

Lj

ii

L

i ij

i ij

ω ωω

ϕϕ ϕ

= → = =− −

− =−∑ ∑2

2

1

multiply both sides by e j− ϕ

∂∂ϕ

ϕ ϕ ϕ ϕ

ϕ ϕ

eL

Y U Me jU Me jU Me Y U Me

Y U Me U U M Y U Me U U M

i ij

ij

ij

i ij

i ij

i i i ij

i i

2

2 2

1 0= − − + − =

− + = − +

− −

−

∑

∑ ∑ ∑ ∑

e je j e je j* * *

* * * *

eY U

Y Uj i i

i i

− = ∑∑

2 ~ *

*ϕthe derivative is zero for

SMpA 2001 n. 24

W(jω) Estimation - 3

∂∂

eM

20=

substitute

~( )W j uy

uuω =

Ψ

Ψ

~( ) ~ ~ *

*

*

*W j MeY U

U U

U W R U

U UWj i i

i i

i i i

i i

ur

uuω ϕ= = =

+= +− ∑

∑∑

∑ΨΨ

Simple substitutions point outthe errorIt is zero for L inf

~ ~ * ~ *

*

*

*MeY U e Y U

U U

Y U

U Uj i i

ji i

i i

i i

i i

−−

=+

=∑∑∑

∑∑

ϕϕ2

2In the same way

Ratio between PSD estimates

SMpA 2001 n. 25

Coherence Function

Substitute the optimal Win e(jω)e

Y U Me Y U Me

L

LY Y W Y U W Y U WW U U

LY Y W Y U

Y U

U UY U W

Y U

U UU U

LY Y

Y U Y U

U U

i ij

i ij

i i i i i i i i

i i i ii i

i ii i

i i

i ii i

i ii i i i

i iyy

uy

2

2

1

1

1

( )

~ ~ ~ ~

~ ~

* *

* * * * * *

* * **

** *

*

**

** *

*

ωϕ ϕ

=− −

=

= − − + =

= − − +LNMM

OQPP

= −LNMM

OQPP

= −

−∑

∑ ∑∑∑

∑ ∑∑

∑∑

∑∑∑

∑ ∑∑

∑ ΨΨ

Ψuuyy≈ −Ψ 1 2γe j

γ ωω

ω ωγ ω2

220 1( )

( )( ) ( )

( )jj

j jjuy

uu yy= ≤ ≤

Φ

Φ Φ

1 no noise

0 only noise

SMpA 2001 n. 26

Coherence Function - 2

γ 22 2

2≈ =~ ~W Wuu

uu yy

uu

yy

ΨΨ Ψ

ΨΨ

Suppose we have estimated the TF as

Ψ

Ψyy

uuW=

12

2

γ~

Ψ

Ψyy

uuIs this wrong?

Write the coherence as:

Therefore

that is correct only if noise is zero

SMpA 2001 n. 27

Feedback Systems

G

H

u xn

y

d

v We cannot assume thatx(t) is uncorrelated with noises

We want G(jω),sometimes also H(jω)

We have

Y GH N G U D1+ = + −a f ( )

X U D HY U D HGH

N U D GHGH

= − − = − −+

− −+

( ) ( )1 1

X GH HN U D1+ = − + −a f ( )

SMpA 2001 n. 28

Wrong Way to Estimate G(jω)

Φ xy

(similarly)

$G G=

Limit and interesting cases

!$

*G G G H

Hxy

xx

uu dd nn

nn dd uu= =

+ −

+ +

Φ

ΦΦ Φ Φ

Φ Φ Φ2

Φ xx GH

HN H N HN U D H N U D U D U D

1 2+ =

= − − − − − − + − −

a fe j e j e j* * * * * * * *( ) ( )

Φ xx

Φ Φ Φ Φxx nn dd uuGH H1 2 2+ = + +a f

Φ Φ Φ Φxy uu dd nnGH G G H1 2+ = + −a f *

Φnn = 0

Φ Φuu dd= =0 0, $G H= −1

To use the previous estimator we need

n,d,u uncorrelated

SMpA 2001 n. 29

Adding More Measures

ΦuyΦuu

=G

1+ GH

Further useful eqs.

1 + GH( )Φ uv = GHΦuu

Φ

ΦΦ

Φuy

uxuu

uu

GGH

GHG=

+⋅

+=

11a f

G

H

u xn

y

d

v1 + GH( )Φ ux = Φuu

ΦΦ

ΦΦ

uv

uyuu

uu

GHGH

GHG

H=+

⋅+

=1

1a f

We need also the measures of u(t) and v(t),easily accessible.

SMpA 2001 n. 30

Example

10-1 100 101 10210-1

100

101

102correct estimate & naive estimate & computed TF

1

s+10

9000

s +30s+2002

T ransfer Fcn

t

T o Workspace3

y

zuSignal

Generator

2

Gain

Clock

Chirp S igna l

SMpA 2001

Funzioni di Correlazione

DefinizioniProprietà

Stima della risposta impulsiva

SMpA 2001 n. 32

F. di correl

Misura la somiglianza di due segnali per diversi ritardi

ϕfg (τ) = limT→∞

12T

f(t) g(t + τ)dt−T

T

∫

per eliminare i valori medi, si fa la correlazione dei processi centrati

ϕ°fg (τ) = limT→∞

12T

f(t) − f [ ] g(t + τ) − g [ ]dt−T

T

∫ = ϕfg (τ) − f g

SMpA 2001 n. 33

Proprietà della f. di corr.

ϕ τ ϕ τ

ϕ τ ϕ τ

ϕ τ ϕ

ϕ

ϕ τ ϕ ϕ

ω ϕ τ

gf fg

ff ff

ff ff

ff RMS

fg f g

fg fg

V f t

( ) ( )

( ) ( )

( ) ( )

( ) ( )

( ) ( ) ( )

( ) ( )

= −

= −

≤

=

≤

=

0

0

0 0

2

2

di

Φ Fc h

simmetria

un segnale è simile a sé stesso

importantissima

0 50000

0.5

1

-5000 0 50000.1

0

0.1

0.2

0.3

SMpA 2001 n. 34

Stime della f. di correlazione

˜ ϕ ' fg (τ) =1

2Tf(t)g(t + τ)dt

−T

T

∫

˜ ϕ "fg (τ) =1

2T − τf(t)g(t + τ)dt

−T

T

∫

ϕfg(d)(k) =

12N +1

g(i) f(i − k)i=−N

N∑

Biased

Unbiased

Digitale

-5000 0 5000-0.5

0

0.5biased

-5000 0 5000-0.6

-0.4

-0.2

0

0.2

0.4

0.6unbiased

0 1000 2000 3000 4000 5000-1

-0.5

0

0.5

1

SMpA 2001 n. 35

Correlazione: c’è un segnale?

0 1000 2000 3000 4000 5000-20

-10

0

10

20

-5000 0 5000-1

-0.5

0

0.5

1unbiased

x1=sin(0.01:0.01:50);x1=x1+5*randn(size(x1));x2=cos(0.01:0.01:50);

Al centro viene impiegatoil massimo numero di campioni

SMpA 2001 n. 36

Stima della Risposta Impulsiva Discreta

u y

Sistema lineare stabileda cui

Semplice da risolvere in h,ma ci sono i rumori di misura

yt = ht−kukk=0

t∑ = h jut− j

j=0

∞∑

y0 = h0u0y1 = h1u0 + h0u1y2 = h2u0 + h1u1 + h0u2

SMpA 2001 n. 37

In forma Matriciale

in forma di C.L.

per k->inf, la Σ è il guadagnodel sistema

se E uk[ ]= u , ∀k, allora Eyk[ ]= u ⋅ hj

0

k∑

yt = hjut− j

j=0

∞∑ ; ui = 0, se i < 0

ytyt−1

M

y0

⎡

⎣

⎢ ⎢ ⎢ ⎢

⎤

⎦

⎥ ⎥ ⎥ ⎥

=

h0 h1 L ht0 h0 h1

O

0 L 0 h0

⎡

⎣

⎢ ⎢ ⎢ ⎢

⎤

⎦

⎥ ⎥ ⎥ ⎥

utut−1

M

u0

⎡

⎣

⎢ ⎢ ⎢ ⎢

⎤

⎦

⎥ ⎥ ⎥ ⎥

y = H u

se cov uk[ ]= σ2 , ∀k, allora cov r y [ ]= H ⋅σ2I ⋅H T = σ2HHT

SMpA 2001 n. 38

Una relazione fondamentale

u y

ruy ruu

nel tempo continuo

yi+k = hj ui+k− jj=0

∞∑

ruy(k) = limN→∞

12N +1

uii=−N

N∑ h j ui+k− j

j=0

∞∑ =

h jj=0

∞∑ lim

N→∞

12N + 1

ui ui+k− ji=−N

N∑

⎡

⎣ ⎢

⎤

⎦ ⎥ = h j ruu(k − j)

j=0

∞∑

cross-correlazione In/Out

autocorrelazione

ruy (τ) = h(t)ruu (τ − t)dt0∞

∫

SMpA 2001 n. 39

Stima di h(t) dalla correlazione

Applicando un rumore bianco come ingresso per la misura

Che succede se ci sono rumori additivi su y ?

Dove v è un rumore di misura e/o l’uscita dovuta ad altre frazioni di ingresso (scorrelate dal resto e a media nulla)

Da cui

y( j) = ˆ y ( j) + v( j)

ruy(k) ≈1

2N +1ui ˆ y i+k + vi+k( )

i=−N

N∑ = ˆ r uy (k) + ruv(k) = ˆ r uy(k)

ruu (k) = σ2 δ(k) ⇒ ruy (k) = h(k) ⋅σ2

h kr kuy( )

( )=

σ 2

SMpA 2001 n. 40

Lo stesso in frequenza

Le trasformate di auto ecross-correlazione sonolo spettro e il cross-spettro di potenza

E si trova il metodo dei periodogrammi

Ruy (jω) = F[ruy(τ)] =1

2TY(− jω)U(jω)

Ruu (jω) =1

2TU(jω) 2

ruy (τ) = h(t)ruu (τ − t)dt0∞

∫ → R uy(jω) = H(jω) ⋅ Ruu( jω)

H( jω) =R uy(jω)R uu(jω)

Giovanni Ulivi:

usato R invece di \Fi

Giovanni Ulivi:

usato R invece di \Fi

SMpA 2001

Cepstrum

DefinizioniProprietà

Alcuni impieghi

SMpA 2001 n. 42

Definizioni

Cepstrum complesso• funzione complessa del tempo

Cepstrum di potenza (Power cepstrum)• funzione reale del tempo

Esistono le trasformazioni inverseREM: lo spettro di potenza di segnali misurati si stima meglio dello spettro di ampiezza, quindi anche il cepstrum

[ ]{ } { }-1 -1( ) log ( ) log( ( ) ) ( )xC x t X jτ = = ω + ϑ ωF F F

[ ]{ } { }2-1 -1( ) log ( ) log( ( ) )xx xxC x tτ = = Φ ωF F F

SMpA 2001 n. 43

Separazione tra Sorgente e Canale

2( ) ( ) ( ) ( ) ( ) ( )Y j X j H j H jyy xxω = ω ω Φ ω = Φ ω ω

( ) ( ) ( )yy xx hhC C Cτ = τ + τ

Gli impieghi del Cepstrum sono legati alla sua natura logaritmica

Un impiego è nel separare il pitch delle vocali (periodo abbastanza lungo) dalle caratteristiche del tratto vocale (breve durata della risposta impulsiva)

SMpA 2001 n. 44

Rimozione di un Eco

l1

l2

Il segnale riflesso è

• ritardato τ = (l2– l1)/c

• attenuato l1 / l2• modificato dalla superficie H(ω)

x(t) y(t)

1

2( ) ( ) 1 ( ) jlY j X j H j e

l− ωτ⎡ ⎤

ω = ω + ω⎢ ⎥⎣ ⎦

*2 2 1 1

2 2log ( ) log ( ) log 1 ( ) log 1 ( )j jl l

Y j X j H j e H j el l

− ωτ − ωτω = ω + + ω + + ω⎡ ⎤ ⎡ ⎤⎢ ⎥ ⎢ ⎥⎣ ⎦ ⎣ ⎦

2log(1 ) / 2x x x+ = − +L

211 1( ) ( ) ( ) ( ) ( )22 2

l lC t C t h t h t h tyy xx l l

= + − τ − − τ ⊗ − τ +⎛ ⎞⎜ ⎟⎝ ⎠

L

segnale originale riflessione

se è un vero eco (τ lungo), si può risalire a H e eliminarlo

approssimando