FPGA Central

Hello All,

There is a physical system which is described well by the 2nd order pole
only model. The poles can be complex (Q up to 10), or real, including
the degenerated cases with one pole or no poles at all. The impulse
response can be measured.

The problem is to find the system parameters from the impulse response
as accurate as possible. There are methods based on the Bode plots,
min/max peaks, slopes, zero crossing position estimations. However they
look pretty much ad-hoc, require the initial decision about the system
type, and prone to errors.

I made the AR estimator: compute the autocorrelation function, then
Levinson-Durbin to get the model coefficients, then solve for the poles.
However it appears that the result is dependent on the sample rate,
especially in the case of low Q poles. This is not a numeric issue. The
result could be completely wrong with the higher oversampling, and this
looks to be the problem of the approach.

What would be the right way to estimate the parameters?


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

On 30 Jun, 19:50, Vladimir Vassilevsky <antispam_bo...@hotmail.com>
wrote:
> Hello All,
>
> There is a physical system which is described well by the 2nd order pole
> * only model. The poles can be complex (Q up to 10), or real, including
> the degenerated cases with one pole or no poles at all. The impulse
> response can be measured.
>
> The problem is to find the system parameters from the impulse response
> as accurate as possible. There are methods based on the Bode plots,
> min/max peaks, slopes, zero crossing position estimations. However they
> look pretty much ad-hoc, require the initial decision about the system
> type, and prone to errors.
>
> I made the AR estimator: compute the autocorrelation function, then
> Levinson-Durbin to get the model coefficients, then solve for the poles.
> However it appears that the result is dependent on the sample rate,
> especially in the case of low Q poles. This is not a numeric issue. The
> result could be completely wrong with the higher oversampling, and this
> looks to be the problem of the approach.
>
> What would be the right way to estimate the parameters?

I would pursue the AR line, with focus on numerical methods.
Use an as long observation record as you dare (you don't
mention whether this is stationary) and get a good estimate
for the signal autocovariance matrix. The matrix ought to
be dimension at least 6 x 6, provided you are right that
there is at most one pair of poles.

You could try the Levinson recursion + order estimators
with succesively longer frames and see if better statistics
help.

My second approach would be to use the SVD to compute the
parameters for the AR model. Set up a data matrix (details
in the 1982 paper by Tufts and Kumaresan) compute the
SVD and look at the singular values to determine the order.

These sorts of inquiries *might* work, but may also be
too computationally expensive for your application.

Is there any chance that you post example data somewhere?
Preferably both easy and tricky ones. Ah, and a control,
where you have modeled the data (both easy and tricky)
and know the answers. If you could post these data
somewhere and not disclose what data set is which, it
would be very interesting to see what might be done with
them.

Rune

On Jun 30, 1:50*pm, Vladimir Vassilevsky <antispam_bo...@hotmail.com>
wrote:

> The problem is to find the system parameters from the impulse response
> as accurate as possible.

How noisy is your impulse response characterization? Prony's Method
does exactly this, but becomes inconsistent in the presence of noise.

Greg

Rune Allnor wrote:
> On 30 Jun, 19:50, Vladimir Vassilevsky <antispam_bo...@hotmail.com>
> wrote:
>
>>There is a physical system which is described well by the 2nd order pole
>> only model.
>>The problem is to find the system parameters from the impulse response
>>as accurate as possible.
> Is there any chance that you post example data somewhere?

> Preferably both easy and tricky ones. Ah, and a control,
> where you have modeled the data (both easy and tricky)
> and know the answers.

Here is the typical one:

http://www.abvolt.com/misc/data.cpp

> If you could post these data
> somewhere and not disclose what data set is which, it
> would be very interesting to see what might be done with
> them.

VLV

On Jun 30, 1:50 pm, Vladimir Vassilevsky <antispam_bo...@hotmail.com>
wrote:
> Hello All,
>
> There is a physical system which is described well by the 2nd order pole
> only model. The poles can be complex (Q up to 10), or real, including
> the degenerated cases with one pole or no poles at all. The impulse
> response can be measured.
>
> The problem is to find the system parameters from the impulse response
> as accurate as possible. There are methods based on the Bode plots,
> min/max peaks, slopes, zero crossing position estimations. However they
> look pretty much ad-hoc, require the initial decision about the system
> type, and prone to errors.
>
> I made the AR estimator: compute the autocorrelation function, then
> Levinson-Durbin to get the model coefficients, then solve for the poles.
> However it appears that the result is dependent on the sample rate,
> especially in the case of low Q poles. This is not a numeric issue. The
> result could be completely wrong with the higher oversampling, and this
> looks to be the problem of the approach.
>
> What would be the right way to estimate the parameters?
>
> Vladimir Vassilevsky
> DSP and Mixed Signal Design Consultanthttp://www.abvolt.com

