Prony's Method is a quasi-least-squares time-domain
IIR filter design method.
First, assume
Hz
H
z
is an "all-pole" system:
Hz=
b
0
1+∑k=1M
a
k
z-k
H
z
b
0
1
k
1
M
a
k
z
k
(1)
and
hn=-∑k=1M
a
k
hn−k+
b
0
δn
h
n
k
1
M
a
k
h
n
k
b
0
δ
n
where
hn=0
h
n
0
,
n<0
n
0
for a causal system.
For
h=0
h
0
,
h0=
b
0
h
0
b
0
.
Let's attempt to fit a desired impulse response (let it be
causal, although one can extend this
technique when it isn't)
h
d
n
h
d
n
.
A true least-squares solution would attempt to
minimize
ε2=∑n=0∞|
h
d
n−hn|2
ε
2
n
0
h
d
n
h
n
2
where
Hz
H
z
takes the form in Equation 1. This is a
difficult non-linear
optimization problem which is known to be plagued by local
minima in the error surface. So instead of solving this
difficult non-linear problem, we solve the deterministic linear
prediction problem, which is related to, but
not the same as, the true least-squares optimization.
The deterministic linear prediction problem is a
linear least-squares optimization, which is
easy to solve, but it minimizes the
prediction error, not the
|desired−actual|2
desired
actual
2
response error.
Notice that for
n>0
n
0
, with the all-pole filter
hn=-∑k=1M
a
k
hn−k
hn
k
1
M
a
k
h
n
k
(2)
the right hand side of
this
equation is a
linear predictor of
hn
h
n
in terms of the
M
M previous samples of
hn
h
n
.
For the desired reponse
h
d
n
h
d
n
, one can choose the recursive filter coefficients
a
k
a
k
to minimize the squared prediction error
ε
p
2=∑n=1∞|
h
d
n+∑k=1M
a
k
h
d
n−k|2
ε
p
2
n
1
h
d
n
k
1
M
a
k
h
d
n
k
2
where, in practice, the ∞
is replaced by an N
N.
In matrix form, that's
h
d
00...0
h
d
1
h
d
0...0⋮⋮⋱⋮
h
d
N−1
h
d
N−2...
h
d
N−M
a
1
a
2
⋮
a
M
≈-
h
d
1
h
d
2⋮
h
d
N
h
d
0
0
...
0
h
d
1
h
d
0
...
0
⋮
⋮
⋱
⋮
h
d
N
1
h
d
N
2
...
h
d
N
M
a
1
a
2
⋮
a
M
h
d
1
h
d
2
⋮
h
d
N
or
H
d
a≈-hd
H
d
a
h
d
The optimal solution is
alp=-
H
d
H
H
d
-1
H
d
Hhd
a
lp
H
d
H
d
-1
H
d
h
d
Now suppose
Hz
H
z
is an
M
th
M
th
-order IIR (ARMA) system,
Hz=∑k=0M
b
k
z-k1+∑k=1M
a
k
z-k
H
z
k
0
M
b
k
z
k
1
k
1
M
a
k
z
k
or
hn=-∑k=1M
a
k
hn−k+∑k=0M
b
k
δn−k=-∑k=1M
a
k
hn−k+
b
n
if0≤n≤M-∑k=1M
a
k
hn−kifn>M
h
n
k
1
M
a
k
h
n
k
k
0
M
b
k
δ
n
k
k
1
M
a
k
h
n
k
b
n
0
n
M
k
1
M
a
k
h
n
k
n
M
(3)
For
n>M
n
M
, this is just like the all-pole case, so we can solve
for the best predictor coefficients as before:
h
d
M
h
d
M−1...
h
d
1
h
d
M+1
h
d
M...
h
d
2⋮⋮⋱⋮
h
d
N−1
h
d
N−2...
h
d
N−M
a
1
a
2
⋮
a
M
≈
h
d
M+1
h
d
M+2⋮
h
d
N
h
d
M
h
d
M
1
...
h
d
1
h
d
M
1
h
d
M
...
h
d
2
⋮
⋮
⋱
⋮
h
d
N
1
h
d
N
2
...
h
d
N
M
a
1
a
2
⋮
a
M
h
d
M
1
h
d
M
2
⋮
h
d
N
or
Ĥda≈ĥd
H
d
a
h
d
and
aopt=ĤdH
H
d
-1
H
d
Hĥd
a
opt
H
d
H
d
-1
H
d
h
d
Having determined the
a
a's, we can use them in
Equation 3 to obtain
the
b
n
b
n
's:
b
n
=∑k=1M
a
k
h
d
n−k
b
n
k
1
M
a
k
h
d
n
k
where
h
d
n−k=0
h
d
n
k
0
for
n−k<0
n
k
0
.
For
N=2M
N
2
M
,
Ĥd
H
d
is square, and we can solve exactly
for the
a
k
a
k
's with no error. The
b
k
b
k
's are also chosen such that there is no error in the
first
M+1
M
1
samples of
hn
h
n
. Thus for
N=2M
N
2
M
, the first
2M+1
2
M
1
points of
hn
h
n
exactly equal
h
d
n
h
d
n
. This is called Prony's Method. Baron de
Prony invented this in 1795.
For
N>2M
N
2
M
,
h
d
n=hn
h
d
n
h
n
for
0≤n≤M
0
n
M
, the prediction error is minimized for
M+1<n≤N
M
1
n
N
, and whatever for
n≥N+1
n
N
1
. This is called the Extended Prony
Method.
One might prefer a method which tries to minimize an
overall error with the numerator coefficients, rather than just
using them to exactly fit
h
d
0
h
d
0
to
h
d
M
h
d
M
.
- Assume an all-pole model and fit
h
d
n
h
d
n
by minimizing the prediction error
1≤n≤N
1
n
N
.
- Compute
vn
v
n
, the impulse response of this all-pole
filter.
- Design an all-zero (MA, FIR) filter which
fits
vn*
h
z
n≈
h
d
n
v
n
h
z
n
h
d
n
optimally in a least-squares sense (Figure 1).
The final IIR filter is the cascade of the all-pole and
all-zero filter.
This is is solved by
min
b
k
{∑n=0N|
h
d
n−∑k=0M
b
k
vn−k|2}
b
k
n
0
N
h
d
n
k
0
M
b
k
v
n
k
2
or in matrix form
v000...0v1v00...0v2v1v0...0⋮⋮⋮⋱⋮vNvN−1vN−2...vN−M
b
0
b
1
b
2
⋮
b
M
≈
h
d
0
h
d
1
h
d
2⋮
h
d
N
v
0
0
0
...
0
v
1
v
0
0
...
0
v
2
v
1
v
0
...
0
⋮
⋮
⋮
⋱
⋮
v
N
v
N
1
v
N
2
...
v
N
M
b
0
b
1
b
2
⋮
b
M
h
d
0
h
d
1
h
d
2
⋮
h
d
N
Which has solution:
bopt=VHV-1VHh
b
opt
V
V
-1
V
h
Notice that none of these methods solve the true
least-squares problem:
mina,b{∑n=0∞|
h
d
n−hn|2}
a
b
n
0
h
d
n
h
n
2
which is a difficult non-linear optimization problem. The true
least-squares problem can be written as:
minα,β{∑n=0∞|
h
d
n−∑i=1M
α
i
ⅇ
β
i
n|2}
α
β
n
0
h
d
n
i
1
M
α
i
β
i
n
2
since the impulse response of an IIR filter is a sum of
exponentials, and non-linear optimization is then used to
solve for the
α
i
α
i
and
β
i
β
i
.
Digital-to-Digital Frequency Transformations
Overview of IIR Filter Design
The DFT, FFT, and Practical Spectral Analysis
