...

Ch.8.2 slides

by user

on
Category: Documents
1

views

Report

Comments

Transcript

Ch.8.2 slides
Digtal Signal Processing (et 4235)
Non-Parametric Spectrum Estimation (chapter 8.2)
1
8 non-parametric spectrum est. October 10, 2011
Spectral analysis
. We would like to show its “spectrum”.
is very long (e.g. an hour of music...) it does not make sense to compute the
If
Given a signal
Fourier transform—the signal is nonstationary.
Typically we break up the signal into shorter segments and compute the Fourier
transform of each segment (time-frequency analysis).


otherwise

How does this relate to the Fourier transform of the original signal?
2
8 non-parametric spectrum est. October 10, 2011
Spectral analysis
Describe the domain restriction by a multiplication with a rectangular window:


otherwise

Multiplication in time domain is convolution in frequency domain:
where
X
This is the Dirichlet kernel (generalized sinc function)
3
8 non-parametric spectrum est. October 10, 2011
Spectral analysis
Rectangular window
1
N=19
0.8
0.6
W(ω)
0.4
0.2
0
−0.2
−2pi
0
ω
is at
pi
2pi
. The highest sidelobe occurs for
The first zero crossing of
−pi
−0.4
and is at
dB of the main lobe.
The window gives distortions: smearing (resulting in loss of resolution) and sidelobe interference (resulting in confusion)
4
8 non-parametric spectrum est. October 10, 2011
Spectral analysis
Example:
cos(ω0 n) w[n]
N=19
Re(X^(ω))
10
5
0
−5
−pi
−pi/2
0
ω
5
pi/2
pi
8 non-parametric spectrum est. October 10, 2011
Spectral analysis
Spectrum via DFT
.
, with
In practice, we use the DFT (or FFT). This corresponds to sampling
) we can also take (
is periodic (period
Because
)
corresponding with
6
8 non-parametric spectrum est. October 10, 2011
Spectral analysis
Example
. What is its frequency?
Consider the following signal
x[n]
1
cos(ω n)
0
0
−1
0
|X(ω)|
10
5
10
n
15
20
−pi/2
0
ω
pi/2
pi
N=19
L=19
5
0
−pi
7
8 non-parametric spectrum est. October 10, 2011
Spectral analysis
additional zeros to length . (The amount of information will not change)
1
cos(ω0 n)
x^[n]
0.5
0
−0.5
−1
0
|X(ω)|
10
5
10
15
20
25
n
30
35
40
45
50
N=19
L=50
5
0
−pi
−pi/2
0
ω
pi/2
pi
Spectral “leak” (distortion) by the side lobes of the window
8
8 non-parametric spectrum est. October 10, 2011
with
To get a smoother spectrum, we can use zero padding, i.e., extend
Spectral analysis
Sum of two frequencies
2
cos(ω n) + cos(ω n)
x^[n]
1
1
2
0
−1
−2
0
20
|X(ω)|
15
10
20
30
n
40
50
60
omega1=0.2 ⋅ 2π
omega =0.22 ⋅ 2π
N=20
L=500
1/N=0.05
2
10
5
0
−pi
0
ω
is at
pi/2
pi
. This determines the resolution.
The first zero crossing of
−pi/2
(or
Frequencies with differences smaller than
) cannot be
resolved.
9
8 non-parametric spectrum est. October 10, 2011
Spectral analysis
For better resolution: enlarge
(number of informative samples):
2
cos(ω n) + cos(ω n)
1
2
^
x [n]
1
0
−1
−2
0
20
^
|X (ω)|
15
10
20
N=40
L=500
1/N=0.025
30
n
40
omega =0.2 ⋅ 2π
1
omega =0.22 ⋅ 2π
50
60
Rectangular windowed
2
10
5
0
−pi
−pi/2
0
ω
pi/2
pi
The resolution is further limited by the height of the side lobes (confusion).
10
8 non-parametric spectrum est. October 10, 2011
Spectral analysis
Side lobe reduction
Other windows can be designed to obtain lower side lobes, at the expense of a
wider main lobe. E.g., Bartlett, Hann, Hamming, Blackman, Kaiser, etc.
1
Bartlett window
0.5
0
0
2
4
6
8
10
12
14
16
18
n
|W(ω)|
1
0.5
0
−pi
N=19
Bartlett window
(rectangular window)
−pi/2
0
ω
11
pi/2
pi
8 non-parametric spectrum est. October 10, 2011
Spectral analysis
Bartlett window
2
cos(ω n) + cos(ω n)
x [n]
1
1
2
^
0
−1
−2
0
8
^
|X (ω)|
6
20
40
N=40
L=500
1/N=0.025
60
n
80
omega =0.2 ⋅ 2π
1
omega =0.22 ⋅ 2π
100
120
Bartlett windowed
2
4
2
0
−pi
−pi/2
0
ω
pi/2
pi
Lower side lobes at the expense of some resolution
12
8 non-parametric spectrum est. October 10, 2011
Spectral analysis
Hamming window
1
Hamming window
0.5
0
0
2
4
6
8
10
12
14
16
18
n
|W(ω)|
1
0.5
0
−pi
N=19
Hamming window
(rectangular window)
−pi/2
0
ω
pi/2
pi
13
8 non-parametric spectrum est. October 10, 2011
Spectral analysis
Hamming window
2
cos(ω n) + cos(ω n)
x [n]
1
1
2
^
0
−1
−2
0
8
^
|X (ω)|
6
20
40
N=40
L=500
1/N=0.025
60
n
80
omega =0.2 ⋅ 2π
1
omega =0.22 ⋅ 2π
100
120
Hamming windowed
2
4
2
0
−pi
−pi/2
0
ω
14
pi/2
pi
8 non-parametric spectrum est. October 10, 2011
Spectral analysis
Wide Sense Stationary random signal
)
Example: communication signal (binary phase shift keying, with random bits
bpsk signaal
2
x[n]
1
0
−1
−2
0
10
20
30
40
50
n
60
70
80
90
100
40
|X(ω)|
30
N=100
20
10
0
−pi
−pi/2
0
ω
pi
gives a chaotic picture. Do the
Applying the DFT directly on the signal
pi/2
peaks correspond to certain frequencies?
15
8 non-parametric spectrum est. October 10, 2011
Spectral analysis
bpsk signaal
2
x[n]
1
0
−1
−2
0
100
200
300
400
500
n
600
700
800
900
1000
200
|X(ω)|
150
N=1000
100
50
0
−pi
gives
A larger
−pi/2
0
ω
pi/2
pi
points in the frequency domain, but the spectrum does not
converge. We need to average.
16
8 non-parametric spectrum est. October 10, 2011
Spectral analysis
2
2
2
2
2
1
1
1
1
1
0
0
0
0
0
−1
−1
−1
−1
−1
−2
0
50
n
100
−2
100
20
40
15
30
10
20
5
10
0
−pi −pi/2
0 pi/2 pi
ω
0
−pi −pi/2
150
n
200
−2
200
250
n
300
−2
300
30
30
20
20
350
n
400
−2
400
450
n
500
40
30
20
10
0 pi/2 pi
ω
10
0
−pi −pi/2
0 pi/2 pi
ω
10
0
−pi −pi/2
(here
Separate the signal into segments, each of length
0
−pi −pi/2
0 pi/2 pi
ω
0 pi/2 pi
ω
). Do a
windowed DFT on each segment. Compute the amplitude-square, and average.
This gives the averaged periodogram (a.k.a. Welch method).
17
8 non-parametric spectrum est. October 10, 2011
Spectral analysis
Averaged periodogram
1
0.8
N=10000
L=100
T =6
s
bpsk spectrum
single pulse (dirichlet)
0.6
Hamming window
0.4
0.2
0
−pi
−pi/2
0
pi/2
pi
The spectrum is seen to converge to that of the time-domain pulse (a rectangular
pulse
Dirichlet-kernel)
18
8 non-parametric spectrum est. October 10, 2011
Random process
After this introduction, we now follow the theory in Hayes, Ch. 8.2.
Given a random sequence, how can we meaningfully estimate its spectrum?
Definition
X
For an autocorrelation ergodic process:
X
for all , to compute the power spectrum,
determine
Hence: given
compute its Fourier transform
)?
19
But what if we have data only over a finite interval (
8 non-parametric spectrum est. October 10, 2011
Periodogram
),
For finite data (
Use a finite sum,
X
Exclude data outside the interval (equivalent to setting it to zero):
X
(as we see later, this is a biased estimate)
for
.
, and set
Use also conjugate symmetry:
Then the periodogram is defined as
X
20
8 non-parametric spectrum est. October 10, 2011
Periodogram
Expressed in terms of the original data: define


otherwise
is a rectangular window. Then
, where
so that

X
Taking the Fourier transform gives
where
X
X
.
21
Thus, the periodogram can be computed from the DFT of
8 non-parametric spectrum est. October 10, 2011
Periodogram
Performance of the periodogram
Ideally, a (spectrum) estimate is consistent:
asymptotically unbiased:
Var
:
variance goes to 0 as
Bias of the periodogram

is a biased estimate.
22
8 non-parametric spectrum est. October 10, 2011
is a Bartlett window:
.
where


for
Therefore,
Thus
X
and
X
Periodogram
It follows for the periodogram that
so that
where
is the Fourier transform of the Bartlett window,
times the square of the Dirichet window we saw earlier;
(Note: it is
)
converges to an impulse,
,
The periodogram is biased. But as
and the periodogram is asympotically unbiased:
23
8 non-parametric spectrum est. October 10, 2011
X
X
X
Periodogram
Effect of the convolution by
Spectral smoothing (smearing): the response to a signal consisting of a single
, with a bandwidth of approx.
is not an impulse, but
frequency
.
, one might think there are
at
Confusion: due to sidelobes of
additional frequencies.
This limits the resolution: ability to resolve closely spaced sinusoids. The resolution
is defined by the half-power width of the main lobe,
Res
24
8 non-parametric spectrum est. October 10, 2011
Periodogram
Variance of the periodogram
The computation of the variance is involved. One finds
Var
(exact for white Gaussian noise)
. the variance does not go to 0 as
This is not a function of
and the
periodogram is not a consistent estimate of the power spectrum.
bpsk signaal
2
x[n]
1
0
−1
−2
0
100
200
300
400
500
n
600
700
800
900
1000
200
|X(ω)|
150
N=1000
100
50
0
−pi
−pi/2
0
ω
25
pi/2
pi
8 non-parametric spectrum est. October 10, 2011
Periodogram
Alternative derivation of the same
Describe the domain restriction by a multiplication with a rectangular window:


otherwise

Multiplication in time domain is convolution in frequency domain:
where
X
This is the Dirichlet kernel (generalized sinc function)
26
8 non-parametric spectrum est. October 10, 2011
Modified Periodogram
Recall
for which (derivation see Hayes p.408)
X
(note that this is consistent to the earlier result
)
where
, e.g. a Bartlett window,
Instead of a rectangular window, try other windows
X
is a suitable normalization constant (energy of the window),
X
This is called the modified periodogram.
27
8 non-parametric spectrum est. October 10, 2011
Modified Periodogram
Bias
A derivation similar as before gives
.
is the Fourier transform of
where
Since (by Parseval)
Z
X
it follows that
Z
will converge to a unit impulse.
For
so that
, the modified periodogram is asympotically unbiased.
28
8 non-parametric spectrum est. October 10, 2011
Modified Periodogram
Variance
The window does not change the situation,
Var
Thus, the modified periodogram is not a consistent estimate of the power spectrum.
Effect of the window
The window gives a way to trade-off the spectral resolution (width of main lobe) with
the confusion (sidelobe amplitude).
Sidelobe (dB)
3dB BW
Bartlett
-27
Hann
-32
Hamming
-43
Blackman
-58
29
8 non-parametric spectrum est. October 10, 2011
-13
Rectangular
Periodogram averaging
. Recall
Our aim is to find a consistent estimate for
. This can
It is sufficient to find a consistent estimate for the mean,
be done by averaging several periodograms (based on independent realizations
of the same random process).
uncorrelated realizations of
. Compute the periodogram of each realization:
X
Then compute the average:
over the interval
. We have
be
,
Let
X
30
8 non-parametric spectrum est. October 10, 2011
Periodogram averaging
The resulting averaged periodogram estimate is
X X
Bias
based on a rectangular window with
with
samples.
is asympotically unbiased (as
As before, the estimate
).
Variance
Var
Var
as
.
and
is a consistent estimate if both
Thus,
This goes to
go to infinity.
31
8 non-parametric spectrum est. October 10, 2011
Periodogram averaging
Bartlett’s method
nonoverlapping segments
, split it into
,
where
,
of length
For a given sequence
then compute the averaged periodogram,
X X
2
2
2
2
2
1
1
1
1
1
0
0
0
0
0
−1
−1
−1
−1
−1
−2
0
50
n
100
−2
100
20
40
15
30
10
20
5
10
0
−pi −pi/2
0 pi/2 pi
ω
0
−pi −pi/2
150
n
200
−2
200
250
n
300
−2
300
30
30
20
20
10
10
350
n
400
−2
400
450
n
500
40
30
20
0 pi/2 pi
ω
0
−pi −pi/2
0 pi/2 pi
ω
32
0
−pi −pi/2
10
0 pi/2 pi
ω
0
−pi −pi/2
0 pi/2 pi
ω
8 non-parametric spectrum est. October 10, 2011
Periodogram averaging: Bartlett’s method
Bias
X X
Precisely as before,
based on
where
samples. The estimate is asympotically unbiased.
The resolution is now determined by :
Res
which is
times worse than for the periodogram.
Variance
Var
Var
and
is a consistent estimate if both
Thus,
go to infinity.
33
8 non-parametric spectrum est. October 10, 2011
Periodogram averaging
Welch method: averaged modified periodograms
Two modifications:
Use overlapping segments (each of length , but spaced by
),
The total number of data points that is used is
Often, 50% overlap is used:
.
instead of simple rectangular windowing
Use data windows
The resulting Welch’s estimate is given by
34
X
X X
8 non-parametric spectrum est. October 10, 2011
Periodogram averaging: Welch Method
Bias
(length )
is the Fourier transform of the window
where
Variance
)
Var
One shows that, for 50% overlap (some correlation among the sequences
worse than for Bartlett’s estimate, for the same
This is a factor
. But, twice
doubles), so in fact the variance is a factor
as many segments can be averaged (
better:
Var
35
8 non-parametric spectrum est. October 10, 2011
Fly UP