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

# System Identification

Started by ●June 30, 2008

Reply by ●June 30, 20082008-06-30

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

Reply by ●June 30, 20082008-06-30

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

Reply by ●June 30, 20082008-06-30

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

Reply by ●June 30, 20082008-06-30

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.comSample rate dependency implies you might have some MA in the system. Additive measurement noise also makes AR systems ARMA.

Reply by ●June 30, 20082008-06-30

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

Reply by ●June 30, 20082008-06-30

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

Reply by ●July 1, 20082008-07-01

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

Reply by ●July 1, 20082008-07-01

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

Reply by ●July 1, 20082008-07-01

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