The FDLS worked wonderfully.
Pretty accurate results from the first try.

VLV





Rune Allnor wrote:

> On 9 Jul, 23:05, Vladimir Vassilevsky <antispam_bo...@hotmail.com>
> wrote:
> 
>>Rune Allnor wrote:
>>
>>>>This is what puzzled me. If one pair of poles is clearly dominating,
>>>>then why AR(2) couldn't find it in one shot.
>>
>>>The data you posted were 100
> 
> 
> Ouch! Should be 1000, not 100...
> 
> 
>>>samples long (or theerabout). The order
>>>estimator I used was of order 10. Computing the normal equations
>>>means that it uses the average of 100 independent spectra.
>>
>>>The impulse transient dominates the first ~200 samples, wheras the
>>>spectrum lines (power line noise?) remains for the duration of the
>>>recording. So all in all, the energy seems to be 'similar' even
>>>if the amplitudes are vastly different.
>>>If you want the f=0.2 lines, let the recording go on for 10 times
>>>as long and the AR(2) model is more likely to home in on them.
>>>If you want the start transient, truncate the recording earlier
>>>and the AR(2) model might home in on the dominant poles.
>>
>>This is not it.
> 
> 
> Yes it is, but there is more to it than mere order. For shorter
> data sequences the estimator for the covariance matrix plays a part,
> in addition to the data not cmplying perfectly to an AR(P) model.
> If you use the biased estimator, known from the Yule-Walker equations
> and the 100 first samples, the 4th order poles are
> 
>   0.90483768854503
>   0.78327057533512
>  -0.45960883639968 - 0.36140355237584i
>  -0.45960883639968 + 0.36140355237584i
> 
> which are way off the mark. If one uses the 4th order Tufts/Kumaresan
> method on the 100 first samples, the poles are
> 
>   0.95228809422338 - 0.11576998695318i
>   0.95228809422338 + 0.11576998695318i
>  -0.66220147395015
>  -0.16710812933497.
> 
> If you use the 10th order YW method on 1000 samples you will find
> a pair of poles at
> 
>   0.95387008130196 - 0.11751711909314i
>   0.95387008130196 + 0.11751711909314i
> 
> and
> 
>   0.95649011702604 - 0.11662905371121i
>   0.95649011702604 + 0.11662905371121i
> 
> with the T&K method.
> 
> You get well within 2 significant digits on the three latter methods
> (T&K, 100 samples / 1000 samples; YW 1000 samples). Which result one
> chooses to use is as much a matter of personal preference as anything
> else.
> 
> Rune

On 9 Jul, 23:05, Vladimir Vassilevsky <antispam_bo...@hotmail.com>
wrote:
> Rune Allnor wrote:
> >>This is what puzzled me. If one pair of poles is clearly dominating,
> >>then why AR(2) couldn't find it in one shot.
>
> > The data you posted were 100

Ouch! Should be 1000, not 100...

> > samples long (or theerabout). The order
> > estimator I used was of order 10. Computing the normal equations
> > means that it uses the average of 100 independent spectra.
>
> > The impulse transient dominates the first ~200 samples, wheras the
> > spectrum lines (power line noise?) remains for the duration of the
> > recording. So all in all, the energy seems to be 'similar' even
> > if the amplitudes are vastly different.
> > If you want the f=0.2 lines, let the recording go on for 10 times
> > as long and the AR(2) model is more likely to home in on them.
> > If you want the start transient, truncate the recording earlier
> > and the AR(2) model might home in on the dominant poles.
>
> This is not it.

Yes it is, but there is more to it than mere order. For shorter
data sequences the estimator for the covariance matrix plays a part,
in addition to the data not cmplying perfectly to an AR(P) model.
If you use the biased estimator, known from the Yule-Walker equations
and the 100 first samples, the 4th order poles are

  0.90483768854503
  0.78327057533512
 -0.45960883639968 - 0.36140355237584i
 -0.45960883639968 + 0.36140355237584i

which are way off the mark. If one uses the 4th order Tufts/Kumaresan
method on the 100 first samples, the poles are

  0.95228809422338 - 0.11576998695318i
  0.95228809422338 + 0.11576998695318i
 -0.66220147395015
 -0.16710812933497.

If you use the 10th order YW method on 1000 samples you will find
a pair of poles at

  0.95387008130196 - 0.11751711909314i
  0.95387008130196 + 0.11751711909314i

