Second Order Processes : AR and ARMA Representations

Second Order Processes : AR and ARMA Representations

The Autocorrelation Function

Let be a continuous-timewss process with autocorrelation

. (1)

The reason for using the above notation is that (1) can be associated with a second order random process modeled by white noise passing through a linear system with a transfer function of the form

.(2)

The system (2) is an underdamped second order system with damping ratio ξ<1, undamped natural frequency , and damped natural frequency . The general form of the numerator of (2) admits the broadest class of such processes. The purpose of this work is to investigate the exact relationship between second order processes that can be represented by fictitious white noise passing through (2), and their associated autocorrelation functions. We will begin by investigating the relationship between (2) and (1), since, as noted above, (1) is, in a sense, the paradigm model for the autocorrelation function associated with second order processes.

The Transfer Function Associated with

To arrive at the specific expression for (2) we will utilize the following Laplace transform pair relation:

.(3)

From this relation, it follows that the psd is

.(3)

This can be simplified to

.(4)

In turn, (4) can be expressed as

where

.(5)

From (5), we see that the autocorrelation function (1) corresponds to a second order process that can be represented as unit variance white noise passing through a transfer function of an underdamped second order system that has a zero at . We now proceed to determine the form of the correlation function for a process that can be modeled by white noise passing through a second order underdamped system whose transfer function has no zeros.

The Autocorrelation Function Associated with

The psd of a second order process, represented by white noise with variance passing through the system

(6)

is

(7)

The transfer function (6) is, in a sense, the classic form of the second order underdamped system. It has no zeros. To recover express (7) as

.(8)

By cross-multiplying the terms in (8) to obtain a single fraction, and equating its numerator to the term

in (7), we obtain , which gives . Hence,

.(9)

And so, for , we have

.(10)

Hence, we see that a second order process modeled by white noise passing through the classic underdamped second order system has an autocorrelation function that is distinctly different from the classic autocorrelation function (1).

A Detailed Comparison of the Two Autocorrelation Functions

In this section we will investigate the details of the differences between (1) and (10). To begin, from (10) we have

.(11)

We will assume that the processes and have one and the same variance, . In this case, (1) is

(12)

and (10) becomes

.(13)

From (12) and (13) we have

.(14)

Hence, for small values of the damping ratio , the relative difference between (12) and (13) will be negligible. This is illustrated in Figure 1.

Figure 1. Plots of (12) (BLUE) and (13) (RED) for different damping ratios, and .

From this figure we note that (12) and (13) are nearly identical for smaller damping ratios. This lead one to ask whether there should be any real concern for whether a second order process with the classic autocorrelation function (1) has an all-pole or a pole-zero filter representation. Similarly, in relation to a second order process with the classic all-pole filter representation (6), in the case of smaller damping ratios one can reasonably assume that it has the classic autocorrelation function (1).

A Detailed Comparison of the Two Power Spectral Density Functions

In this section we will investigate the details of the differences between the power spectral density functions (psd’s) (4) and (7). We repeat them here for convenience. For X(t) we have

.(15)

In relation to Y(t), from (11) we have , and so (7) becomes

.(16)

Clearly, there is a significant difference between (15) and (16) at high frequencies, since (15) rolls off at a rate of 20 dB/decade, while (16) rolls of at a rate of 40 dB/decade. Notice that the dc gain of (15) is 3 dB lower than (16). This is to be expected, since (15) has a smaller high frequency attenuation, and the processes have equal total power. Perhaps most important is the fact that at , (15) and (16) are equal. This, combined with the difference in dc gain, leads to the observation that if one were to use the relative (to dc) amplification at resonance (i.e. the system Q-factor) to estimate the damping ratio in a second order model for X(t), it would be over-estimated by 3 dB. Plots of (15) and (16) in relation to the autocorrelations given in Figure 1 are given in Figure 2.

Figure 1. Plots of (15) (BLUE) and (16) (GREEN) for different damping ratios, and .

The Sampled Second Order Processes and

In this section we develop the difference equations corresponding to the discrete-time processes and modeled by white noise passing through the sampled system transfer functions and , corresponding to the continuous-time transfer functions given by (5) and (6), respectively. To this end, we will use a table of Laplace to z transforms obtained by sampling the system impulse response. The discrete-time system transfer functions are found to be

(15)

and

.(16)

The difference equation associated with (15) is

(17a)

where

;

and

.(17b)

