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')
System Identification
Started by ●June 30, 2008
Reply by ●July 2, 20082008-07-02
Reply by ●July 2, 20082008-07-02
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.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 ConsultanMatrix 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
Reply by ●July 2, 20082008-07-02
Rune Allnor <allnor@tele.ntnu.no> wrote in news:5d66b163-db36-43c5-926e- f7d6b9e2971c@a70g2000hsh.googlegroups.com:> > 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
Reply by ●July 2, 20082008-07-02
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/HarminvYou 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
Reply by ●July 2, 20082008-07-02
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
Reply by ●July 8, 20082008-07-08
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
Reply by ●July 9, 20082008-07-09
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
Reply by ●July 9, 20082008-07-09
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.> > RuneVladimir Vassilevsky DSP and Mixed Signal Design Consultant http://www.abvolt.com
Reply by ●July 9, 20082008-07-09
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.comIf 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.
Reply by ●July 9, 20082008-07-09
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