Sample rate dependency implies you might have some MA in the system.
Additive measurement noise also makes AR systems ARMA.

On Jun 30, 1:50*pm, Vladimir Vassilevsky <antispam_bo...@hotmail.com>
wrote:
> Hello All,
>
> There is a physical system which is described well by the 2nd order pole
> * only model. The poles can be complex (Q up to 10), or real, including
> the degenerated cases with one pole or no poles at all. The impulse
> response can be measured.
>
> The problem is to find the system parameters from the impulse response
> as accurate as possible. There are methods based on the Bode plots,
> min/max peaks, slopes, zero crossing position estimations. However they
> look pretty much ad-hoc, require the initial decision about the system
> type, and prone to errors.
>
> I made the AR estimator: compute the autocorrelation function, then
> Levinson-Durbin to get the model coefficients, then solve for the poles.
> However it appears that the result is dependent on the sample rate,
> especially in the case of low Q poles. This is not a numeric issue. The
> result could be completely wrong with the higher oversampling, and this
> looks to be the problem of the approach.
>
> What would be the right way to estimate the parameters?
>
> Vladimir Vassilevsky
> DSP and Mixed Signal Design Consultanthttp://www.abvolt.com

Hello Vlad,

An LPC method should yield your two coefs to be r*cos(theta) and -r^2.
The sample rate will certainly affect the both the "r" (decay rate per
sample) and cos(theta) "angular step per sample) values. Are you
getting something different?

Clay

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

Greg already suggested Prony's method, and I agree with that.
I have a version lying around which was based on stationary
data (whereas your data is a transient.) My version computes
this 10'th order prediction operator:

gv = 1.0000
0.4854
-5.7894
5.1015
1.1272
-5.1854
4.3738
1.5878
-3.8126
0.7763
0.3344

If you plot the root locations of this polynomial
on top of the spectrum of your data,

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% d is the posted data,
% gv is the polynomial above

N=length(d);
thv=(0:N-1)/N*2*pi-pi;
clf
plot(thv,20*log10(abs(fftshift(fft(d)))),'b')
R=roots(gv);
[rth,r]=cart2pol(real(R),imag(R));
hold on
ax=axis;
axis([-pi pi ax(3:4)])
plot([1;1]*rth',ax(3:4)'*ones(1,length(rth)),'r');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

you can see that the roots line up pretty
well with the features in the data, if not
quite, at low frequencies.

I would go one step further and use a version of
Prony's method that work with transients. I don't
have an implementation for that right now, and
the one recipe I know of was out of print when
I hunted for it in 1994-95.

Anyway, it seems a Transient version of Prony's
method might be worth investigating.

Rune

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
>
>
> Greg already suggested Prony's method, and I agree with that.

I guess what I am trying to do is called Prony method, is it?

> 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?

> My version computes
> this 10'th order prediction operator:
>
> gv = 1.0000
> 0.4854
> -5.7894
> 5.1015
> 1.1272
> -5.1854
> 4.3738
> 1.5878
> -3.8126
> 0.7763
> 0.3344
>

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.


> I would go one step further and use a version of
> Prony's method that work with transients. I don't
> have an implementation for that right now, and
> the one recipe I know of was out of print when
> I hunted for it in 1994-95.
>
> Anyway, it seems a Transient version of Prony's
> method might be worth investigating.

I searched on Transient Prony Method and didn't get any direct hits. Can
you suggest a reading ?



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

On 1 Jul, 05:14, Vladimir Vassilevsky <antispam_bo...@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
>
> > Greg already suggested Prony's method, and I agree with that.
>
> I guess what I am trying to do is called Prony method, is it?

The term 'Prony's method' relates to a certain divide
and conquer strategy. While I used the Prony method
for half a decade during the work on my PhD thesis,
I never found a 100% authorative definition of the term.

After I posted yesterday I found some discussions in

Therrien: "Discrete Random Data and Statistical Signal
Processing", Prentice Hall 1992.

Therrien suggestes that the term 'Prony method' applies to
ARMA estimation where the A() and B() coefficents are
estimated in separable steps. When I used the Prony method
I saw a statement somewhere that the term applies to
situations where AR (not ARMA) models are handled by
estimating root locations and amplitudes in separate
steps.

The Levinson recursion for AR models estimate one predictor
polynomial where both amplitude and root locations are
embedded once and for all.

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).

> > 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

P
x[n] = sum a_p exp (j w_p n + theta_p)
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.

