System Identification

On Jul 1, 6:02=A0am, Rune Allnor <all...@tele.ntnu.no> wrote:

> Therrien suggestes that the term 'Prony method' applies to
> ARMA estimation where the A() and B() coefficents are
> estimated in separable steps.

At the risk of embarrassing myself, since it's been literally 25 years
since I worked with this stuff, I recall Prony's Method to model the
impulse response of an ARMA system as the sum of damped sinusoids.
There was also Pade-Prony, which did the same thing but in a z-
transform formulation.  Either way, the assumption of damped sinusoids
would indicate that this method might be suitable for modeling of
transients.

It's all just a series of foggy memories.

Greg

On Jul 1, 6:02&#4294967295;am, Rune Allnor <all...@tele.ntnu.no> wrote:

> Therrien suggestes that the term 'Prony method' applies to
> ARMA estimation where the A() and B() coefficents are
> estimated in separable steps.

At the risk of embarrassing myself, since it's been literally 25 years
since I worked with this stuff, I recall Prony's Method to model the
impulse response of an ARMA system as the sum of damped sinusoids.
There was also Pade-Prony, which did the same thing but in a z-
transform formulation.  Either way, the assumption of damped sinusoids
would indicate that this method might be suitable for modeling of
transients.

It's all just a series of foggy memories.

Greg

On 1 Jul, 13:15, Greg Berchin <gberc...@sentientscience.com> wrote:
> On Jul 1, 6:02&#4294967295;am, Rune Allnor <all...@tele.ntnu.no> wrote:
>
> > Therrien suggestes that the term 'Prony method' applies to
> > ARMA estimation where the A() and B() coefficents are
> > estimated in separable steps.
>
> At the risk of embarrassing myself, since it's been literally 25 years
> since I worked with this stuff, I recall Prony's Method to model the
> impulse response of an ARMA system as the sum of damped sinusoids.

That's correct. The original works by Prony in the 1790s was
on damped sinusoids, oscillations in bubbles. He formulated
his solution in terms of roots of polynomial 20 years before
Abel proved that roots of 5th order polynomoials could not be
found by means of arithmetics, 50 years before Babbage's
rudimentary computer (and 150 years before the first working
computer), 15 years before Gauss came up with the least squares
formulation of parameter estimation. Ah, and Prony (or Gaspard
Riche, as his birthname was) was a baron in the midst of the
French revolution.

There are lots of reasons why Prony's work could have been
forgotten over the past two centuries. His divide-and-conquer
trick is the reason why his name is still remembered.

> There was also Pade-Prony, which did the same thing but in a z-
> transform formulation. &#4294967295;Either way, the assumption of damped sinusoids
> would indicate that this method might be suitable for modeling of
> transients.

Sure.

Rune

On 1 Jul., 13:48, Rune Allnor <all...@tele.ntnu.no> wrote:
> On 1 Jul, 13:15, Greg Berchin <gberc...@sentientscience.com> wrote:
>
> > On Jul 1, 6:02&#4294967295;am, Rune Allnor <all...@tele.ntnu.no> wrote:
>
> > > Therrien suggestes that the term 'Prony method' applies to
> > > ARMA estimation where the A() and B() coefficents are
> > > estimated in separable steps.
>
> > At the risk of embarrassing myself, since it's been literally 25 years
> > since I worked with this stuff, I recall Prony's Method to model the
> > impulse response of an ARMA system as the sum of damped sinusoids.
>
> That's correct. The original works by Prony in the 1790s was
> on damped sinusoids, oscillations in bubbles. He formulated
> his solution in terms of roots of polynomial 20 years before
> Abel proved that roots of 5th order polynomoials could not be
> found by means of arithmetics, 50 years before Babbage's
> rudimentary computer (and 150 years before the first working
> computer), 15 years before Gauss came up with the least squares
> formulation of parameter estimation. Ah, and Prony (or Gaspard
> Riche, as his birthname was) was a baron in the midst of the
> French revolution.

<nit>
Prony published his algorithm 1795 [1], which is exactly the year that
Gauss used least-squares to estimate the trajectory of Ceres. Gauss
himself never published this, but only commented in a dismissive
manner on Legendre's publication of the least-squares algorithm in
1805. Babbage proposed his Difference Engine in 1822 (Pascal's
calculator was built in 1642 :-).
</nit>