The difference equation associated with (16) is

(18a)

where

.(18b)

Notice that (17) is an ARMA(2,1) process, whereas (18) is (effectively) an AR(2) process.

Numerical Example 1: and

Numerical Example 1: and

APPENDIX

After scouring Google, the closest Fourier transform pair I could find was the following:

Since is an even function of τ, from the book relation (A.3) on p. 464 we have:

. The term in the brackets can be simplified to a real-valued function of ω as follows:

.

Combining these two terms into a single ratio gives:

.

The denominator is:

Hence,

. And so

Write . Then where

The Impulse Response:

. Hence,

, or

Write

Figure 1a. Impulse responses for 3 values of ξ.

For the Matlab command impulse(sys) is shown below.

Figure 1b. Impulse responses for 3 values of ξ, computed vie Matlab impulse(sys).

This figure is exactly the same as Figure 1a, and so verifies the derived expression for h(t).

Summary: For :

=

and

We now consider a system with transfer function

.(1)

and with impulse response

(2)

Then

(3)

Our goal here is to recover . To this end, we will ‘speculate’ that we can write:

.(4)

By cross-multiplying the terms in (4) to obtain a single fraction, and equating its numerator to the term

in (3), we obtain , which gives . Hence,

.(5)

And so, for , we have

.(6)

Plots of (6) for different damping ratios are shown below.

Figure 2a. Plots of (6) for different damping ratios .

From the third plot it is not clear that we have a valid autocorrelation function. That it is valid is confirmed by the zoomed region below:

Figure 2b. Plot of the first portion of the third plot in Figure 2a.

Conclusion: The ‘standard second order underdamped system (without a zero) has an autocorrelation function that is not a decaying cosine. That form of autocorrelation belongs to a system with a zero.

Including the Factor of (1/T) in Going from the Laplace to Z-Transform

Consider the Laplace to z-transform pairs below.

Commentary on the Above Pairs

Example 1.

Consider the impulse response with transfer function . For a sampling interval, T, where . In this case, we have

where .

Hence, the above transforms 1. and 2. are not entirely correct, as the z-transforms are missing a factor of T.

Now, suppose that we have a random process X(t) with autocorrelation function . We can express this as:

where is the Kronecker delta function. It is not the Dirac delta function. The psd is

.

.

Letting , we then have

Now, since , The sampled process has psd:

where

is the z-transform of ; which is what the table that gave 1. and 2. above would give. Even though continuous-time white noise does not exist, the above psd for the sampled process can be viewed a discrete-time white noise (which does exist) passing through the filter . This fictitious white noise, call it U(kT), has variance .

Conclusion 1: It is essential that the factor (1/T) be included when going from Laplace to z-transforms.

Using Matlab’s fft and ifft Codes in Relation to Autocorrelations & PSDs

For a continuous-time process X(t), the Fourier transform pair relation between and is:

The sampled process has the frequency range where the sampling frequency is . Hence, the sampled process transform pair is:

Now, suppose that we wish to evaluate at N uniformly spaced frequencies over the interval . Then where . In this case, we have

.

Note that , and so , and . It follows that . Define the normalized frequency . Then if we assume that for , we have

.

The Matlab fft.m and ifft.m commands refer to the following pairs:

.

Hence, once again we see that the factor of T is missing!

Example 1 continued.

Now, suppose that we have computed the psd via . Then to recover we cannot simply use the ifft.m command applied to . Instead, we must first multiply by 1/T.

CONCLUSION: In going from a continuous-time to a discrete-time random process modeled as white noise passing through a filter, the sampling interval T must be accounted for in relation to (i) computing the white noise variance associated with the normalized filter obtained from a Laplace to z-transform table that does not include that factor, and (ii) in using the ifft.m command to recover the the autocorrelation function from a given expression for the psd.

The following continuous-time impulse response to be sampled are:

, which can be expressed as

and.

, or

<- WRONG!

.

Numerical Example. Set . Then .

The continuous-time transfer functions are:

and

and .

#20:

Matlab Code in folder aeem/ee573/s2010/exams

%PROGRAM NAME: arma21.m

%This code address a wss process with R(t)=s^2 *exp(-zw_n|t|)cos(w_d t)

%------

% REUQIRED INPUT: z = damping ratio & sampling ratio ws/wn

z=0.707; rsamp=100;wn=1;

ws=wn*rsamp;

T=2*pi/ws;

%======

