COMPGROUPS.NET | Search | Post Question | Groups | Stream | About | Register

### AR coefficients for complex numbers

• Follow

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.

 0
Reply sharathreddy 1/16/2004 7:57:00 PM

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

 0
Reply allnor 1/17/2004 12:27:52 PM

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

 0
Reply sharathreddy 1/19/2004 1:46:59 AM

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

 0
Reply allnor 1/19/2004 10:58:10 AM

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

 0
Reply sharathreddy 1/19/2004 10:25:08 PM

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

 0
Reply allnor 1/20/2004 6:21:04 AM

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

 0
Reply sharathreddy 1/20/2004 10:49:40 PM

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

 0
Reply allnor 1/21/2004 8:15:12 AM

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

 0
Reply allnor 1/21/2004 8:17:14 AM

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

 0
Reply sharathreddy 1/22/2004 5:04:30 AM

9 Replies
237 Views

(page loaded in 0.16 seconds)

Similiar Articles:

7/26/2012 11:22:35 PM