Classical psd

Fourier Methods of Spectral Estimation C.S.Ramalingam Department of Electrical Engineering IIT Madras

C.S.Ramalingam

Fourier Methods of Spectral Estimation

Outline

Definition of Power Spectrum Deterministic signal example

Power Spectrum of a Random Process The Periodogram Estimator The Averaged Periodogram Blackman-Tukey Method

Use of Data Windowing in Spectral Analysis Spectrogram: Speech Signal Example

C.S.Ramalingam

Fourier Methods of Spectral Estimation

What is Spectral Analysis?

Spectral analysis is the estimation of the frequency content of a random process By “frequency content” we mean the distribution of power over frequency Also called Power Spectral Density, or simply spectrum

What frequency components are present? What is the intensity of each component?

C.S.Ramalingam

Fourier Methods of Spectral Estimation

The Earliest Spectral Analyzer

All colours are present with equal intensity

C.S.Ramalingam

Fourier Methods of Spectral Estimation

What is Frequency?

Our notion of frequency comes from sin(2πf0 t) and cos(2πf0 t) both are called “sinusoids”

Frequency ≡ Sinusoidal Frequency: f0 cycles/sec (Hz) exp(j2πf0 t) is the basis function needed for representing a component with frequency f0 for an arbitrary frequency component, it becomes exp(j2πft)

As f varies from −∞ to ∞, we get the Fourier basis set! f is sometimes called “Fourier frequency”

C.S.Ramalingam

Fourier Methods of Spectral Estimation

Spectral Analysis ≡ Expansion Using Fourier Basis

Spectral analysis is nothing but expanding a signal x(t) using the Fourier basis Z ∞ X (f ) = x(t) exp(−j2πft) dt ← inner product! −∞

“X (f ) is the continuous-time Fourier transform of x(t)”

|X (f1 )| large ⇒ dominant frequency component at f = f1 |X (f1 )| = 0 ⇒ no frequency component at f = f1

Plot of |X (f )|2 as a function of f is called the “power spectrum”

C.S.Ramalingam

Fourier Methods of Spectral Estimation

Deterministic Signal Example Gated sinusoid with f0 = 1 kHz 1 0 −1 −2

0

2

4

6 8 Time (ms)

10

12

14

16

Magnitude

60 40 20

Log Magnitude

0 −4000

−3000

−2000

−1000

0 1000 Frequency (Hz)

2000

3000

4000

−3000

−2000

−1000

0 1000 Frequency (Hz)

2000

3000

4000

40 20 0 −20 −4000

C.S.Ramalingam

Fourier Methods of Spectral Estimation

What About PSD of a Random Process? “Spectral analysis is the estimation of the frequency content of a random process” An ensemble of sample waveforms constitute a random process Z ∞ ? X (f ) = x(t) exp(−j2πft) dt −∞

Does it exist? Even if it does, is it meaningful? We’ll focus on discrete-time random processes, i.e., ensemble of x[n], where n ∈ Z

C.S.Ramalingam

Fourier Methods of Spectral Estimation

PSD of a WSS Random Process

Let x[n] be a complex wide-sense stationary process Its autocorrelation sequence (ACS) is defined as rxx [k] = E{x ∗ [n] x[n + k]} Wiener-Khinchine Theorem: Pxx (f ) =

∞ X

rxx [k] exp(−j2πfk)

−∞

−

1 1 ≤f ≤ 2 2

DTFT

That is, ACS ←→ PSD

C.S.Ramalingam

Fourier Methods of Spectral Estimation

An Alternative Definition for PSD

If the ACS decays sufficiently rapidly, 2  M  X 1 x[n] exp(−j2πfn) Pxx (f ) = lim E  M→∞  2M + 1  

−M

The so-called “Direct Method” is based on the above formula

C.S.Ramalingam

Fourier Methods of Spectral Estimation

Why is the Problem Difficult?

ACS is not available Finite number of samples from one realization We are only given x[0], x[1], . . . , x[N − 1]

No “best” spectral estimator exists Many practical signals, such as speech, are non-stationary Pxx (f ) obtained from given data is a random variable Bias versus Variance trade-off

C.S.Ramalingam

Fourier Methods of Spectral Estimation

The Periodogram Estimator Recall 2  M  X 1 x[n] exp(−j2πfn) Pxx (f ) = lim E  M→∞  2M + 1  

−M

In practice we drop lim because data are finite

M→∞

the expectation operator E since we have only one realization

