Hello, Can anyone point me to a MATLAB function/program that calculates the AR coefficients for a complex series of numbers, similar to 'aryule' for a real numbers? Thanks a lot all. Sharath.
AR coefficients for complex numbers
Started by ●January 16, 2004
Reply by ●January 17, 20042004-01-17
sharathreddy@hotmail.com (Sharath) wrote in message news:<3c62a40e.0401161157.817a14d@posting.google.com>...> Hello, > > Can anyone point me to a MATLAB function/program that calculates the > AR coefficients for a complex series of numbers, similar to 'aryule' > for a real numbers? > > Thanks a lot all. > > > Sharath.Unfortunately I don't have the ARYULE function available so I can't look into it and see what it does. Still, I find it hard to believe that it does not handle complex data. the matlab "'" (apostrophe) operator transposes *and* conjugates the vector. The matlab ".'" (point and apostrophe) operator *only* transposes the vector. If you find the ".'" anywhere in ARYULE, copy ARYULE.M to MYARYULE.M and change ".'" to "'". Now, the problematic part is if you find any calls to FLIPUD im MYARYULE. Chances are that at least some (although not necessarily all) FLIPUDs should be wrapped in a CONJ. So, in effect, make your own function more or less as follows: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function car=cplxar(x,p) % car=cplxar(x,p) % % Computes complex-valued AR(p) coefficients from data. % % Input parameters: % % x - Data column vector % p - Number of AR coefficients % % Output: % % car - Complex AR coefficients n=length(x); X=zeros(n+p-1,p+1); for m=1:p+1 X(m:m+n-1)=x; end Rxx=X'*X; c=-inv(Rxx(2:p+1,2:p+1))*Rxx(2:p+1,1); car=[1;c]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Or something like that. You could also check table 8.2 in Therrien: Discrete Random Signals and Statistical Signal Processing Prentice-Hall, 1992 if you want to work your way through all the nitty gritty details. Rune
Reply by ●January 18, 20042004-01-18
Rune- Thank you very much for your assistance. Now I have it working for complex numbers too. The code was indeed helpful. I would appreciate it if you/anyone could give me some pointers on another issue I am facing. Now, I need to use the AIC to try and find out the optimal order of a complex AR time series. Matlab literature suggests I use the 'arxstruc' function, but I am having trouble understanding the actual implementation. The data I have is: x = I/P time series AR coefficients for different orders using aryule y = Constructed model using these different order coefficients. The function seems to have an argument called 'delay', which I do not know what value to set. All I am trying to do is trying to find the optimal order of a given time series of type AR(n) using AIC. Any help would be greatly appreciated. Thank you very much. Regards, Sharath. allnor@tele.ntnu.no (Rune Allnor) wrote in message news:<f56893ae.0401170427.714ea053@posting.google.com>...> sharathreddy@hotmail.com (Sharath) wrote in message news:<3c62a40e.0401161157.817a14d@posting.google.com>... > > Hello, > > > > Can anyone point me to a MATLAB function/program that calculates the > > AR coefficients for a complex series of numbers, similar to 'aryule' > > for a real numbers? > > > > Thanks a lot all. > > > > > > Sharath. > > Unfortunately I don't have the ARYULE function available so I can't > look into it and see what it does. Still, I find it hard to believe > that it does not handle complex data. the matlab "'" (apostrophe) > operator transposes *and* conjugates the vector. The matlab ".'" > (point and apostrophe) operator *only* transposes the vector. > If you find the ".'" anywhere in ARYULE, copy ARYULE.M to MYARYULE.M > and change ".'" to "'". > > Now, the problematic part is if you find any calls to FLIPUD im MYARYULE. > Chances are that at least some (although not necessarily all) FLIPUDs > should be wrapped in a CONJ. > > So, in effect, make your own function more or less as follows: > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > function car=cplxar(x,p) > % car=cplxar(x,p) > % > % Computes complex-valued AR(p) coefficients from data. > % > % Input parameters: > % > % x - Data column vector > % p - Number of AR coefficients > % > % Output: > % > % car - Complex AR coefficients > > n=length(x); > X=zeros(n+p-1,p+1); > for m=1:p+1 > X(m:m+n-1)=x; > end > > Rxx=X'*X; > c=-inv(Rxx(2:p+1,2:p+1))*Rxx(2:p+1,1); > car=[1;c]; > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > Or something like that. > > You could also check table 8.2 in > > Therrien: Discrete Random Signals and Statistical Signal Processing > Prentice-Hall, 1992 > > if you want to work your way through all the nitty gritty details. > > Rune
Reply by ●January 19, 20042004-01-19
sharathreddy@hotmail.com (Sharath) wrote in message news:<3c62a40e.0401181746.678993b2@posting.google.com>...> Rune- > > Thank you very much for your assistance. Now I have it working for > complex numbers too. The code was indeed helpful. > > I would appreciate it if you/anyone could give me some pointers on > another issue I am facing. Now, I need to use the AIC to try and find > out the optimal order of a complex AR time series. Matlab literature > suggests I use the 'arxstruc' function, but I am having trouble > understanding the actual implementation.Ouch! I just hate it when this happens! Matlab has chosen a very bad way of implementing these types of functions. As you can see from the code snippet I posted, you have already chosen the model order when you call my function CPLXAR. The same applies for the matlab codes. Now, if you want to do these sorts of things "for real", you need to implement "the real thing", i.e. the Levinson recursion. If I knew that the Therrien (reference below) book was easy to find, I'd say "get your hands on a copy and implement the pseudocode of table 8.2". Unfortunately, I know that the book is very hard to find, so I can't tell you to do that [*]. In this hypothetical situation where you *did* have Therrien available, I would advice you to insert the AIC somewhere near the end of the loop over order but before the order update, and use the reflection coefficients that pop out of the Levinson recursion to check whether the AIC allows the order to actually increase, or if the next order is a worse model fit than the current. Therrien or no Therrien, I guess the basic message is "ditch the canned matlab routines and implement the stuff yourself". Rune [*] I just heard from the teacher who uses the Therrien book in a statistichal SP course here, that the publisher (Prentice-Hall) don't have any plans to reprint the book. The reason is, alledgedly, that they have registered no demand for the book. If you (or any others) try to get it through your local bookstores, make sure that the bookstores actually contact the publisher and ask what plans there are for a re-issue, and not only searches a database where they will find that the book is out of print.> The data I have is: > x = I/P time series > AR coefficients for different orders using aryule > y = Constructed model using these different order coefficients. > > The function seems to have an argument called 'delay', which I do not > know what value to set. All I am trying to do is trying to find the > optimal order of a given time series of type AR(n) using AIC. > > Any help would be greatly appreciated. Thank you very much. > > Regards, > Sharath. > > allnor@tele.ntnu.no (Rune Allnor) wrote in message news:<f56893ae.0401170427.714ea053@posting.google.com>... > > sharathreddy@hotmail.com (Sharath) wrote in message news:<3c62a40e.0401161157.817a14d@posting.google.com>... > > > Hello, > > > > > > Can anyone point me to a MATLAB function/program that calculates the > > > AR coefficients for a complex series of numbers, similar to 'aryule' > > > for a real numbers? > > > > > > Thanks a lot all. > > > > > > > > > Sharath. > > > > Unfortunately I don't have the ARYULE function available so I can't > > look into it and see what it does. Still, I find it hard to believe > > that it does not handle complex data. the matlab "'" (apostrophe) > > operator transposes *and* conjugates the vector. The matlab ".'" > > (point and apostrophe) operator *only* transposes the vector. > > If you find the ".'" anywhere in ARYULE, copy ARYULE.M to MYARYULE.M > > and change ".'" to "'". > > > > Now, the problematic part is if you find any calls to FLIPUD im MYARYULE. > > Chances are that at least some (although not necessarily all) FLIPUDs > > should be wrapped in a CONJ. > > > > So, in effect, make your own function more or less as follows: > > > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > > > function car=cplxar(x,p) > > % car=cplxar(x,p) > > % > > % Computes complex-valued AR(p) coefficients from data. > > % > > % Input parameters: > > % > > % x - Data column vector > > % p - Number of AR coefficients > > % > > % Output: > > % > > % car - Complex AR coefficients > > > > n=length(x); > > X=zeros(n+p-1,p+1); > > for m=1:p+1 > > X(m:m+n-1)=x; > > end > > > > Rxx=X'*X; > > c=-inv(Rxx(2:p+1,2:p+1))*Rxx(2:p+1,1); > > car=[1;c]; > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > > > Or something like that. > > > > You could also check table 8.2 in > > > > Therrien: Discrete Random Signals and Statistical Signal Processing > > Prentice-Hall, 1992 > > > > if you want to work your way through all the nitty gritty details. > > > > Rune
Reply by ●January 19, 20042004-01-19
Rune- When you get the chance, could you check this website out. http://www.gps.caltech.edu/~tapio/arfit/ It's got an interesting function called 'arfit', which gives the AIC and SBC values for real values inputs. The issue really is with complex inputs, and there's only one function called 'arord' that gives me an error because of this complex input. It's because of a function called 'chol', which hiccups on me. When I googled on how to rectify this, the closest solution I got was to take the force the passed argument into a hermitian matrix (or inverse i dont remember) by performing (Rp+Rp')/2. Ofcourse I don't really know the details of why its being done, which is probably why its not working. Any ideas on this? Thanks a ton! I remain, Sharath. allnor@tele.ntnu.no (Rune Allnor) wrote in message news:<f56893ae.0401190258.114aa50f@posting.google.com>...> sharathreddy@hotmail.com (Sharath) wrote in message news:<3c62a40e.0401181746.678993b2@posting.google.com>... > > Rune- > > > > Thank you very much for your assistance. Now I have it working for > > complex numbers too. The code was indeed helpful. > > > > I would appreciate it if you/anyone could give me some pointers on > > another issue I am facing. Now, I need to use the AIC to try and find > > out the optimal order of a complex AR time series. Matlab literature > > suggests I use the 'arxstruc' function, but I am having trouble > > understanding the actual implementation. > > Ouch! I just hate it when this happens! Matlab has chosen a very > bad way of implementing these types of functions. As you can see > from the code snippet I posted, you have already chosen the model > order when you call my function CPLXAR. The same applies > for the matlab codes. > > Now, if you want to do these sorts of things "for real", you need > to implement "the real thing", i.e. the Levinson recursion. > If I knew that the Therrien (reference below) book was easy to find, > I'd say "get your hands on a copy and implement the pseudocode > of table 8.2". Unfortunately, I know that the book is very hard to > find, so I can't tell you to do that [*]. > > In this hypothetical situation where you *did* have Therrien available, > I would advice you to insert the AIC somewhere near the end of the > loop over order but before the order update, and use the reflection > coefficients that pop out of the Levinson recursion to check whether > the AIC allows the order to actually increase, or if the next order > is a worse model fit than the current. > > Therrien or no Therrien, I guess the basic message is "ditch the > canned matlab routines and implement the stuff yourself". > > Rune > > [*] I just heard from the teacher who uses the Therrien book in > a statistichal SP course here, that the publisher (Prentice-Hall) > don't have any plans to reprint the book. The reason is, alledgedly, > that they have registered no demand for the book. If you (or any > others) try to get it through your local bookstores, make sure > that the bookstores actually contact the publisher and ask what > plans there are for a re-issue, and not only searches a database > where they will find that the book is out of print. > > > The data I have is: > > x = I/P time series > > AR coefficients for different orders using aryule > > y = Constructed model using these different order coefficients. > > > > The function seems to have an argument called 'delay', which I do not > > know what value to set. All I am trying to do is trying to find the > > optimal order of a given time series of type AR(n) using AIC. > > > > Any help would be greatly appreciated. Thank you very much. > > > > Regards, > > Sharath. > > > > allnor@tele.ntnu.no (Rune Allnor) wrote in message news:<f56893ae.0401170427.714ea053@posting.google.com>... > > > sharathreddy@hotmail.com (Sharath) wrote in message news:<3c62a40e.0401161157.817a14d@posting.google.com>... > > > > Hello, > > > > > > > > Can anyone point me to a MATLAB function/program that calculates the > > > > AR coefficients for a complex series of numbers, similar to 'aryule' > > > > for a real numbers? > > > > > > > > Thanks a lot all. > > > > > > > > > > > > Sharath. > > > > > > Unfortunately I don't have the ARYULE function available so I can't > > > look into it and see what it does. Still, I find it hard to believe > > > that it does not handle complex data. the matlab "'" (apostrophe) > > > operator transposes *and* conjugates the vector. The matlab ".'" > > > (point and apostrophe) operator *only* transposes the vector. > > > If you find the ".'" anywhere in ARYULE, copy ARYULE.M to MYARYULE.M > > > and change ".'" to "'". > > > > > > Now, the problematic part is if you find any calls to FLIPUD im MYARYULE. > > > Chances are that at least some (although not necessarily all) FLIPUDs > > > should be wrapped in a CONJ. > > > > > > So, in effect, make your own function more or less as follows: > > > > > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > > > > > function car=cplxar(x,p) > > > % car=cplxar(x,p) > > > % > > > % Computes complex-valued AR(p) coefficients from data. > > > % > > > % Input parameters: > > > % > > > % x - Data column vector > > > % p - Number of AR coefficients > > > % > > > % Output: > > > % > > > % car - Complex AR coefficients > > > > > > n=length(x); > > > X=zeros(n+p-1,p+1); > > > for m=1:p+1 > > > X(m:m+n-1)=x; > > > end > > > > > > Rxx=X'*X; > > > c=-inv(Rxx(2:p+1,2:p+1))*Rxx(2:p+1,1); > > > car=[1;c]; > > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > > > > > Or something like that. > > > > > > You could also check table 8.2 in > > > > > > Therrien: Discrete Random Signals and Statistical Signal Processing > > > Prentice-Hall, 1992 > > > > > > if you want to work your way through all the nitty gritty details. > > > > > > Rune
Reply by ●January 20, 20042004-01-20
sharathreddy@hotmail.com (Sharath) wrote in message news:<3c62a40e.0401191425.53b828d9@posting.google.com>...> Rune- > > When you get the chance, could you check this website out. > http://www.gps.caltech.edu/~tapio/arfit/ > > It's got an interesting function called 'arfit', which gives the AIC > and SBC values for real values inputs. The issue really is with > complex inputs, and there's only one function called 'arord' that > gives me an error because of this complex input. It's because of a > function called 'chol', which hiccups on me. When I googled on how to > rectify this, the closest solution I got was to take the force the > passed argument into a hermitian matrix (or inverse i dont remember) > by performing (Rp+Rp')/2. Ofcourse I don't really know the details of > why its being done, which is probably why its not working. > > Any ideas on this? Thanks a ton! > > I remain, > Sharath.Well, according to the help file of CHOL, there are no restrictions on the matrix being only real-valued. The Cholesky decomposition does, however, require the matrix to be positive definite. Which means you have to be quite careful with how you construct your data when testing those functions. An AR-fitting method that uses CHOL would, unless very careful precautions were taken, fail on numerical grounds if it tried to fit an AR(2) model to noise-free AR(1) data. Furthermore, my first objection about methods that require pre-determined orders still hold with this method. I'm not sure I want to comment much further than that. There is an Arnold Neumaier who is very active on the NG sci.math.num-analysis. I think it's a safe guess that this is the same person who co-wrote the ARfit package. You should really check with him personally if you want to pursue the ARfit package. Rune
Reply by ●January 20, 20042004-01-20
Rune, I agree. The inputs need to be very carefull chosen else I get very ugly results. However, after digging more into model identification in Matlab I came across this interesting time series course website. http://www.control.lth.se/~kurspi/local/ It has a pdf file called extra exercises and solutions. On page 16, it has worked out an implementation of system identification (the question is on page 5). Although its not the exact issue I am facing, it seems quite close. Unfortunately I still am not able to get it working to my satisfaction. A lack of my understanding perhaps? I'm interested in your comments about the usage of the specified functions. Thanks a ton for your time and effort in replying. I really appreciate it. I remain, Sharath. allnor@tele.ntnu.no (Rune Allnor) wrote in message news:<f56893ae.0401192221.768a6366@posting.google.com>...> sharathreddy@hotmail.com (Sharath) wrote in message news:<3c62a40e.0401191425.53b828d9@posting.google.com>... > > Rune- > > > > When you get the chance, could you check this website out. > > http://www.gps.caltech.edu/~tapio/arfit/ > > > > It's got an interesting function called 'arfit', which gives the AIC > > and SBC values for real values inputs. The issue really is with > > complex inputs, and there's only one function called 'arord' that > > gives me an error because of this complex input. It's because of a > > function called 'chol', which hiccups on me. When I googled on how to > > rectify this, the closest solution I got was to take the force the > > passed argument into a hermitian matrix (or inverse i dont remember) > > by performing (Rp+Rp')/2. Ofcourse I don't really know the details of > > why its being done, which is probably why its not working. > > > > Any ideas on this? Thanks a ton! > > > > I remain, > > Sharath. > > Well, according to the help file of CHOL, there are no restrictions on > the matrix being only real-valued. The Cholesky decomposition does, > however, require the matrix to be positive definite. Which means you > have to be quite careful with how you construct your data when testing > those functions. An AR-fitting method that uses CHOL would, unless very > careful precautions were taken, fail on numerical grounds if it tried to > fit an AR(2) model to noise-free AR(1) data. Furthermore, my first > objection about methods that require pre-determined orders still hold > with this method. > > I'm not sure I want to comment much further than that. There is an > Arnold Neumaier who is very active on the NG sci.math.num-analysis. > I think it's a safe guess that this is the same person who co-wrote > the ARfit package. You should really check with him personally if > you want to pursue the ARfit package. > > Rune
Reply by ●January 21, 20042004-01-21
sharathreddy@hotmail.com (Sharath) wrote in message news:<3c62a40e.0401201449.5f766d67@posting.google.com>...> Rune, > > I agree. The inputs need to be very carefull chosen else I get very > ugly results. However, after digging more into model identification in > Matlab I came across this interesting time series course website. > > http://www.control.lth.se/~kurspi/local/ > > It has a pdf file called extra exercises and solutions. On page 16, it > has worked out an implementation of system identification (the > question is on page 5).Eh, I'm not sure I find the correct pdf file.> Although its not the exact issue I am facing, > it seems quite close. Unfortunately I still am not able to get it > working to my satisfaction. A lack of my understanding perhaps? I'm > interested in your comments about the usage of the specified > functions. Thanks a ton for your time and effort in replying. I really > appreciate it. > > I remain, > Sharath.Since I can't point you to other sources to the levinson recursion than Therrien's book, I'll try to reproduce the infamous table 8.2. Make sure you view this in a fixed-width font. I use LaTeX notation to denote maths symbols. Superscript "H" denotes the conjugates and transposed of vectors, superscript "*" denotes complex conjugation only, while superscript "T" denotes transposed only. Capital letters denote column vector quantities, and the product JA where A is a vector corresponds to the matlab function "flipud". (J is the reflection matrix where all entries are zero, except for on the anti diagonal, where they are 1). If you want to dig into the gory details, note that the quantities that are marked with a prime "'" relate to the backward prediction. ----------------------------------------------------------------------------- Initialization: R_0 = r_xx(1) A_0 =B_0 =1 \sigma_0^2 = \sigma'_0^2 = r_xx(0) I ) \Delta_p=R^H_{p-1}JA_{p-1} \Delta'_p=R^H_{p-1}JB_{p-1} \Delta_p II) \gamma_p = --------------- \sigma^2_{p-1} \Delta'_p \gamma'_p = ----------------- \sigma'^2_{p-1} | A_{p-1} | | 0 | III) A_p = | | -\gamma_p | | | 0 | | JB_{p-1} | | B_{p-1} | | 0 | B_p = | | -\gamma'_p | | | 0 | | JA_{p-1} | IV) \sigma^2_p = \sigma^2_{p-1} - \gamma_p\Delta'_p \sigma'^2_p = \sigma'^2_{p-1} - \gamma'_p\Delta_p -------------------------------------------------------------------- The only vector/matrix quantities here are A, B, R and J. A is the vector of forward AR coefficients, B is the vector of backwards AR coefficients, R is the vector of complex-valued autocorrelation coefficients from the signal, and J is the reflection matrix. All other quantities are scalars. The place for the AIC order estimator would be between steps II and III in this scheme. If the AIC, which test the order estimates on the basis of reflection coefficients \gamma, (or is it |\gamma|^2) finds that there is a penalty in increasing the order from p-1 to p, the algorithm breaks after step II and pop out A_{p-1} at the output. If there is no penalty in increasing the order, the iteration takes one more turn. Rune
Reply by ●January 21, 20042004-01-21
sharathreddy@hotmail.com (Sharath) wrote in message news:<3c62a40e.0401201449.5f766d67@posting.google.com>...> Rune, > > I agree. The inputs need to be very carefull chosen else I get very > ugly results. However, after digging more into model identification in > Matlab I came across this interesting time series course website. > > http://www.control.lth.se/~kurspi/local/ > > It has a pdf file called extra exercises and solutions. On page 16, it > has worked out an implementation of system identification (the > question is on page 5).Eh, I'm not sure I find the correct pdf file.> Although its not the exact issue I am facing, > it seems quite close. Unfortunately I still am not able to get it > working to my satisfaction. A lack of my understanding perhaps? I'm > interested in your comments about the usage of the specified > functions. Thanks a ton for your time and effort in replying. I really > appreciate it. > > I remain, > Sharath.Since I can't point you to other sources to the levinson recursion than Therrien's book, I'll try to reproduce the infamous table 8.2. Make sure you view this in a fixed-width font. I use LaTeX notation to denote maths symbols. Superscript "H" denotes the conjugates and transposed of vectors, superscript "*" denotes complex conjugation only, while superscript "T" denotes transposed only. Capital letters denote column vector quantities, and the product JA where A is a vector corresponds to the matlab function "flipud". (J is the reflection matrix where all entries are zero, except for on the anti diagonal, where they are 1). If you want to dig into the gory details, note that the quantities that are marked with a prime "'" relate to the backward prediction. ----------------------------------------------------------------------------- Initialization: R_0 = r_xx(1) A_0 =B_0 =1 \sigma_0^2 = \sigma'_0^2 = r_xx(0) I ) \Delta_p=R^H_{p-1}JA_{p-1} \Delta'_p=R^H_{p-1}JB_{p-1} \Delta_p II) \gamma_p = --------------- \sigma^2_{p-1} \Delta'_p \gamma'_p = ----------------- \sigma'^2_{p-1} | A_{p-1} | | 0 | III) A_p = | | -\gamma_p | | | 0 | | JB_{p-1} | | B_{p-1} | | 0 | B_p = | | -\gamma'_p | | | 0 | | JA_{p-1} | IV) \sigma^2_p = \sigma^2_{p-1} - \gamma_p\Delta'_p \sigma'^2_p = \sigma'^2_{p-1} - \gamma'_p\Delta_p -------------------------------------------------------------------- The only vector/matrix quantities here are A, B, R and J. A is the vector of forward AR coefficients, B is the vector of backwards AR coefficients, R is the vector of complex-valued autocorrelation coefficients from the signal, and J is the reflection matrix. All other quantities are scalars. The place for the AIC order estimator would be between steps II and III in this scheme. If the AIC, which test the order estimates on the basis of reflection coefficients \gamma, (or is it |\gamma|^2) finds that there is a penalty in increasing the order from p-1 to p, the algorithm breaks after step II and pop out A_{p-1} at the output. If there is no penalty in increasing the order, the iteration takes one more turn. Rune
Reply by ●January 22, 20042004-01-22
Rune, The ARfit package works if you change the complex input from z(t)=x(t)+i*y(t) to z(t)=(x(t)^T,y(t)^T)^T; I will try and implement your method in Matlab too. Thanks a ton! Regards, Sharath. allnor@tele.ntnu.no (Rune Allnor) wrote in message news:<f56893ae.0401210017.16425499@posting.google.com>...> sharathreddy@hotmail.com (Sharath) wrote in message news:<3c62a40e.0401201449.5f766d67@posting.google.com>... > > Rune, > > > > I agree. The inputs need to be very carefull chosen else I get very > > ugly results. However, after digging more into model identification in > > Matlab I came across this interesting time series course website. > > > > http://www.control.lth.se/~kurspi/local/ > > > > It has a pdf file called extra exercises and solutions. On page 16, it > > has worked out an implementation of system identification (the > > question is on page 5). > > Eh, I'm not sure I find the correct pdf file. > > > Although its not the exact issue I am facing, > > it seems quite close. Unfortunately I still am not able to get it > > working to my satisfaction. A lack of my understanding perhaps? I'm > > interested in your comments about the usage of the specified > > functions. Thanks a ton for your time and effort in replying. I really > > appreciate it. > > > > I remain, > > Sharath. > > Since I can't point you to other sources to the levinson recursion than > Therrien's book, I'll try to reproduce the infamous table 8.2. > > Make sure you view this in a fixed-width font. I use LaTeX notation > to denote maths symbols. > > Superscript "H" denotes the conjugates and transposed of vectors, > superscript "*" denotes complex conjugation only, while superscript > "T" denotes transposed only. Capital letters denote column vector > quantities, and the product JA where A is a vector corresponds to the > matlab function "flipud". (J is the reflection matrix where all entries > are zero, except for on the anti diagonal, where they are 1). > > If you want to dig into the gory details, note that the quantities > that are marked with a prime "'" relate to the backward prediction. > > ----------------------------------------------------------------------------- > Initialization: > > R_0 = r_xx(1) > A_0 =B_0 =1 > \sigma_0^2 = \sigma'_0^2 = r_xx(0) > > I ) \Delta_p=R^H_{p-1}JA_{p-1} > \Delta'_p=R^H_{p-1}JB_{p-1} > > \Delta_p > II) \gamma_p = --------------- > \sigma^2_{p-1} > > \Delta'_p > \gamma'_p = ----------------- > \sigma'^2_{p-1} > > > | A_{p-1} | | 0 | > III) A_p = | | -\gamma_p | | > | 0 | | JB_{p-1} | > > > | B_{p-1} | | 0 | > B_p = | | -\gamma'_p | | > | 0 | | JA_{p-1} | > > > IV) \sigma^2_p = \sigma^2_{p-1} - \gamma_p\Delta'_p > \sigma'^2_p = \sigma'^2_{p-1} - \gamma'_p\Delta_p > > -------------------------------------------------------------------- > > The only vector/matrix quantities here are A, B, R and J. > > A is the vector of forward AR coefficients, B is the vector of > backwards AR coefficients, R is the vector of complex-valued > autocorrelation coefficients from the signal, and J is the > reflection matrix. All other quantities are scalars. > > The place for the AIC order estimator would be between steps II > and III in this scheme. If the AIC, which test the order estimates > on the basis of reflection coefficients \gamma, (or is it |\gamma|^2) > finds that there is a penalty in increasing the order from p-1 to p, > the algorithm breaks after step II and pop out A_{p-1} at the output. > If there is no penalty in increasing the order, the iteration takes > one more turn. > > Rune