....
> > I would go one step further and use a version of
> > Prony's method that work with transients. I don't
> > have an implementation for that right now, and
> > the one recipe I know of was out of print when
> > I hunted for it in 1994-95.
>
> > Anyway, it seems a Transient version of Prony's
> > method might be worth investigating.
>
> I searched on Transient Prony Method and didn't get any direct hits. Can
> you suggest a reading ?

Anything by Lawrence Marple. The book I tried to hunt
down 15 years ago (and never found except in the
university library) was his 1987 book on spectrum
estimation. It contained a lot of those half-weird
variations. A 2nd hand copy was offered on amazon.com
for >$300 yesterday.

The Therrien book where I found the discussion on ARMA
models contains a section on matcing impulse responses.
The basic recipes are straight-forward, I havent tried
to implement them so I don't know how many devils lurk
in the details. However, I got that book myself in -94
but it was out of stock by -98 when I wanted some of my
students to buy it.

Anyway, I tracked down a paper by Therrien and a student
(IEEE Trans. Sig. Proc. 1996) where they match impulse
responses by starting out from a 'crude' model (probably
one from his book) and improved on it by applying some
sort of iterative adaption routine. I haven't had the
time to look at it in detail yet, but the improvment
demonstrated in the paper was impressive.

HTH,

Rune

On Jul 1, 6:02*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*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. *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*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, “Essai expérimental et analytique,” Ann. É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 <[email protected]> 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.s...1f1639d5253b1c

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 <[email protected]> wrote in news:8f2a3224-efa0-461e-904b-
[email protected]:

> 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
>
> * * * * P
> x[n] = sum a_p exp (j w_p n + theta_p)
> * * * *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

[email protected] wrote:

> Vladimir Vassilevsky <[email protected]> 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.s...1f1639d5253b1c

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

On 1 Jul, 23:45, Rune Allnor <all...@tele.ntnu.no> wrote:
> 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
>
> > * * * * P
> > x[n] = sum a_p exp (j w_p n + theta_p)
> > * * * *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.

Yep it did indeed. So why not post some of the results.
The matlab script starts out with Vladimir's data in the
variable d and the prediction coefficients from the T&K
variation of Prony's method in the vector av.

Rune
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%

av = [1.00000000000000;
-1.45804313154822;
0.29976575577196;
1.16543790050656;
-1.53927552893326;
0.16385197694451;
0.96374627878772;
-0.74226923864693;
-0.08460783603794;
0.26971433410238;
0.03199109349585];