The Periodogram estimator is defined as ˆPER (f ) = 1 P N def

2 N−1 X x[n] exp(−j2πfn) n=0

“Direct Method”, since it deals with the data directly C.S.Ramalingam

Fourier Methods of Spectral Estimation

Example: Two Sine Waves + Noise x[n] =

√

10 exp(j 2π 0.15n) +

√

20 exp(j 2π 0.2n) + z[n]

z[n] ∼ complex N (0, 1), N = 20 50 40

Magnitude (dB)

30 20 10 0 −10 −20 −30 −0.5

−0.4

−0.3

−0.2

−0.1

0

0.1

0.2

0.3

0.4

0.5

Frequency C.S.Ramalingam

Fourier Methods of Spectral Estimation

Periodogram is a Biased Estimator For Finite Data For finite N, periodogram is a biased estimator Bias is the difference between the true and expected values

50 averaged noiseless

40

Magnitude (dB)

30 20

N = 20

10 0 −10 −20 −30 −0.5

−0.4

−0.3

−0.2

−0.1

0

0.1

0.2

0.3

0.4

0.5

Frequency C.S.Ramalingam

Fourier Methods of Spectral Estimation

Periodogram: Bias Decreases With Increasing N If data length is increased, bias decreases: n o ˆxx (f ) = Pxx (f ) lim E P N→∞

60 50

Magnitude (dB)

40 30 20

N = 100

10 0 −10 −20 −30 −0.5

−0.4

−0.3

−0.2

−0.1

0

0.1

0.2

0.3

0.4

0.5

Frequency C.S.Ramalingam

Fourier Methods of Spectral Estimation

The More (Samples) the Merrier?

For most estimators, bias and variance decrease with increasing N An estimator is said to be consistent if   lim Pr θˆ − θ >  = 0 N→∞

where θˆ is the estimate of θ This implies that, as N → ∞, bias → 0 variance → 0

C.S.Ramalingam

Fourier Methods of Spectral Estimation

Is the Periodogram Consistent?

Consider white noise sequence for various N True Pxx (f ) = constant If the Periodogram estimator were consistent, ˆxx (f ) → constant as N increases P Consider noise sequences of length 32, 64, 128, and 256 N = 32; x = randn(N,1);

% white noise sequence % of length 32

ˆxx (f ) tend to a constant as N increases? Does P

C.S.Ramalingam

Fourier Methods of Spectral Estimation

White Noise Example PSD of White Noise 40

N=32

0 −0.5 40

0

0.5

0

0.5

0

0.5

0 −0.5 40 0 −0.5 40

N=256

0 −0.5

−0.4

−0.3

−0.2

−0.1

0

0.1

0.2

0.3

0.4

0.5

Frequency

As N increases, variance of the estimate does not decrease Periodogram is an inconsistent estimator C.S.Ramalingam

Fourier Methods of Spectral Estimation

What Went Wrong?

“In practice we drop lim because data are finite

M→∞

the expectation operator E since we have only one realization”

For white noise, increasing the data length did not help What can be done to capture the benefits of E{·} ?

C.S.Ramalingam

Fourier Methods of Spectral Estimation

Averaging: The Poor Man’s Expectation Operator Expectation operator can be approximated by averaging Averaged Periodogram: M 1 X ˆ (m) ˆ PPER (f ) PAVPER (f ) = M m=1

(m) ˆPER where P (f ) is periodogram of m-th segment of length N

For independent data records n n o o ˆPER (f ) ˆAVPER (f ) = 1 var P var P M

C.S.Ramalingam

Fourier Methods of Spectral Estimation

Averaged Periodogram of White Noise Result of averaging 8 periodograms Averaged Periodogram 40 35 30

Magnitude (dB)

25 20 15 10 5 0 −5 −10 −0.5

−0.25

0 Frequency

C.S.Ramalingam

0.25

0.5

Fourier Methods of Spectral Estimation

Variance Decreases, But Bias Increases! Two Sines + Noise example Averaged Periodogram for Two Sines + Noise 20

10

N=256, M=1 N=64, M=4 N=16, M=16

Magnitude (dB)

0

−10

−20

−30

−40

−50 −0.5

−0.25

0 Frequency

C.S.Ramalingam

0.25

0.5

Fourier Methods of Spectral Estimation

Welch’s Method Overlapping blocks by 50% Reduces variance without worsening bias Welch’s Method 15

10

Block Length = 64 No. of blocks = 7 Overlap = 50%

Magnitude (dB)

5

0

−5

