[comp.dsp] AR and MA questions

rob@kaa.eng.ohio-state.edu (Rob Carriere) (11/09/89)

In article <1389@mrsvr.UUCP> kohli@gemed.ge.com (Mr. Bad Judgment) writes:
>
>I've been working a little bit with some ARMA stuff
>and have run into some questions which I can't
>generalize an answer to from my experience:
>
>1. Are the AR "parameters" the poles of the waveform?

Not quite.  The AR parameters are the coefficients of a polynomial whose roots
are the poles.  Note that for an AR(1) process it doesn't matter (root = coeff)

>	If so, does this mean that an AR parameter of magnitude
>	> 1 indicates an unstable system (assuming no zero to
>	counteract it)?

No, an _AR_coefficient_ >1 is quite innocent.  It is the poles you have to
worry about. 

>2. (And the MA parameters *are* the zeroes?)

Same thing here: the zeros are the roots of the polynomial formed by the MA
coefficients.  Again, for MA(1), it doesn't matter.

>3. Is the noise process usually neglected when
>	using the AR and MA parameters to model the
>	time series?  If not, is a Gaussian model

The underlying model is:

                  ______________
                 |              |
white noise ---> |  ARMA Filter | ---> observed signal
                 |              |
                  --------------

Here ``ARMA Filter'' means the filter with the transfer function given by the
estimated ARMA parameters.  The white noise is usually assumed to be Gaussian.

What you seem to be looking for is:              Noise
                                 _______________   |
                                |               |  |
Deterministic white signal ---> |  Prony Filter | -+-> observed signal
                                |               |
                                 ---------------

This is handled in the literature as ``Prony Estimation''  (and is very
similar to ARMA stuff).  The noise model is indeed usually assumed to be iid
Gaussian.

>	used?  If the noise process is calculated
>	using the calculated ARMA parameters and
>	the original sampled data, doesn't that mean
>	that your model will be identical to the
>	sampled data?

Yes.  Your model will (ideally) match the spectral shape of the process you
observed.  That is is the meaning and purpose of the whole exercise, to
estimate the spectral content of the signal.

>
>Thanks for your thoughts!

Such as they are at 3:41 am :-)
If you want some background literature, I'm personally rather fond of 

@book{Mar87,
      author = "Marple, S.~Lawrence,~Jr.",
      title  = "Digital Spectral Analysis with applications",
      publisher = "Prentice-Hall",
      year   = 1987
}

It discusses everything mentioned here (and quite a bit more) and comes with a
disk full of FORTRAN routines that implement all the common algorithms.

SR

jbuck@epimass.EPI.COM (Joe Buck) (11/11/89)

In article <1389@mrsvr.UUCP> kohli@gemed.ge.com (Mr. Bad Judgment) writes:
>1. Are the AR "parameters" the poles of the waveform?
>	If so, does this mean that an AR parameter of magnitude
>	> 1 indicates an unstable system (assuming no zero to
>	counteract it)?

There are many different families of AR parameters.  One could
consider the poles the parameters; then the requirement for stability
is that they be inside the unit circle.  But the most common choices
of parameters to use are either the filter coefficients, the
coefficients of the A(z) polynomial in the model, where
           1
   H(z) = ----
          A(z)

so the poles are the roots of A(z).  But you're probably thinking about
the reflection coefficient parameters, which can be obtained from the
autocorrelations or the filter coefficients using the Levinson recursion.
The requirement for stability is that all reflection coefficients have
magnitude <= 1.

>2. (And the MA parameters *are* the zeroes?)

It's sufficient to specify the poles and zeros, though people usually
specify the coefficients of the numerator and denominator polynomial
instead.

>3. Is the noise process usually neglected when
>	using the AR and MA parameters to model the
>	time series?

Not at all.  It's usually explicitly assumed that the process is the
result of a white noise process being filtered by a pole-zero filter.
In some methods, we basically look for a filter that will convert the
data into white noise; the model is simply the inverse of this filter.

-- 
-- Joe Buck, just visiting/consulting at Entropic
-- write me at: jbuck@janus.berkeley.edu 
		...!{uunet,ucbvax}!janus.berkeley.edu!jbuck

smit@.ucalgary.ca (Theo Smit) (11/15/89)

In article <1389@mrsvr.UUCP> kohli@gemed.ge.com (Mr. Bad Judgment) writes:
>
>I've been working a little bit with some ARMA stuff
>and have run into some questions which I can't
>generalize an answer to from my experience:
>
>1. Are the AR "parameters" the poles of the waveform?
>	If so, does this mean that an AR parameter of magnitude
>	> 1 indicates an unstable system (assuming no zero to
>	counteract it)?

