Chapter 3: Parametric spectral estimation
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 54 / 82
3. Parametric spectral estimation Introduction
Motivation
Non-parametric methods
quite general: applicable to any kind of signals
neglect the specific properties of the signals under analysis
Parametric methods
take into account the signal properties ! signal model, defined by model parameters
estimate the model parameters
estimate the PSD from the estimated model parameters
Benefits
typically reduced number of model parameters
! allows sample size reduction for equivalent performance
improved performance (resolution, variance)...
... if the assumed model is correct.
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 55 / 82
3. Parametric spectral estimation Introduction
General concept
1 The signal model is built from prior information about the process ! model parameters
2 The model parameters are estimated from the observed process
3 The PSD is constructed from the estimated model parameters.
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 56 / 82
3. Parametric spectral estimation Signal modeling
Properties of the PSD
P1) The PSD is defined as the Fourier transform of the ACS
S(!) = F{r(k)} =
+1X
k=1
r(k)e|!k
We can also define the PSD in the z-domain as:
S(z) = Z{r(k)} =
+1X
k=1
r(k)zk with S(!) = S(z)

z=e|!
P2) Since the ACS is conjugate symmetric, i.e., r(k) = r⇤(k), the PSD is real valued:
S(!) = S⇤(!) and S(z) = S⇤(1/z⇤)
P3) If the process is real valued, the PSD is an even function:
S(!) = S(!) and S(z) = S⇤(z⇤)
P4) The PSD is nonnegative:
S(!) 0, 8! 2 R
P5) Total power:
PX = E{|x(n)|2} = r(0) = 1
2⇡
Z ⇡
⇡
S(!)d!
! S(!) represents the power density of the process in the frequency domain (= PSD).
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 57 / 82
3. Parametric spectral estimation Signal modeling
Filtering random processes
WSS random process X (n) with mean µX , ACS rX (k) and PSD SX (!) is filtered with a linear
time-invariant (LTI) filter with transfer function H(z):
Q: What is the mean, ACS and PSD of the output process Y (n)?
We can prove (see tutorial) that
µY = µXH(e
|0)
rY (k) = rX (k) ⇤ h(k) ⇤ h⇤(k)
SY (!) = SX (!)|H(!)|2
SY (z) = SX (z)H(z)H
⇤(1/z⇤)
Special case: white noise input with variance 2X
SY (!) =
2
X |H(!)|2
SY (z) =
2
XH(z)H
⇤(1/z⇤)
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 58 / 82
3. Parametric spectral estimation Signal modeling
Spectral factorization theorem
Let S(!) be the PSD of a WSS random process X (n).
If S(!) is a continuous function of !, then it can be factored as
S(z) = 2"H(z)H
⇤(1/z⇤)
where H(z) is a causal, stable, minimum phase LTI filter, i.e., with no poles or zeros outside the
unit circle.
Any process that admits this factorization is called a regular process.
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 59 / 82
3. Parametric spectral estimation Signal modeling
Properties of regular processes
They can be realized as the output of a causal, stable filter driven by white noise with
variance 2" ! innovations representation of the process.
If filtered with 1/H(z), the output is a white noise with variance 2" ! whitening
I 1/H(z): whitening filter
I white noise output "(n): innovation process
I 2": innovation variance or modeling error.
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 60 / 82
3. Parametric spectral estimation Signal modeling
Spectral factorization (cont’d)
Wold decomposition theorem
Any WSS random process may be decomposed into the sum of two orthogonal processes:
X (n) = Xr (n) + Xp(n)
with E{Xr (n)X⇤p (n k)} = 0, 8k 2 Z, where
Xr (n): regular process
Xp(n): predictable process.
A predictable process can be predicted without error from a linear combination of its previous
values:
Xp(n) =
1X
k=1
a(k)Xp(n k).
Its PSD consists of impulses:
SXp (!) =
KX
k=1
↵k(! !k ).
Corollary: the general form of the PSD of a WSS process X (n) is given by
SX (!) = SXr (!) +
KX
k=1
↵k(! !k ).
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 61 / 82
3. Parametric spectral estimation Signal modeling for rational spectra
Stochastic process models for rational spectra
Important special case: S(z) is a rational function
S(z) =
N(z)
D(z)
=
"
spectral factorization
2"H(z)H
⇤(1/z⇤) = 2"
B(z)B⇤(1/z⇤)
A(z)A⇤(1/z⇤)
R1 < |z| < R2
with
H(z)
def
=
B(z)
A(z)
=
Pq
k=0 bkz
k
1 +
Pp
k=1 akz
k |z| > R1. (12)
By the spectral factorization theorem, H(z) and 1/H(z) are causal, stable, minimum phase
! B(z) and A(z) have all their roots inside the unit circle.
Power spectral density
S(!) = 2"
|B(!)|2
|A(!)|2 =
2
"
Pq
k=0 bke
|!k 21 +Ppk=1 ake|!k 2
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 62 / 82
3. Parametric spectral estimation Signal modeling for rational spectra
Equivalent time series representation
X (z) = H(z)"(z) =
B(z)
A(z)
"(z)
A(z)X (z) = B(z)"(z) ) Z1{A(z)X (z)} = Z1{B(z)"(z)}
X (n)+
pX
k=1
akX (nk) =
qX
k=0
bk"(nk) (13)
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 63 / 82
3. Parametric spectral estimation Signal modeling for rational spectra
Autoregressive moving average (ARMA) processes
If p > 0 and q > 0
! autoregressive moving average (ARMA) process of order (p, q): ARMA(p, q)
X (n) +
pX
k=1
akX (n k) =
qX
k=0
bk"(n k) ) S(!) = 2"
Pq
k=0 bke
|!k 21 +Ppk=1 ake|!k 2
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 64 / 82
当 纳 8时 上仙 1→脉冲响应。
3. Parametric spectral estimation Signal modeling for rational spectra
ARMA processes — examples
0 50 100 150 200 250
−5
0
5
0 50 100 150 200 250
−0.5
0
0.5
0 50 100 150 200 250
−5
0
5
Figure: (Top) White noise input.
(Middle) a = [1,0.975, 0.95], b = [1, 1, 1, 1]/4.
(Bottom) a = [1,0.975, 0.95], b = [1, 1, . . . , 1]| {z }
100 entries
/100.
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 65 / 82
3. Parametric spectral estimation Signal modeling for rational spectra
Link between ARMA parameters and ACS
Establishing this link is useful for estimating the model parameters.
Multiplying eqn. (13) by X⇤(n k) and taking expectations:
X (n)X⇤(n k) +
pX
m=1
amX (n m)X⇤(n k) =
qX
m=0
bm"(n m)X⇤(n k)
E{X (n)X⇤(n k)}+
pX
m=1
amE{X (n m)X⇤(n k)} =
qX
m=0
bmE{"(n m)X⇤(n k)}.
Hence:
r(k) +
pX
m=1
amr(k m) =
qX
m=0
bmr"X (k m).
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 66 / 82
3. Parametric spectral estimation Signal modeling for rational spectra
To compute the cross-correlation r"X (k), we note that:
The innovation processe is white: r"(k) = 2"(k).
The model output is given by the convolution of the filter impulse response
h(n)
def
= Z1{H(z)} and the filter input:
X (n) = h(n) ⇤ "(n) =
+1X
m=1
h(m)"(n m)
H(z) is causal (spectral factorization theorem), i.e., h(n) = 0, n < 0, and then:
+1X
m=1
h(m)"(n m) =
+1X
m=0
h(m)"(n m).
Accordingly:
r"X (k) = E{"(n)X⇤(n k)} = E
n 1X
m=0
h⇤(m)"⇤(n k m)"(n)
o
=
1X
m=0
h⇤(m)E{"⇤(n k m)"(n)}| {z }
r"(k+m)=2"(k+m)
= 2"
1X
m=0
h⇤(m)(k +m) = 2"h⇤(k).
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 67 / 82
3. Parametric spectral estimation Signal modeling for rational spectra
Link between ARMA parameters and ACS (cont’d)
Again, because h(n) is causal (spectral factorization theorem), then
qX
m=0
bmr"X (k m) ="
r"X (k)=
2
"h
⇤(k)
2"
qX
m=0
bmh
⇤(m k) = 2"c(k)
with
c(k)
def
=
qkX
m=0
bm+kh
⇤(m).
Hence, for k 0:
r(k) +
pX
m=1
amr(k m) =
8<:
2
"c(k) 0  k  q
0 k > q.
(14)
For k < 0, we just enforce the conjugate symmetry of the ACS: r(k) = r⇤(k).
These are the Yule-Walker equations, linking the ARMA parameters with the ACS.
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 68 / 82
3. Parametric spectral estimation Signal modeling for rational spectra
Yule-Walker equations — matrix form
Casting the equations for k = 0, 1, . . . , p + q in matrix form, we have:
2666666666664
r(0) r(1) · · · r(p)
r(1) r(0) · · · r(p + 1)
...
...
. . .
...
r(q) r(q 1) · · · r(q p)
r(q + 1) r(q) · · · r(q p + 1)
...
...
. . .
...
r(q + p) r(q + p 1) · · · r(q)
3777777777775
26664
1
a1
...
ap
37775 = 2"
2666666666664
c(0)
c(1)
...
c(q)
0
...
0
3777777777775
(15)
Dicult solution for model coecients of a general ARMA(p, q) process, because c(k)
depends nonlinearly on them (see previous slide) ! nonlinear system of equations.
The system of equations becomes linear for q = 0, since in that case c(k) = |b0|2(k).
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 69 / 82
3. Parametric spectral estimation Signal modeling for rational spectra
Autoregressive (AR) processes
If q = 0 ! ARMA(p, 0) ! autoregressive (AR) process of order p: AR(p)
X (n) +
pX
k=1
akX (n k) = "(n) ) S(!) =
2
"1 +Ppk=1 ake|!k 2 (16)
Also known as all-pole model.
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 70 / 82
3. Parametric spectral estimation Signal modeling for rational spectra
AR processes — examples
0 50 100 150 200 250
−5
0
5
0 50 100 150 200 250
−10
0
10
20
0 50 100 150 200 250
−10
0
10
Figure: (Top) White noise input.
(Middle) AR process with a = [1, 0.1,0.8].
(Bottom) AR process with a = [1,0.975, 0.95].
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 71 / 82
3. Parametric spectral estimation Signal modeling for rational spectra
AR processes — Yule-Walker equations
For AR(p) processes, assuming that b0 = 1, the Yule-Walker equations (14) simplify to:
r(k) +
pX
m=1
amr(k m) = 2"(k) k 0
and their matrix form (15) becomes
26664
r(0) r(1) · · · r(p)
r(1) r(0) · · · r(p + 1)
...
...
. . .
...
r(p) r(p 1) · · · r(0)
37775
26664
1
a1
...
ap
37775 =
26664
2"
0
...
0
37775 (17)
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 72 / 82
3. Parametric spectral estimation Signal modeling for rational spectra
AR processes — Yule-Walker equations (cont’d)
Using a matrix block formulation, the Yule-Walker equations can compactly be expressed as:26664
r(0) r(1) · · · r(p)
r(1) r(0) · · · r(p + 1)
...
...
. . .
...
r(p) r(p 1) · · · r(0)
37775
26664
1
a1
...
ap
37775 =
26664
2"
0
...
0
37775