−10

−15

−20

−25 −0.5

−0.25

0 Frequency

C.S.Ramalingam

0.25

0.5

Fourier Methods of Spectral Estimation

Why Did The Periodogram Fail? Periodogram was defined as 2 P ˆPER (f ) = 1 N−1 x[n] exp(−j2πfn) P n=0 N Equivalent to ˆPER (f ) = P

N−1 X

ˆrxx [k] exp(−j2πfk)

−(N−1)

where

 N−1−k X    1 x ∗ [n] x[n + k] k = 0, 1, . . . , N − 1 N ˆrxx [k] = n=0    ˆr ∗ [−k] k = −(N − 1), . . . , −1 xx Note that ˆrxx [N − 1] = x ∗ [0]x[N − 1]/N No averaging ⇒ estimate with high variance! C.S.Ramalingam

Fourier Methods of Spectral Estimation

Blackman-Tukey Method Recall Pxx (f ) =

∞ X

rxx [k] exp(−j2πfk)

−∞

−

1 1 ≤f ≤ 2 2

In practice: (a) replace rxx [k] by estimate ˆrxx [k], (b) truncate the summation, and (c) apply “lag window” M X ˆ PBT f ) = w [k] ˆrxx [k] exp(−j2πfk) −M

where 0 ≤ w [k] ≤ w [0] = 1 w [−k] = w [k]

w [k] = 0 for |k| > M W (f ) ≥ 0

“Indirect Method”, since it does not deal with the data directly C.S.Ramalingam

Fourier Methods of Spectral Estimation

Example: Two Sine Waves + Noise Data length N = 100, Correlation Lag M = 10 60 50

Magnitude (dB)

40 30 20 10 0 −10 −20 −30 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 Frequency

C.S.Ramalingam

0.2

0.3

0.4

0.5

Fourier Methods of Spectral Estimation

Periodogram Vs. Blackman-Tukey Periodogram Magnitude (dB)

60 40 20 0 −20 −0.5

−0.25

0

0.25

0.5

0.25

0.5

Magnitude (dB)

Blackman−Tukey 60 40 20 0 −20 −0.5

−0.25

0 Frequency

Blackman-Tukey method: reduction in variance comes at the expense of increased bias Speech Analysis C.S.Ramalingam

Fourier Methods of Spectral Estimation

Data Windowing in Spectral Analysis

Useful for data containing sinusoids + noise Sidelobes of a stronger sinusoid may mask the main lobe of a nearby weak sinusoid We multiply x[n] by data window w [n] before computing periodogram Weaker sinusoid becomes more visible Main lobe of each sinusoid broadens: two close peaks may merge into one

C.S.Ramalingam

Fourier Methods of Spectral Estimation

Example: How Many Sine Waves Are There? How Many Sinusoids Are There? 60

Magnitude (dB)

40 20 0 −20 −40 −60 0.1

0.12

0.14

0.16 Frequency

C.S.Ramalingam

0.18

0.2

0.22

Fourier Methods of Spectral Estimation

Example: Three Sine Waves Three Sinusoids: Rectangular Window 60

Magnitude (dB)

40 20 0 −20 −40 0.15 −60 0.1

0.12

0.14

0.157 0.16 Frequency

C.S.Ramalingam

0.18

0.2

0.22

Fourier Methods of Spectral Estimation

Example: Three Sine Waves Three Sinusoids: Hanning Window 60

Magnitude (dB)

40 20 0 −20 −40 0.15 −60 0.1

0.12

0.14

0.157 0.16 Frequency

C.S.Ramalingam

0.18

0.2

0.22

Fourier Methods of Spectral Estimation

Commonly Used Windows

Name

w [k]

Rectangular

1 |k| M

Bartlett

1−

Hanning

0.5 + 0.5 cos

Hamming

0.54 + 0.46 cos

πk M πk M

Fourier transform sin πf (2M + 1) WR (f ) =  sin πf 1 sin πfM 2 M sin πf
1 0.25 WR f − 2M
+ 0.5 WR (f ) + 1 0.25 WR f + 2M
1 0.23 WR f − 2M
+ 0.54 WR (f ) + 1 0.23 WR f + 2M

w [k] = 0 for |k| > M

C.S.Ramalingam

Fourier Methods of Spectral Estimation

Hamming Vs. Hanning Magnitude (dB)

Fourier Transforms of Hamming and Hanning Windows 0 −50 −100 −0.5

−0.25

0 Frequency

0.25

Magnitude (dB)