r=roots(av).';
[rth,R]=cart2pol(real(r),imag(r));
Vv = zeros(length(d),length(r));
Vw = Vv;
for n=1:length(d)
Vv(n,=(n-1)*rth;
Vw(n,=R.^(n-1);
end
[Vr,Vi]=pol2cart(Vv,Vw);
V = Vr+sqrt(-1)*Vi;
A=pinv(V)*d;
dd=real(V*A);

clf
subplot(3,1,1)
plot(d,'b')
hold on
plot(dd,'r')

thv=(0:length(d)-1)/length(d);
subplot(3,1,2)
plot(thv,20*log10(abs(fft(d))),'b');
hold on
plot(thv,20*log10(abs(fft(dd))),'r');

subplot(3,1,3)
plot(real(r),imag(r),'x')
set(gca,'dataaspectratio',[1,1,1])
hold on
Thv = [0:360]/(2*pi);
Rv = ones(1,361);
[xx,yy]=pol2cart(Thv,Rv);
plot(xx,yy,'k')

On 2 Jul., 00:21, Vladimir Vassilevsky <antispam_bo...@hotmail.com>
wrote:
> pisz_na.mi...@dionizos.zind.ikem.pwr.wroc.pl wrote:
> > Vladimir Vassilevsky <antispam_bo...@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.s...1f1639d5253b1c
>
> 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 Consultan

Matrix pencil and Prony's variants have been suggested. You might also
consider the freely availble Harminv program to tackle this problem.
It is based on yet another approach ("filter diagonalization"), and
available here:

http://ab-initio.mit.edu/wiki/index.php/Harminv

It would be interesting to shoot out all those variants under
different conditions (SNRs, sampling rates, etc.).

Regards,
Andor

Rune Allnor <[email protected]> wrote in news:5d66b163-db36-43c5-926e-
[email protected]:

>
> 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.

There was a next edition of Marple threatened for a while, but I don't
think it ever materialized. Too bad, good book.



--
Scott
Reverse name to reply

On 2 Jul, 09:27, Andor <andor.bari...@gmail.com> wrote:

> Matrix pencil and Prony's variants have been suggested. You might also
> consider the freely availble Harminv program to tackle this problem.
> It is based on yet another approach ("filter diagonalization"), and
> available here:
>
> http://ab-initio.mit.edu/wiki/index.php/Harminv

You forgot your Kalman filter. Or do you count it as a
Prony variant? BTW, do you have anything written down
about that approach?

> It would be interesting to shoot out all those variants under
> different conditions (SNRs, sampling rates, etc.).

That's a very good idea. The best way to do this is in
a blind test scenario, where the analysts know nothing
about the data at the outset. In those sorts of tests
I prefer only to know the data by reference (dataA, dataB,
etc.) and file format.

So if somebody wants to organize a test like that, and
prepare and publish a few data sets, I'm in. Ah, yes,
a time should be set up front when the true parameters
of the data are disclosed, preferably allowing pwople
at least a week, maybe two, to test and tinker before
the true parameters are disclosed.

Rune

On 2 Jul., 14:42, Rune Allnor <all...@tele.ntnu.no> wrote:
> On 2 Jul, 09:27, Andor <andor.bari...@gmail.com> wrote:
>
> > Matrix pencil and Prony's variants have been suggested. You might also
> > consider the freely availble Harminv program to tackle this problem.
> > It is based on yet another approach ("filter diagonalization"), and
> > available here:
>
> >http://ab-initio.mit.edu/wiki/index.php/Harminv
>
> You forgot your Kalman filter. Or do you count it as a
> Prony variant? BTW, do you have anything written down
> about that approach?

Well, the Kalman filter is only part of calculating the AR
coefficients, the rest is the standard Prony method. I would classify
all algorithms that somehow compute the AR coefficients and proceed
from there as Prony variants. A colleague of mine tried the Unscented
Kalman filter as an alternative for the non-converging Extended Kalman
Filter. His results are still pending.

> > It would be interesting to shoot out all those variants under
> > different conditions (SNRs, sampling rates, etc.).
>
> That's a very good idea. The best way to do this is in
> a blind test scenario, where the analysts know nothing
> about the data at the outset. In those sorts of tests
> I prefer only to know the data by reference (dataA, dataB,
> etc.) and file format.
>
> So if somebody wants to organize a test like that, and
> prepare and publish a few data sets, I'm in. Ah, yes,
> a time should be set up front when the true parameters
> of the data are disclosed, preferably allowing pwople
> at least a week, maybe two, to test and tinker before
> the true parameters are disclosed.

That sounds doable. I would suggest a two-stage experiment: during the
first stage nothing is known about the signal. For the second-stage
the order of the model is revealed. This also allows to compare
different order estimation methods. The consistent order estimation
method for Prony's would be the (numerical) rank of the covariance
matrix.

We also need test signals. I can certainly supply a few from the
hundreds of simulations we did. However, the more test signals, the
merirer the competition :-). We also need a target model type, an
error functional and some kind of evaluation procedure to rank the
contestants("send in your program and we'll compute the error") ...

Regards,
Andor

While trying to find the parameters of the pole only system of the 2nd
order from the impulse response, I am running into the following problem:

If the coefficients of the system are calculated using the AR method,
the solution appears to be wrong.

But if the coefficients of the system are found by search for the best
matching impulse response of the model to the data, the solution is correct.

The AR minimizes the sum of errors at every step, whereas the impulse
response matching optimizes the whole thing. They are not equvalent.

Didn't figured out what to do yet.

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

On 8 Jul, 16:01, Vladimir Vassilevsky <antispam_bo...@hotmail.com>
wrote:
> While trying to find the parameters of the pole only system of the 2nd
> order from the impulse response, I am running into the following problem:
>
> If the coefficients of the system are calculated using the AR method,
> the solution appears to be wrong.

The spectrum of the data you posted last week doesn't comply to the
2nd order model. There are at least one set of spectrum lines (f=0.2)
and a 'hill' (f=0.4) in the spectrum. That's three poles right there.
Since the data are real, you need 3 more poles for the complex
conjugates. Try an AR model of order 6 or higher, I'd say 8-10 or so.

> But if the coefficients of the system are found by search for the best
> matching impulse response of the model to the data, the solution is correct.

Sure. One pair of poles completely domminates the impulse response.

> The AR minimizes the sum of errors at every step, whereas the impulse
> response matching optimizes the whole thing. They are not equvalent.
>
> Didn't figured out what to do yet.

You need to approach the data on their own terms. The fact that
the end target is a 2-pole model means nothing if the data don't
comply to the AR(2) model. The best you can hope for is that the
data respond well to an AR(P) order of 'sufficiently high' order.
If that is the case, you can inspect the pairs of roots and see
if any of them match the criteria you look for.

Rune