and

  0.95649011702604 - 0.11662905371121i
  0.95649011702604 + 0.11662905371121i

with the T&K method.

You get well within 2 significant digits on the three latter methods
(T&K, 100 samples / 1000 samples; YW 1000 samples). Which result one
chooses to use is as much a matter of personal preference as anything
else.

Rune

Rune Allnor wrote:

>>This is what puzzled me. If one pair of poles is clearly dominating,
>>then why AR(2) couldn't find it in one shot.
> 
> The data you posted were 100 samples long (or theerabout). The order
> estimator I used was of order 10. Computing the normal equations
> means that it uses the average of 100 independent spectra.
> 
> The impulse transient dominates the first ~200 samples, wheras the
> spectrum lines (power line noise?) remains for the duration of the
> recording. So all in all, the energy seems to be 'similar' even
> if the amplitudes are vastly different.
> If you want the f=0.2 lines, let the recording go on for 10 times
> as long and the AR(2) model is more likely to home in on them.
> If you want the start transient, truncate the recording earlier
> and the AR(2) model might home in on the dominant poles.


This is not it. I tried to use only the meaningful data at the 
beginning; I also tried to make the unbiased WSS data by combining the 
the shifted versions of the input with random amplitude and phase. The 
AR(2) result is no different. However if the data is decimated by the 
factor of 5..15, then AR(2) works quite well.

Here is why:

AR tries to minimize the sum of errors at each step:

E[i] = x[i] - A1*x[i-1] - A2*x[i-2] - ....

This is the different problem from finding the best approximation of the 
impulse response! Those problems are identical only in the case of the 
ideal data.

Besides, the mean difference between x[i] and x[i-1] is at the order of 
Amplitude * Fc/Fsa. This makes the model very sensitive to imperfections.

> In between, well, use a 'reasonably' high AR(P) model and
> see what happens.

Higher order models appear to be not too diferent from AR(2):
The dominating poles are still real, however the next pair of poles is 
close to what is expected.


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

On 9 Jul, 16:41, Vladimir Vassilevsky <antispam_bo...@hotmail.com>
wrote:
> Rune Allnor wrote:
> > 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.
>
> NOW I understand what you are saying! Thank you for pointing me.

Y're welcome.

> > Sure. One pair of poles completely domminates the impulse response.
>
> This is what puzzled me. If one pair of poles is clearly dominating,
> then why AR(2) couldn't find it in one shot.

The data you posted were 100 samples long (or theerabout). The order
estimator I used was of order 10. Computing the normal equations
means that it uses the average of 100 independent spectra.

The impulse transient dominates the first ~200 samples, wheras the
spectrum lines (power line noise?) remains for the duration of the
recording. So all in all, the energy seems to be 'similar' even
if the amplitudes are vastly different.

If you want the f=0.2 lines, let the recording go on for 10 times
as long and the AR(2) model is more likely to home in on them.
If you want the start transient, truncate the recording earlier
and the AR(2) model might home in on the dominant poles.

In between, well, use a 'reasonably' high AR(P) model and
see what happens.

> > 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.
>
> This implies the numeric solving of the polynomial, and estimation of
> the amplitudes. I.e. the "full" Prony method. I wonder if Greg Berchin
> method of FDLS can be a quicker way to solution; it seems like the
> selection of the reference points is clear in this case.

Don't know; haven't tried it. For what it's worth, the full
Prony's method is trivially implemented in matlab.

Rune

On Jul 8, 10:01 am, 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.
>
> 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 Consultanthttp://www.abvolt.com

If I recall correctly, assuming you used the impulse matching routine
in the Signal processing toolbox in matlab, the impulse matching
method due to Friedlander is full ARMA.

Rune Allnor wrote:

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

NOW I understand what you are saying! Thank you for pointing me.

The LPC model can be surprisingly wrong even if the data mismatch is 
small. This happens due to the nature of the LPC method. The way to deal 
with the problem is increasing the order of the model. It took me a 
whole day to realize that.

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

This is what puzzled me. If one pair of poles is clearly dominating, 
then why AR(2) couldn't find it in one shot.

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

Yes, yes, and yes.

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

This implies the numeric solving of the polynomial, and estimation of 
the amplitudes. I.e. the "full" Prony method. I wonder if Greg Berchin 
method of FDLS can be a quicker way to solution; it seems like the 
selection of the reference points is clear in this case.

> 
> Rune

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

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

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