Yes. Common ways to 'create' a stable system include shrinking all of the
pole radii by some factor (alpha). Then the AR polynomial becomes:

                             -1        2     -2        3     -3
	A(z) = 1 + alpha a1 z   + alpha  a2 z   + alpha  a3 z   + ...

alpha will be slightly less than 1.0. (0.99, or 0.999; it doesn't take much).
Mostly this is used to counter arithmetic errors; if your system is unstable
the poles will likely be well outside the unit circle and you have no hope in
heck of doing much about it.

>2. (And the MA parameters *are* the zeroes?)

Right again. In general you can set up an ARMA process as separable AR and
MA processes, i.e.

	        B(z)
	H(z) = ------
	        A(z)

where B(z) is the MA process, and A(z) is the AR process.

>3. Is the noise process usually neglected when
>	using the AR and MA parameters to model the
>	time series?  If not, is a Gaussian model
>	used?  If the noise process is calculated
>	using the calculated ARMA parameters and
>	the original sampled data, doesn't that mean
>	that your model will be identical to the
>	sampled data?

As long as there is a noise process, the process described by the model cannot
be identical to the sample process. At the 'ideal' model order (whatever that
is), the noise process will (should, anyway) be perfectly white, i.e. totally
uncorrelated, since the model has predicted all there is to predict.
>
>Thanks for your thoughts!

No trouble!

>
>Jim Kohli
>GE Medical Systems
>ge.crd.com!gemed!hal!kohli

Theo Smit (smit@enel.ucalgary.ca)

rob@kaa.eng.ohio-state.edu (Rob Carriere) (11/19/89)

In article <2074@cs-spool.calgary.UUCP> smit@enel.ucalgary.ca (Theo Smit)
writes: 
>In article <1389@mrsvr.UUCP> kohli@gemed.ge.com (Mr. Bad Judgment) writes:
>>1. Are the AR "parameters" the poles of the waveform?
>>	If so, does this mean that an AR parameter of magnitude
>>	> 1 indicates an unstable system (assuming no zero to
>>	counteract it)?
>
>Yes. 

Wrong.  The AR parameters are the coefficients of a polynomial whose _roots_
are the poles.

[stability force deleted]
>Mostly this is used to counter arithmetic errors; if your system is unstable
>the poles will likely be well outside the unit circle and you have no hope in
>heck of doing much about it.

If your system is unstable and you force the estimate to be stable you have a
meaningless estimate.

>>2. (And the MA parameters *are* the zeroes?)
>
>Right again. [...]

Wrong again.  The zeroes are the _roots_ of the polynomial described by the
MA parameters.

SR

kohli@gemed (Mr. Bad Judgment) (11/20/89)

rob@kaa.eng.ohio-state.edu (Rob Carriere) writes:
<In article <2074@cs-spool.calgary.UUCP> smit@enel.ucalgary.ca (Theo Smit)
<writes: 
<>In article <1389@mrsvr.UUCP> kohli@gemed.ge.com (Mr. Bad Judgment) writes:
<>>1. Are the AR "parameters" the poles of the waveform?
<>>	If so, does this mean that an AR parameter of magnitude
<>>	> 1 indicates an unstable system (assuming no zero to
<>>	counteract it)?
<
<[AR parameters are the coefficients of a polynomial whose _roots_
<are the poles]
<
<[stability force deleted]
<>Mostly this is used to counter arithmetic errors; if your system is unstable
<>the poles will likely be well outside the unit circle and you have no hope in
<>heck of doing much about it.
<
<If your system is unstable and you force the estimate to be stable you have a
<meaningless estimate.
<
	That is true.  On the other hand, the system I am
modeling is stable, consisting exclusively of damped sinusoids
of various frequencies and additive noise (so why not use
Prony's method?  I'm doing that, too).  Marple's book suggests
(and rightfully so) that if the results of doing AR analysis
generates poles which have magnitude greater than one, attempts
to model the system will be disappointing because of the
instability in the model.  This may be easily confirmed upon
inspection of an AR(1) system, for example, where

	X(t) = A(t)X(t-1) + n(t).

	X(t) may be expected to grow without bound if |A(t)| >
1.0.  If the AR coefficients are the coefficients of
polynomials whose roots are the poles are the system (we know
they are), what are the consequences of scaling the
coefficients?  It seems to me that the frequency components
should be preserved, and if |variance of n(t)| were also scaled
down, the modeled X(t) sequence would be stable, but scaled
incorrectly.  Is this wrong?  I would like to have stable
estimates for this known-stable system, but most of the
recommended AR techniques in Marple generate unstable
parameters (Marple makes note of this possibility in a couple
of places, but I didn't see any recommendations about how to
fix the problem.  I *think* there's an implication that the
"next better" algorithm should be used).

Opinions?  Facts?
Either are welcome, just please be kind.

Thanks again,
Jim Kohli
GE Medical Systems

P.S.  Thank you Theo, Rob, and Stephen, the discussion has been
interesting and enlightening (so far!).

rob@kaa.eng.ohio-state.edu (Rob Carriere) (11/22/89)

In article <1456@mrsvr.UUCP> kohli@gemed.ge.com (Mr. Bad Judgment) writes:
>rob@kaa.eng.ohio-state.edu (Rob Carriere) writes:
><In article <2074@cs-spool.calgary.UUCP> smit@enel.ucalgary.ca (Theo Smit)
><writes: 
><[stability force deleted]
><>Mostly this is used to counter arithmetic errors; if your system is unstable
><>the poles will likely be well outside the unit circle and you have no hope
   in 
><>heck of doing much about it.
><
><If your system is unstable and you force the estimate to be stable you have a
><meaningless estimate.
><
>	That is true.  On the other hand, the system I am
>modeling is stable, consisting exclusively of damped sinusoids
>of various frequencies and additive noise (so why not use

Then a stability force can make sense.

>Prony's method?  I'm doing that, too).  Marple's book suggests

Sounds like a good idea.  Especially comparing results can be helpful.

>(and rightfully so) that if the results of doing AR analysis
>generates poles which have magnitude greater than one, attempts
>to model the system will be disappointing because of the
>instability in the model.  This may be easily confirmed upon

Not only that, but many algorithms that play with these models are numerically
unstable if the model is not stable.  So you get hosed from two sides :-)