0

0.5

Hamming Hanning

−20 −40 −60 −0.05

−0.025

0 Frequency

C.S.Ramalingam

0.025

0.05

Fourier Methods of Spectral Estimation

Three Sine Waves Rectangular Vs. Hamming Vs. Hanning 60

Magnitude (dB)

40 20 0 −20 −40 0.15 −60 0.1

0.12

0.14

0.157 0.16 Frequency

C.S.Ramalingam

0.18

0.2

0.22

Fourier Methods of Spectral Estimation

How Can We Analyze Non-Stationary Signals? Consider a “linear chirp”, i.e., a signal whose frequency increases linearly from f1 Hz to f2 Hz over a time interval T What is its magnitude spectrum?

Amplitude

1

0

−1 0

0.01

0.02

0.03

0.04

0.05

400 600 Frequency

800

1000

Time Magnitdue (dB)

60 40 20 0 0

200

C.S.Ramalingam

Fourier Methods of Spectral Estimation

Need a More Useful Representation

In Fourier analysis, even if a signal is non-stationary, it is still represented using stationary sinusoids An unsatisfactory approach

Power spectrum is identical to x(−t), whose frequency decreases from f2 to f1 x(t) and x(−t) differ only in the phase of the Fourier transform

What we really want to know is how frequency varies with time Can it still be called “frequency” ?

C.S.Ramalingam

Fourier Methods of Spectral Estimation

Spectrogram

Plot of power spectrum of short blocks of a signal as a function of time Over each short block, signal is considered to be stationary Speech is a classic example of a commonly occurring non-stationary signal Voiced sounds: /a/, /e/, /i/, /o/, /u/ (quasi-periodic) Unvoiced sounds: /s/, /sh/, /f/ (noise-like) Plosives: /p/, /t/, /k/ (transient sounds)

C.S.Ramalingam

Fourier Methods of Spectral Estimation

Spectrogram of Linear Chirp 600 500

Frequency

400 300 200 100 0 0

0.1

0.2

0.3

0.4

Time

C.S.Ramalingam

Fourier Methods of Spectral Estimation

Non-stationarity in Speech Signal 0.2 0 /k/ −0.2 1.16 1.17 1

1.18

1.19

1.2

1.21

1.22

1.23

1.24

0 /ow/ −1 1.24 1.26 0.05

1.28

1.3

1.32

1.34

1.36

1.38

0 /s/ −0.05 1.75

1.8

1.85 Time

C.S.Ramalingam

Fourier Methods of Spectral Estimation

1.9

Application to Speech Analysis Should We Chase Those Cowboys? Amplitude

1 0.5 0 −0.5 −1 0

0.25

0.5

0.75

1

1.25

1.5

1.75

2

0.25

0.5

0.75

1 Time

1.25

1.5

1.75

2

Frequency

5000 4000 3000 2000 1000 0 0

C.S.Ramalingam

Fourier Methods of Spectral Estimation

Application to Speech Analysis Should We Chase Those Cowboys? Amplitude

1 0.5 0 −0.5 −1 0

0.25

0.5

0.75

1

1.25

1.5

1.75

2

0.25

0.5

0.75

1 Time

1.25

1.5

1.75

2

Frequency

5000 4000 3000 2000 1000 0 0

C.S.Ramalingam

Fourier Methods of Spectral Estimation

Application to Speech Analysis Should We Chase Those Cowboys? Amplitude

1 0.5 0 −0.5 −1 0

0.25

0.5

0.75

1

1.25

1.5

1.75

2

0.25

0.5

0.75

1 Time

1.25

1.5

1.75

2

Frequency

5000 4000 3000 2000 1000 0 0

C.S.Ramalingam

Fourier Methods of Spectral Estimation

Application to Speech Analysis Should We Chase Those Cowboys? Amplitude

1 0.5 0 −0.5 −1 0

0.25

0.5

0.75

1

0.25

0.5

0.75

1 Time

1.25

1.5

1.75

Frequency

5000 4000 3000 2000 1000 0 0

C.S.Ramalingam

1.25

1.5

1.75

Fourier Methods of Spectral Estimation

2

Summary Definition of Power Spectrum Deterministic signal example

Power Spectrum of a Random Process The Periodogram Estimator The Averaged Periodogram Bias versus Variance

Blackman-Tukey Method

Use of Data Windowing in Spectral Analysis Spectrogram: Speech Signal Example

C.S.Ramalingam

Fourier Methods of Spectral Estimation