System Identification

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
> &#4294967295; 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&#4294967295;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&#4294967295;pm, Vladimir Vassilevsky <antispam_bo...@hotmail.com>
wrote:
> Hello All,
>
> There is a physical system which is described well by the 2nd order pole
> &#4294967295; 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=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