>[...]  If the AR coefficients are the coefficients of
>polynomials whose roots are the poles are the system (we know
>they are), what are the consequences of scaling the
>coefficients?  

Scaling the coefficients (multiply the k-th coeff by alpha^k for some
0<alpha<1) has the effect of moving all the poles in by alpha (ie:

    i phi                   i phi
 r e        ---->  alpha r e

So geometric relations between the poles are preserved, but everything is
fitted into a smaller circle.

Another possibility is to calculate all the poles and move the unstable ones
just enough to stabilize them.  This sounds nice, but I have slightly horrible
experiences with it (and would be interested to hear from people who have good
ones).

Then there's calculating the partial covariances (or partial autocorrelations,
or whatever you want to call them today) and setting the magnitudes there.
Never tried this one.  Any takers?

>It seems to me that the frequency components
>should be preserved, and if |variance of n(t)| were also scaled
>down, the modeled X(t) sequence would be stable, but scaled
>incorrectly.  Is this wrong?  I would like to have stable
>estimates for this known-stable system, but most of the
>recommended AR techniques in Marple generate unstable
>parameters (Marple makes note of this possibility in a couple
>of places, but I didn't see any recommendations about how to
>fix the problem.  I *think* there's an implication that the
>"next better" algorithm should be used).

The general attitude seems to be that "as along as you futz only a little, it
doesn't matter".  Of course that doesn't help you much if you have to futz by
a lot...  I do not know of any good, general, rigorous treatment of this (very
real) problem (again, anybody out there know better?  Speak up please!)

I would say that all these ad-hoc methods are highly suspect if you have to
move any pole by more than a little.  At least, my experience has been that in
such cases you might as well throw away the estimate and try again.