r(0) rHp
rp Rp

1
✓p

=

2"
0p

where
rp
def
= [r(1), r(2), . . . , r(p)]T Rp
def
=
264 r(0) r(1) · · · r(p + 1)... ... . . . ...
r(p 1) r(p 2) · · · r(0)
375
✓p
def
= [a1, a2, . . . , ap ]
T 0p
def
= [0, 0, . . . , 0| {z }
p
]T
and (·)H denotes the Hermitian (conjugate-transpose) operator.
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 73 / 82
3. Parametric spectral estimation Signal modeling for rational spectra
AR processes — Yule-Walker equations (cont’d)
Hence, the Yule-Walker equations can be divided into two parts:
r(0) rHp
rp Rp

1
✓p

=

r(0) + rHp ✓p
rp + Rp✓p

=

2"
0p

From the bottom part, we can obtain the AR coecients ✓p by solving the linear system:
Rp✓p = rp (18)
Then, from the top part, we can compute the innovation variance 2" as:
2" = r(0) + r
H
p ✓p (19)
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 74 / 82
3. Parametric spectral estimation Signal modeling for rational spectra
Yule-Walker method for AR models
Yule-Walker (or autocorrelation) method
1 Select the model order p.
2 Compute suitable ACS estimates {rˆ(k)}pk=0 from the available samples {x(n)}N1n=0 .
3 Build autocorrelation matrix Rp and autocorrelation vector rp .
4 Solve Yule-Walker linear system (18) to estimate the AR model coecients ✓ˆp .
5 Estimate the innovation variance (or modeling error) ˆ2" through eqn. (19).
6 Compute the AR PSD estimate (16).
Remarks
Choosing the model order is a dicult problem ! seen later in the chapter
Biased ACS estimates are preferred: Rp positive definite ! lower variance in PSD estimate
Inverting matrix Rp ! O(p3) products-divisions (Gaussian elimination)...
... but Rp is a Hermitian Toeplitz matrix
! computationally ecient algorithm with O(p2) prod.-div.: Levinson-Durbin recursion
ACS estimates rˆ(k): O(Np) operations
! cost reduction of Levinson-Durbin algorithm may be negligible if N p.
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 75 / 82
3. Parametric spectral estimation Signal modeling for rational spectra
Link with linear prediction
Linear prediction: estimate or predict x(n) from a linear combination of previous samples
xˆ(n) =
pX
k=1
akx(n k)
Criterion: minimize the mean square error (MSE)
JMSE = E{|e(n)|2} (20)
e(n)
def
= xˆ(n) x(n): linear prediction error
{ak}pk=1: linear prediction coecients
A(z) = 1 +
Pp
k=1 akz
k : linear prediction filter
Linear prediction and AR modeling
The linear prediction coecients minimizing the MSE criterion (20) are given by the solution of
Yule-Walker equations (18) for AR modeling.
The linear prediction error is the innovation process: e(n) = "(n).
The minimum MSE is given by the innovation variance (19): min JMSE = 2".
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 76 / 82
3. Parametric spectral estimation Signal modeling for rational spectra
Moving average (MA) processes
If p = 0 ! ARMA(0, q) ! moving average (MA) process of order q: MA(q)
x(n) =
qX
k=0
bk"(n k) ) S(!) = 2"

