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

### AR coefficients for complex numbers

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

See related articles to this posting

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

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

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

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

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

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

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

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

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

9 Replies
305 Views

Similar Articles

12/9/2013 4:48:05 AM
[PageSpeed]

Similar Artilces:

Checking for invalid floating point numbers
I'm using Borland Turbo C++ 3.0 to develop an embedded system to shift data around a network. At the moment we receive a string of bytes over a serial line and reassemble them into floating point values. If the bytes are not assembled correctly then it's possible to produce some floating point values that aren't 'genuine' numbers. Does anyone have any suggestions on how to check a float/double is valid (i.e. checking for NANs, infinity etc.)? Turbo C++ doesn't seem to have much support for this. Thanks, -- L. LSW wrote: > I'm using Borland Turbo C++ 3.0 to develop an embedded system to shift > data around a network. At the moment we receive a string of bytes > over a serial line and reassemble them into floating point values. > > If the bytes are not assembled correctly then it's possible to produce > some floating point values that aren't 'genuine' numbers. Does anyone > have any suggestions on how to check a float/double is valid (i.e. > checking for NANs, infinity etc.)? Implementation-specific stuff. Ask in the newsgroup for your compiler. The IEEE 754 specification leaves very little

Convert words to numbers and back?
I was wondering if somebody could give me some insight and help on how to convert words to numbers and back again. Maybe a function? -- Posted via http://www.ruby-forum.com/. Note: parts of this message were removed by the gateway to make it a legal Usenet post. Check http://www.deveiate.org/projects/Linguistics/ I think numwords is what you are looking for. On 11/24/07, Jordon Bedwell <jordon@envygeeks.com> wrote: > > I was wondering if somebody could give me some insight and help on how > to convert words to numbers and back again. Maybe a function? > -- > Posted... are ignored. If there is not a valid number at the start of str, 0.0 is returned. This method never raises an exception. "123.45e1".to_f #=> 1234.5 "45.67 degrees".to_f #=> 45.67 "thx1138".to_f #=> 0.0 bye sala----- Original Message ----- From: "Andrei Maxim" <andrei@andreimaxim.ro> To: "ruby-talk ML" <ruby-talk@ruby-lang.org> Sent: Saturday, November 24, 2007 2:05 PM Subject: Re: Convert words to numbers and back? > Check http://www.deveiate.org/projects/Linguistics/ > I think numwords

complex matrix inversion
matrix takes more than 2ms) Can any o= ne > tell me free intel library or any other related library to do my task mor= e > fastly? You are in the right neighbourhood. My matlab inverts one complex-valued 52x52 matrix (looking at the other numbers, do you mean 32x32?) in 0.5 ms on the 2.6 GHz PC. It *might* be that this is a parallel inversion routine which takes advantage of the fact that my PC has two cores (I haven't explicitly configured matlab to do that, though), but it seems you are within a factor 5-10 of the libraries matlab uses for these things, depending on the system...Hello I need to invert square complex matrixes of sizes 52x52,64x64 and 128x128 in MS VC++6.I need to implement it as fast as possible. I have written a code for complex matrix inversion but it takes lot of time. (e.g. a complex matrix inversion of 52x52 matrix takes more than 2ms) Can any one tell me free intel library or any other related library to do my task more fastly? THANKS IN ADVANCE On Nov 21, 11:43=A0am, "asimmasud" <asimbhatti...@hotmail.com> wrote: > Hello > > =A0 =A0 =A0 I need to invert

Re: missing numeric 'smallest numeric'? #13 660777
On Wed, 8 Feb 2006 21:06:29 -0800, Jack Hamilton <jfh@ACM.ORG> wrote: >on 2/7/2006 5:33 AM Ed Heaton said the following: > >> SAS numerics take advantage of the fact that some bit combinations in >> the IEEE eight-byte, double-precision standard are not numbers. > >It's a bit more complicated than that - neither VAX machines nor OS/360 >machines use IEEE representation, and SAS has missing values on those >platforms. Also, a missing value can be stored in 2-7 bytes as well as >8 bytes. And probably 16 bytes when we make the next leap forward. > >It may be more accurate to say that missing values are stored as bit >patterns which do not represent floating point numbers. > >Many years ago, the missing value representation was described in the >SAS/Toolkit documentation, but if it's online I can't find it. Not that I would want to suggest using the internal representation fo missing values, but here is a short demo showing how you might discover, by using peekc. 880 data _null_; 881 length n 8 real $8 ; 882 format real$hex16. ; 883 do n = ._, ., .a, .b, .z, 0, 2/32767; 884 real = peekc( addr

where to buy decent laptop without complex rebate ?
My mother bought a PC from best buy and I think the rebate paper work was worse than fileing your tax returns. Any recommendations on how you can get a decent laptop at a good price without going through that ? I'm mainly looking to go wireless on the internet, do some perl, java programming, web development, and install perhaps install linux, perhaps debian. Any recommendations on a brand or machine ? I saw this, a presario, for \$500 something like that sounds not bad with 40 gig, but other systems have more disk, not sure if more disk might be better down the road. http://www.comp