wd = wn*(1-z^2).^0.5;

zwn = z*wn;

c=exp(-zwn);

d=z/((1-z^2)^0.5);

tau = zwn^-1; % max times of plots of impulse responses

dt = .01*tau; % Sampling intervals

t=0:T:1000*T; %plot times

%======

%Autocorrelation Functions

th = atan(d^-1);

RX = exp(-zwn*t).*cos(wd*t);

RY = RX + d*exp(-zwn*t).*sin(wd*t);

figure(1)

plot(t,RX,t,RY,'r')

grid

title('R(t)')

pause

%======

% PSDs via Filter Expressions

wmin = .01*wn;

wmax = wn*rsamp/2;

nw = 2000; % number of frequency points

dw = (wmax-wmin)/nw;

w = wmin:dw:wmax-dw;

%------

nH=(2*zwn)^.5 *[1 wn];

dH=[1 2*zwn(1) wn^2];

H=tf(nH,dH);

nG=(4*zwn*wn^2)^0.5;

G=tf(nG,dH);

SX(1,:)=bode(H,w);SX=SX.^2; SXdB=10*log10(SX);

SY(1,:)=bode(G,w);SY=SY.^2; SYdB=10*log10(SY);

figure(2)

semilogx(w,SXdB,w,SYdB,'r')

title('Continuous-Time PSDs')

grid

pause

%======

% Construction of ARMA & AR Models

%======

%AR coefficients

z22=(1-z^2)^0.5;

wdT=wd*T;

zwnT=zwn*T;

a1=2*exp(-zwnT).*cos(wdT);

a2=-exp(-2*zwnT);

% MA coefficient

b1= (((1-z)/(1+z))^.5*sin(wdT)-cos(wdT))*exp(-zwn*T);

sigux = T*(2*zwn)^.5;

siguy=2*T*exp(-zwnT)*sin(wdT)*zwn^.5 / z22;

% Compute PSDs

nfft=10000;

a = [1 ; -a1 ; -a2];

A=fft(a,nfft);

bx = [1; b1];

BX=fft(bx,nfft);

by=[0 ; 1];

BY=fft(by,nfft);

SSx=(sigux*abs(BX./A)).^2; SSX=SSx(1:nfft/2);SSXdB=10*log10(SSX);

SSy=(siguy*abs(BY./A)).^2; SSY=SSy(1:nfft/2);SSYdB=10*log10(SSY);

wvec=1:nfft; wvec=(ws/nfft)*wvec; wvec=wvec(1:nfft/2);

figure(3)

semilogx(wvec,SSXdB,wvec,SSYdB,'r')

grid

title('ARMA (blue) and AR (red) PSDs')

xlabel('Frequency')

ylabel('dB')

pause

% Comparison of Continuous & Discrete Process PSDs

figure(4)

semilogx(w,SXdB,wvec,SSXdB,'r')

xlabel('Frequency (rad/sec)')

ylabel('dB')

title('Comparison of PSDs for X(t) and X(kT)')

grid

pause

figure(5)

semilogx(w,SYdB,wvec,SSYdB,'r')

xlabel('Frequency (rad/sec)')

ylabel('dB')

title('Comparison of PSDs for Y(t) and Y(kT)')

grid

pause

%Comparison of Cont and Discrete Process Acorrs

figure(6)

plot(RX)

nn=length(RX);

hold on

RXX=(1/T)*real(ifft(SSx)); RXX=RXX(1:nn);

plot(RXX,'r')

xlabel('Lag')

title('Comparison of Autocorrelations for X(t) and X(kT)')

pause

figure(7)

plot(RY)

hold on

RYY=(1/T)*real(ifft(SSy)); RYY=RYY(1:nn);

plot(RYY,'r')

xlabel('Lag')

title('Comparison of Autocorrelations for Y(t) and Y(kT)')

pause

%Generate realization of both processes

u=(1/T)*randn(1,1500);

x=zeros(1,1500); y=zeros(1,1500);

ux=sigux*u;

uy=siguy * u;

for k=3:1500

x(k)=a1*x(k-1)+a2*x(k-2)+ux(k)+b1*ux(k-1);

y(k)=a1*y(k-1)+a2*y(k-2)+uy(k-1);

end

figure(8)

plot(x(501:1500));

hold

plot(y(501:1500),'r')

xlabel('Time')

title('Samples of X(kT) (BLUE) and Y(kT) (RED)')