qX
k=0
bke
|!k

2
Also known as all-zero model.
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 77 / 82
3. Parametric spectral estimation Signal modeling for rational spectra
MA processes — examples
0 50 100 150 200 250
−5
0
5
0 50 100 150 200 250
−10
−5
0
5
0 50 100 150 200 250
0
Figure: (Top) White noise input.
(Middle) MA process with b = [1, 1, 1, 1]/4.
(Bottom) MA process with b = [1, 1, . . . , 1]| {z }
100 entries
/100.
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 78 / 82
3. Parametric spectral estimation Signal modeling for rational spectra
MA processes — Yulker-Walker equations
For MA(q) processes, eqn. (12) reduces to H(z) =
Pq
k=0 bkz
k . Then:
h(k)
def
= Z1{H(z)} =
8<: bk 0  k  q0 otherwise.
Hence, the ACS (14) becomes
r(k) =
8>><>>:
2"
Pqk
m=0 bm+kb
⇤
m 0  k  q
r⇤(k) k < 0
0 |k| > q.
Because r(k) = 0 for |k| > q, a natural PSD estimate is:
Sˆ(!) =
qX
k=q
rˆ(k)e|!k
where rˆ(k) is a suitable ACS estimate.
Equivalent to Blackman-Tuckey with rectangular window of length 2q + 1 [Chap. 2]
Biased PSD estimate if process X (n) is not actually governed by an MA(q) model.
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 79 / 82
3. Parametric spectral estimation Model order selection
Model order selection
Important step in rational spectrum modeling: estimate model order (p, q) from observed data.
Selecting model order p in AR modeling
if too small ! smoothed spectrum, poor resolution
if too large ! spurious peaks: spectral line splitting.
General approach
increase model order until modeling error is minimized
problem: error is a monotonically nonincreasing function of p
idea: incorporate a penalty term that increases with model order p
! select p minimizing the criterion:
C(p) = N log 2"(p) + f (N)p
N : data record length
2"(p) : modeling error for model order p
f (N) : constant that may depend on N.
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 80 / 82
3. Parametric spectral estimation Model order selection
Model order selection (cont’d)
Akaike information criterion (AIC)
AIC(p) = N log 2"(p) + 2p
estimated p typically too small for non-AR processes
it tends to overestimate p as N increases
Minimum description length (MDL)
MDL(p) = N log 2"(p) + (logN)p
consistent model order estimator: pˆMDL !
N!1
p
Akaike’s final prediction error (FPE)
FPE(p) = 2"(p)
N + p + 1
N p 1
Remarks
no criterion works particularly well for short data sequences
generally they should just be used as ‘indicators’ of the model order
prediction error 2"(p) depends on modeling technique.
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 81 / 82
3. Parametric spectral estimation
Summary
Exploit prior knowledge about the random process to be analyzed
General approach
1 Select an appropriate model for process under analysis
2 Estimate the model parameters from the available data
3 Estimate the PSD by incorporating the estimated parameters into the parametric form for the PSD
PSD with rational function structure — 3 models:
I autoregressive moving average (ARMA): spectrum with poles and zeros
I autoregressive (AR): all-pole model
I moving average (MA): all-zero model
Model parameters can be estimated using Yule-Walker equations
I estimate ACS ! build autocorrelation matrix ! solve a linear system of equations
Several criteria for model order selection
I AIC, MDL, FPE, ...
Performance of parametric approach depends on fitness of model to the process being
analyzed.
V. Zarzoso PNS - ELEC4 - Spectral Analysis 2022–2023 82 / 82