Prony's method did not propose to use least-squares solutions for the
linear problems in his three step programme to find amplitudes,
damping factors, frequencies and phases of sums of sinusoids. It is a
really neat algorithm that works perfectly under ideal (=noise free)
conditions. Many people believe that Prony's method is in fact only
the first step of his algorithm. The first step is the AR modeling,
which essentially computes the linear prediction coefficients from the
impulse response (or any other finite excitation). Prony's approach
can easily be generalized to include least-squares solutions for the
linear parts (AR coefficient estimation and phase/amplitude
estimation). The major drawback is that even the extend Prony's method
is super-ultra-sensitive to additive measurement noise (this comes
from the polynomial factoring part of the algorithm). Sometimes, this
is compensated with over-modelling (hint for Vlad) - however, use this
with caution (what happens to the predictor polynomial?). As Rune
said, Lawrence has a nice chapter devoted to this.

If the first step in Prony's algorithm is formulated as a Kalman
filter, both the state itself and the state transmission matrix (built
from the AR coefficients) are unkown, leading to a non-linear Kalman
representation. The extended Kalman filter in is known not to converge
well in that case, the "best" approach for estimating AR signals under
additive measurement noise is still open. I used a dual-linear
approach in a recent project, combining a Kalman filter with a robust
(against additive measurement noise) AR estimation algorithm in an
iterative fashion. Messy, but works - especiall tuned for estimating
the parameters from responses to impulsive excitations, where the
decay of the response means that one does not have an arbitrary amount
of data available for the parameter estimation.

I do believe our very own Julius Kusuma's PhD thesis uses a modified
version of Prony's to find parameters of sparse signals (Dirac deltas
with Poisson distribute inter-arrival times). Perhaps he has more to
add here.

Regards,
Andor

[1] R. Prony, &#4294967295;Essai exp&#4294967295;rimental et analytique,&#4294967295; Ann. &#4294967295;cole
Polytechnique, vol. 1,
no. 2, p. 24, 1795.

On 1 Jul, 14:24, Andor <andor.bari...@gmail.com> wrote:

> Prony's method did not propose to use least-squares solutions for the
> linear problems in his three step programme to find amplitudes,
> damping factors, frequencies and phases of sums of sinusoids. It is a
> really neat algorithm that works perfectly under ideal (=noise free)
> conditions.

There are lots of variations. The basic Prony's method
formulation says that the data set contains N samples and
the sum-of-sines contains N/2 (complex) sinusoidals.

I don't remember the details off the top of my head, but
there is also an 'extended Prony's method' and a 'modified
Prony's method.' One takes into account that there are
more samples than sinusoidal (e.g. dealing with over-
determined systems) and the other accounts for the fact
that the exact number of sinusoidals usually is not known,
and thus introduces a 'model mismatch' term.

And then we can start discussing exactly what terms are
included in the model; frequency, amplitude, phase,
damping...

> Many people believe that Prony's method is in fact only
> the first step of his algorithm. The first step is the AR modeling,
> which essentially computes the linear prediction coefficients from the
> impulse response (or any other finite excitation).

Careful. There is a trick here if you deal with the standard
sum-of-sines models. In order to make the maths work, one
needs to use the modified covariance equations, not the usual
covariance equations (details in the 1982 Proc IEEE paper by
Kay and Marple). This has a number of consequences:

- The covariance matrix is no longer Toeplitz
- The covariance matrix is no longer unconditionally stable
- The method can be used for spectrum line estimation

> Prony's approach
> can easily be generalized to include least-squares solutions for the
> linear parts (AR coefficient estimation and phase/amplitude
> estimation). The major drawback is that even the extend Prony's method
> is super-ultra-sensitive to additive measurement noise

Correct.

> (this comes
> from the polynomial factoring part of the algorithm).

Wrong. There was a different approach suggested in 1982 by
Tufts and Kumaresan. They used the same signal model as the
'usual' variations of Prony's method, but they chose a
different strategy for solving for the *coefficients*, not
roots, of the predictor polynomials. I tested the T&K
method against one of Marple's routines in my PhD thesis.
The T&K methods stil worked at SNRs 20dB lower than
where the Marple method had given in. The only difference
is how they compute the coefficients of the predictor
polynomial from the data. Everything else is the same.