There are methods that are guaranteed to give you stable estimates (they are
in Marple, I'm pretty sure), but the price you pay is bias in the estimates.
As your data-length increases, the bias goes down, but then, so does the
chance of getting an unstable estimate out of the other methods.  In short,
getting a good (or even reasonable) estimate out of short data records is
_hard_. 

If you are indeed using short records, you may want to try modifying your
algorithm to incorporate a SVD.  In both ARMA and Prony algorithms, you are
solving an equation of the form:

Ya = -z

Create the matrix [ z Y ] and do an SVD on it.  Say you get USV.  Now if you
are expecting the "correct" order to be n, keep n singular values.  Call this
partially zero-ed S matrix S'.  Compute a to be the LS solution of

Y'a = -z'

where 

[ z' Y' ] = US'V

The way you do this is start out with forming your Y and z for a much higher
order (even up to Number-of-data-points/2) and then use the SVD to go down to
order you want.  Of course your a polynomial will still be high order, ie have
a whole slew of spurious roots.  You need to get rid of those.  If the system
is known stable you can just throw away al unstable poles.  For backward
prediction that should take care of the spurious guys.

SR

smit@.ucalgary.ca (Theo Smit) (11/28/89)

In article <3556@quanta.eng.ohio-state.edu> rob@kaa.eng.ohio-state.edu (Rob Carriere) writes:
>In article <2074@cs-spool.calgary.UUCP> smit@enel.ucalgary.ca (Theo Smit)
>>Yes. 
>
>Wrong.  The AR parameters are the coefficients of a polynomial whose _roots_
>are the poles.
[...]
>>Right again. [...]
>
>Wrong again.  The zeroes are the _roots_ of the polynomial described by the
>MA parameters.
> SR

OK, so I goofed. Anyway, what Rob says is what I _meant_. Most of the time
the AR and MA polynomials are more useful than the explicit pole locations
(i.e. doing an FFT on the coeffs to get the periodogram of the model), so
I usually don't bother finding the actual poles.

On a different note, doing FFT's on the AR or MA parameters is computationally
much less intensive if you use the interpolating FFT I described earlier.
(that would be the Screenivas-Rao algorithm, not mine; mine works well for
interpolating time sequences when you start out in the frequency domain).
For all of you who responded to my previous posting on the EZFFT, I will be
posting the source and some more explanation shortly. Stay tuned.

Theo Smit (smit@enel.ucalgary.ca)

smit@.ucalgary.ca (Theo Smit) (11/28/89)

In article <1456@mrsvr.UUCP> kohli@gemed.ge.com (Mr. Bad Judgment) writes:
[...previous discussion deleted]
>	That is true.  On the other hand, the system I am
>modeling is stable, consisting exclusively of damped sinusoids
>of various frequencies and additive noise (so why not use
>Prony's method?  I'm doing that, too).  Marple's book suggests
[]
>	X(t) = A(t)X(t-1) + n(t).
>
>	X(t) may be expected to grow without bound if |A(t)| >
>1.0.  If the AR coefficients are the coefficients of
>polynomials whose roots are the poles are the system (we know
>they are), what are the consequences of scaling the
>coefficients?  It seems to me that the frequency components
>should be preserved, and if |variance of n(t)| were also scaled
>down, the modeled X(t) sequence would be stable, but scaled
>incorrectly.  Is this wrong?  I would like to have stable
>Opinions?  Facts?
>Thanks again,
>Jim Kohli
>GE Medical Systems
>P.S.  Thank you Theo, Rob, and Stephen, the discussion has been
>interesting and enlightening (so far!).

I assume by damped sinusoids you mean sinusoids with exponentially decaying
amplitude? (Are you by any chance doing MRI data? We've done a lot of
experimentation on modeling MRI)
Let's assume that we have an exponentially decaying unit step. Then the time
series of the output is:
			        n
			x(n) = a  , a < 1.0, 0 <= n < inf.

The AR(1) model is:	       -1
			1 + a z

and if a < 1.0, the thing is stable. However, if we look at the thing 
backwards, we see an exponentially increasing series. If we model this
sequence, we get an AR coefficient of 1/a*, that is, outside the unit circle.
The coefficient is also conjugated from the forward model coefficient (I
think. I'm doing this from memory), so the overall angle of the pole in the
Z-plane stays the same, it's just been reflected outside the unit circle.
Obviously, this situation is unstable.
If your algorithm applies forward _and_ backward prediction to decaying data,
your results will not be correct, since the model cannot provide a stable
solution in both directions. The coefficients will be the result of a 
best-fit in both directions, and will not mean much of anything. This
is a problem I found in using the Burg algorithm; anything that's not
stationary doesn't get modeled worth &*@Q$.

About the pole radius scaling - if you don't care what the absolute amplitude
is (ie in spectral analysis, you're usually only interested in relative
magnitude), there is no problem. The pole scaling results in reduced
resolution in the frequency domain, similar to what windowing the data
does. If you're trying to model the data exactly enough to be able to
recreate it, a different approach will be necessary.

Hope this helps,

Theo Smit (smit@enel.ucalgary.ca)