Interestingly, with Vladimir's data the situation was
the oposite: The results from T&K made no sense whereas
the Marple algorithm apparently hit quite close to home.

> Sometimes, this
> is compensated with over-modelling (hint for Vlad) - however, use this
> with caution (what happens to the predictor polynomial?).

That's where the art lies; choosing the correct model orders.
That single step makes or brakes the analysis. With Vladimir's
data there are no problems, but in some applications the
number of data samples impose very limiting bounds on what
model orders can be chosen.

Rune

Vladimir Vassilevsky <antispam_bogus@hotmail.com> wrote:
> 
> 
> Rune Allnor wrote:
>> On 30 Jun, 21:09, Vladimir Vassilevsky <antispam_bo...@hotmail.com>
>> wrote:
>> 
>>>Rune Allnor wrote:
>> 
>> 
>>>Here is the typical one:
>>>
>>>http://www.abvolt.com/misc/data.cpp

Full segment analysis can give two problems:
- tail - almost only noise visible
- head (about 10-20 samples) - nonlinearity or undamped higher frequency
  componets

>> 
>> 
>> Greg already suggested Prony's method, and I agree with that.

Prony at low SNR gives bad results.

My best fit  for Vladimir's data  (full segment - 1000 points) suggests that
SNR 2-pole-model/residuum is only ~17dB.

You can improve estimation using only samples 20:200 - SNR 40dB

> The correct answer for this example is two complex conjugate poles at 
> F=Fc/50 with Q = 1.7 . However I am getting wildly different results 
> depending on the sample rate.

My favourite is Sarkar's matrix pencil method (SVD based), 
an implementation example:
http://groups.google.pl/group/comp.soft-sys.matlab/msg/c71f1639d5253b1c

My results with this function are below.

Some details:
- analysed segment data(20:n:200) - segment choosen by visual inspection
  of a plot log(abs(hilbert(signal)))
- n-downsampling factor in range 1-25, for n=25 only seven samples analysed ;)
- downsampling done without filtering 
- estimated model with 4-poles (with only 2 you can get spurious results 
  in this case)

n= 1 Q=1.647 f=0.01929
n= 2 Q=1.645 f=0.01929
n= 3 Q=1.645 f=0.01929
n= 4 Q=1.646 f=0.01930
n= 5 Q=1.644 f=0.01929
n= 6 Q=1.644 f=0.01929
n= 7 Q=1.644 f=0.01931
n= 8 Q=1.647 f=0.01930
n= 9 Q=1.649 f=0.01934
n=10 Q=1.638 f=0.01926
n=11 Q=1.642 f=0.01928
n=12 Q=1.650 f=0.01930
n=13 Q=1.660 f=0.01930
n=14 Q=1.640 f=0.01928
n=15 Q=1.637 f=0.01925
n=16 Q=1.667 f=0.01942
n=17 Q=1.643 f=0.01929
n=18 Q=1.673 f=0.01938
n=19 Q=1.673 f=0.01941
n=20 Q=1.632 f=0.01928
n=21 Q=1.701 f=0.01934
n=22 Q=1.763 f=0.01927
n=23 Q=1.477 f=0.01916
n=24 Q=1.671 f=0.01929
n=25 Q=1.766 f=0.02000 

Mirek Kwasniak

Rune Allnor <allnor@tele.ntnu.no> wrote in news:8f2a3224-efa0-461e-904b-
ded21d38864d@a1g2000hsb.googlegroups.com:

> As far as I can tell, 'Prony's method' applies to dividing
> the multiple-parameter estimation problem into separate
> steps where each parameter is estimated independently
> from the others (from a numerical sense).
> 

Isn't this also known as "curve peeling"?

-- 
Scott
Reverse name to reply

On 1 Jul, 12:02, Rune Allnor <all...@tele.ntnu.no> wrote:
> On 1 Jul, 05:14, Vladimir Vassilevsky <antispam_bo...@hotmail.com>
> wrote:

> > > I have a version lying around which was based on stationary
> > > data (whereas your data is a transient.)
>
> > What is fundamentally different between the cases of impulse response
> > and the stationary data?
>
> The *theoretical* discussions in Therrien doesn't make the
> distinction. Practical implementations are very different
> in that they are based on different data models. The one
> I used yesterday was based on a sum-of-sinusoidals model
> on the form
>
> &#4294967295; &#4294967295; &#4294967295; &#4294967295; P
> x[n] = sum a_p exp (j w_p n + theta_p)
> &#4294967295; &#4294967295; &#4294967295; &#4294967295;p=1
>
> The model does not include noise nor exponential damping
> terms. When I implemented that processor I had a book
> available where a number of variations was available,
> including one which included the exponential damping terms.

This thread got me thinking. First of all, the Marple book,
published a mere 20 years ago  has been out of reach for
over 15 years already. The Tufts&Kumaresan methods which
beat Marple hands down didn't work.

So I went back to look over the T&K methods, and sure
enough, the methods I used for my PhD thesis were configured
for the forward-backward prediction problem. Once I adjusted
them for forward-only prediction, Prony's method based on
the T&K implementation worked like a charm.

Vladimir, if you are intersted in discussing this in more
detail in private, drop a hint here. If you post from a valid
email address I'll send you a valid email address of mine
via my google posting account; the email address in my
posts has been obsolete for years.

Rune

pisz_na.mirek@dionizos.zind.ikem.pwr.wroc.pl wrote:

> Vladimir Vassilevsky <antispam_bogus@hotmail.com> wrote:
> Here is the typical one:
>
>http://www.abvolt.com/misc/data.cpp
> 
> 
> Full segment analysis can give two problems:
> - tail - almost only noise visible
> - head (about 10-20 samples) - nonlinearity or undamped higher frequency
>   componets
> 
> Prony at low SNR gives bad results.
> 
> My best fit  for Vladimir's data  (full segment - 1000 points) suggests that
> SNR 2-pole-model/residuum is only ~17dB.
> 
> You can improve estimation using only samples 20:200 - SNR 40dB

I checked; this is not the SNR issue. The approach was incorrect: I was 
trying to build Wiener filter by AR method without having the meaningful 
data sufficiently long. So it worked very well for the longer impulse 
responses (Q ~ 10) and goofed for the shorter ones.

> My favourite is Sarkar's matrix pencil method (SVD based), 
> an implementation example:
> http://groups.google.pl/group/comp.soft-sys.matlab/msg/c71f1639d5253b1c

What book would you suggest on this method?

> Some details:
> - analysed segment data(20:n:200) - segment choosen by visual inspection
>   of a plot log(abs(hilbert(signal)))

You have my respect for log|Hilbert| envelope. It is nice to see the 
professional approach.

> - n-downsampling factor in range 1-25, for n=25 only seven samples analysed ;)
> - downsampling done without filtering 
> - estimated model with 4-poles (with only 2 you can get spurious results 
>   in this case)
> 
> n= 1 Q=1.647 f=0.01929

[...]

Neat. Right on target.

> Mirek Kwasniak

Thank you, Dr. Kwasniak


Vladimir Vassilevsky
DSP and Mixed Signal Design Consultant
http://www.abvolt.com

Rune Allnor wrote:

> This thread got me thinking. First of all, the Marple book,
> published a mere 20 years ago  has been out of reach for
> over 15 years already. The Tufts&Kumaresan methods which
> beat Marple hands down didn't work.

Just returned from the local univ. library. Marple is not there.

> So I went back to look over the T&K methods, and sure
> enough, the methods I used for my PhD thesis were configured
> for the forward-backward prediction problem. Once I adjusted
> them for forward-only prediction, Prony's method based on
> the T&K implementation worked like a charm.

The problem was because the impulse response was too short for using the 
AR method as a brute force.

> Vladimir, if you are intersted in discussing this in more
> detail in private, drop a hint here.

I would certainly be interested.

> If you post from a valid
> email address I'll send you a valid email address of mine
> via my google posting account; the email address in my
> posts has been obsolete for years.

My valid email and phone contact are at the web site in the signature; 
please send from some account other then gmail: the gmail is likely to 
be filtered out by spamkillers.

> 
> Rune

Vladimir Vassilevsky
DSP and Mixed Signal Design Consultant
http://www.abvolt.com