Impulse Response Determination in the Time Domain ... - IEEE Xplore

1 downloads 0 Views 612KB Size Report
Impulse Response Determination in the Time Domain-Theory. Abstract-Two methods are presented for determining the impulse response of an object in the ...
657

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. AP-30, NO. 4, JULY 1 9 8 2

Impulse Response Determination in the Time Domain-Theory

Abstract-Two methods are presented for determining the impulse response of an object in the time domain when both the input and output time domain waveforms are specified. Pole extraction features can then be applied to the nonimpulsive portion of the impulse responseto determine the singularity expansion method (SEM) parameters of the object. The first method involves synthetic division and is comparatively straightforward.The second technique is a least squares method which is computationally more stable. The effect of measurement errors in the input and output waveforms is evaluated for each method.An investigation is made as to what formof the input minimizesthenoise variances in the computation. Finally a generalized least squarestechniqne is presented, which yields a minimum variance unbiased estimate for the impulse response when the noise covariance matrix is known.

I. BACKGROUND

I

N ORDER to characterize an unknown object by remote sensing, certainelectricalproperties of the object may be used [ I ] -[9]. In particular the impulseresponse, ramp response, and the natural resonantfrequencies of the object have been proposed. In order to probe the object accurately, one mustilluminate it with a band of frequencies whose wavelengths are approximately the same dimensions as the overall length of the object. Since, in general, the size of the object is notknownbeforehand,it is necessary to illuminate the object by a broad-bandsignal such as a waveform which is an approximation to an impulse.The transmitted pulse induces electrical currents on the object and reradiates the incidental energy. The problem is then to obtain theelectrical properties of the object from the transmitted and received time domain waveforms. This area has been investigated by researchers at The Ohio StateUniversity [ 11-151. In recent years with the advent of the singuiarity expansion method (SEM), efforts have been directed toward characterizing an object by its natural frequencies-poles and zeros [ 6 ]18 1 . For this purpose Prony’s method has been applied with much success where analyzing measured impulse responses with high signal-to-noise ratio. However it is difficult to extend Prony’s methodtoarbitraryinput and arbitrary output waveforms. Also the SEM representationdoesnot account for theimpulsive portion of the system time behavior. Hence the complete impulse response contains more information than the SEM poles. One of the objectives of this paper is to obtain the impulseresponse in the time domainwhen arbitrary inputand output waveforms are given.

More recently the pencil of function method [ 101-[ 141 has been applied withmuch success to obtainthenatural measured arbitraryinput and frequencies of objectsfrom outputtime waveforms. Both Prony’s and the pencil of functionmethods are frequencydomain techniques. The accuracies of determining impulse responses byeither of these techniques dependdirectly on thetimedomaindata available. Also asignificantproblemwith anysuch parametric scheme is the need to estimate the number of dominant poles and zeros. A poor estimate can lead to large errors in the results. Finally as discussed in the next section, it is not easy to obtain the impulseresponse from given input and output timedomain waveforms using fastFourier transform techniques [ 91 , The objective of this paper is to present two techniques for obtainkg the impulseresponse of an object directly in the time domain from measured input and output waveforms. 11. INTRODUCTION

Let x ( t ) be theinput to a time invariant causal linear system which is characterized by its impulseresponse h ( t ) . The corresponding output y ( t ) is given by

One approachfor solving (1)forthe impulseresponse involves the Laplace transform. The result is

Y(s) = H(s) ’ X(s) where X ( s ) , Y(s), and H(s) are the Laplace transforms of and h ( t ) , respectively. The impulseresponse h ( t ) is obtained by taking the inverse Laplace transform of (2) to yield

x(t), y(t),

where I‘ is the Bromwich contour. Numerical problems are frequentlyencounteredwith this approach when dealing with real data [ 91. For example, when using fast Fourier transformtechniques,the division blows up if X G w ) = 0 in the pass band. Manuscript received September29,1980; revised November 9, 1981. Twotimedomainmethodsare presented in this paper. This work was supported by the Office of Naval Research under ConOne is themethod of synthetic division while the other is tract N00014-79-C-0598. T. K. Sarkar and S. A. Dianat are with the Department of Electrical the method of leastsquares. The former is computationally Engineering, Rochester Institute of Technology, One Lomb Memorial straightforward whereas thelatteriscomputationallymore Drive, Rochester, NY 14623. stable. D. D. Weirwr is with the Department of Electrical Engineering, SyraIn this presentation y ( t ) and x ( t ) are assumed to be discuse University, Syracuse,NY 12310. V. K. Jain is with the Department of Electrical Engineering, Univer- crete time signals with identical sampling rates. The problem sity of South Florida, Tampa,FL 33620. i s to process the discrete datapoints so as toestimatethe 001 8-926X/S2/0700-0657$00.75 0 1982 IEEE

.

.-

x

.

-

..

~.

658

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION,VOL. AP-30, NO. 4,JULY 1982

correspondingdiscreteimpulseresponse without introducing additionalinformation like the augmented high frequency response [ 9 1 . Let the ith samples of y ( t ) , x ( t ) , and h ( t ) be denoted by y i , x i , and hi, respectively. Also assume that y ( t ) and x ( t ) are available as the f ~ t time e series } = bo,y l , y 2 , YN }and { x } = { x g , XI, x 2 , X M 1 The object is to estimate the finite impulse response time series {h}= {ho, hl , h2, --,

b

b}

-*e,

*e*,

hL

Example 1 : As a simple example, consider the rectangular pulse { x } = { 1, 1, l} to be an input to a system whose output is the triangular pulse = (1, 2, 3, 2, 1). Using synthetic division, as implied by (4), theunknowntransferfunction (h} is given by k(z) =

1

+ 2z-l + 3z-2 + 2 z - 3 + z 1+2-1+z-2

1.

For convenience, the following assumptions are made.

- ~

= 1+ z - ' + z - 2 .

a) hi = 0 f o r i > L . b) x i is observed for 0 < j < M and is not identically zero. C ) X j = OforO>jandj>M. is observed for 0 < k < N (=M L ) . d) y k

+

(9)

This implies {h} = { 1, 1, 1). Alternatively if it is assumed Ar = 1, application of (8) resultsin {h} = (1, 1, l}. This agrees with the solution obtained by synthetic division. IV. EFFECT OF MEASUREMENT ERRORS

111. BA HLI'S TECHNIQUE (SYNTHETIC DIVISION) 1151, [I61 This method is applicable for all positive integer choices of N , M, and L . The sequence corresponding t o h ( t ) is written symbolically as

where the division of sequence { y } by sequence {.}is carried out by synthetic division. Justification of the synthetic division operation is obtained from the theory of z-transforms. The input and output z-transforms can be written as

The effect of measurementerrors on the synthetic division scheme is nowconsidered. Assume thattheinput measurement' is corruptedby additivenoise p ( t ) so that the total recorded input is x p . Also assume that the output measurement is corrupted by noise q ( t ) such that the total measured output is y q. Using samples of the corrupted input and output meaprements, an estimate of the impulseresponse, denoted by h(z), is given by

+

+

In this section the relationship of the estimate @}to the {h} is exploredwhere it is assumed true impulseresponse that the measurement errors are small and that the elements of { q } and (i.} which pertain to two different measurements are statistically independent zero mean random variables. The expression in (IO) can be rewritten as s(z) 1+-

-

h ( z ) is then the ratio of the two polynomials given by (5) and ( 6 ) . Since two z-polynomials are divided inthesame manner as algebraic polynomials, the operation of synthetic division is justified. Alternatively a discrete approximationto (1) is i

hY(Z) ( z ) = Y(Z) x@>

(7)

P(Z)

1+-

x(z)

Forthose values of z for which I p ( z ) / x ( z ) I < 1, [ 1 4@ ( z ) / x ( z ) ) ]- l may be expanded in a binomial series to yield

p o t p1z-

I= 0

~.

+ -.*+ PMZ

I-'

[;+xo+xlz-'+~-~+xMz-~

where Ar is the sampling interval. A computationally compact form forhi is then given by

i y i - Ar h. = I

It follows for these values of z that

xihi- i i= 1

x.

Ar

(8)

where it is assumed x . # 0. The values for hi obtained in this way differ from those obtained by synthetic division by the constant factor l/Ar. Note that both of these methodsare applicable to discontinuous inputs and outputs.

X(Z)

=-

4

659

DETERMINATION SARKAR RESPONSE e t al. : IMPULSE

The expected value of a random sequence is defined to be the sequence whose elements are expectatiop of the original elements. Since each element represents a sample value of the waveform, the expectation of the element .$orresponds to an ensemble average of the waveform at the instantof the sample value. Then the expectedvalue of the estimate E(.) is given by

is s t i l l associated with the output.However if p i f 0 then there is a bias in the solutionwhich is given as

UP0

=h(z) (x0

+

UP1

2 -2

+ ...

+ xlz- + + X M z - M ) 2

By definition, the variance of a sequence is the sequence whose elements are thevariance of the original elements. Each variance is the variance of the original random waveform at the corresponding sampling instant. Note that

where

and

E[q(z)]=E[qo

+ 412-1 + q2z-2 + -*I As a first approximation

+-.

=E[qo] +E[q,]z-l +E[42]Z-2

=o

(15)

where use hasbeen made of the assumptions of zero mean and statistical independence of the random variables p i and 4i. Similarly mean sequence Notethat ( 2 1 ) corresponds tothezero (co, c1, c 2 , -}. It follows that the variance of this sequence, results in which also equals the variance of

x@),

and

= UP0

2

+ op12z-2 + u

~

~

+ E [ c 1 2 ] z - 1 + q c 2 2 3 z - 2 + ..= uo2 + u12z-1 + (T22z-2 + (22)

02[Z(Z)] = E [ c o 2 ]

+ 2pop1z-I + 2P0P2.-2 + p 1 2 z - 2 + 2 ~ , p ~ z - ..~+ + p 2 ~ - 4+ 2p2p3z-’ + ...I

+ --‘

E[p(z)] =E[$

-.*

+~

z

- (17) ~

where ( ~ i ’ = E [ c : ] . Use of ( 2 2 ) and(16),inconjunctionwith inequality [ 171 , yields the inequality

Chebyshev’s

1 + kUi] 2-1 - 7 k

(23)

P[hi’ - k U i < i;i < h i

where

up; = E [ & 2 1 .

&

Consequently

E [ K(z)] = h(z) = ho’

where P [ ] denotes the probability of the argument in the brackets and is the element in {E}. Equation ( 2 3 ) is valid for small measurement errors. In summary, if the first- and second-order statistic of the random sequences b ) a n d {q)are known, then ( 2 2 ) and (16) yield approximate expressions for the variance and mean of the sample values of the estimate of the time domain impulse reponse obtained from (1 0). Example 2: Consider a system with impulse response { h ) = { 1, l}, for which the input is (x} = {1,1), and the output is b}= { 1, 2, 1). Assume, for the purpose of illustration, that the input and output are each contaminated by zero mean additive measurementerrorswhich are uniformly distributed between the values +0.0176 and -0.0176.Notethatthe measurement errors have zeromeansand variances u2 =

[

2

1+

+ h l’z-

OPO

2z-2 +

UP1

+ x 1z- + + +h2’Z4 +

(x0

**.

*.-

+ ...

X M Z - y

1 (1 9)

where h ( z ) = y ( z ) / x ( z ) . Thus if the varia&es of theinput measurement errors are known,approxima values forthe expectation of the estimate of the time domainimpulse response $(z) can beevaluated. Notice that.if pi 0 for all i (i.e., the input measurement is free from any additive noise), then the estimateis unbiased even though a measurement error

4f

~-

660

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. AP-30, NO. 4, JULY 1982

l o u 4 . Consider a particular measurement error which results in theobservations {x} = { l . O l , 0.99)and { y } = C0.99, 2.01,1.0). In this exampletheeffects of themeasurement errors on the estimate { h } are explored. From (10) the estimate of the impulse responseis obtained as

L(z) =

%z) x(z)

0.99 -

+ 2.012-' + z-'

1.01 -t 0.992-1 = 0.9802

+ 1.02932-1

+

-0.0190~-~ 0.0190~-~

(24)

{z} ={0.9802,1.0293,-0.0190,0.0190,

*-}.

(25)

The measurement errors are seen to cause the element in to differfrom those in { h } , and to cause to be an infinite sequence whereas { h } was finite. By using (24) the expected value of the estimate is found to be

{x}

{z}

Observe that the estimates in (25) are within the 30 of the expected values. A disadvantage with the synthetic division scheme is that the error tends to build up with the lengthof the sequence.

V. A LEAST SQUARES APPROACH TO THE ESTIMATION PROBLEM If X ( s ) [which is the Laplace transform of x(t)] contains a zero in the right half-plane (RHF') [i.e., nonminimum phase], a serious computationaldifficulty may arise with the synthetic division scheme because H(s) may then contain a pole in the RHP [ 181. This results in numerical instability. With a stable systemthe RHP zero of X($)is cancelled by theidentical RHP in Y ( s ) However, in practice,round-off and/or truncation errors in computing X ( s ) and Y(s) prevent this cancellation. The synthetic division approach does not then converge. Should this occur, an alternate approach must be used. This sectiondevelopsthe least squaresestimation of thetime domainimpulseresponsefrommeasuredinputand output data. With reference to (7), theresponse ofadiscretetime system may be writtenas

E[K(z)] = ( 1 + O ~ ) + ( l - u ~ ) z - 1 + 2 0 2 z - - - 2 u z - 3 + - - '

= 1.0001

+ 0.99992-1 + 0 . 0 0 0 2 ~ - ~

- 0.00022-

+ -.

The convolution in (30) may be expressed in matrix form (26)

as

It is seen that the meanvalue of the estimate is biased. This is due to measurement errors in the input. The variances of the estimate are obtained from (21) and (22). Observe that

It follows that

where Y' and H are N X 1and L x 1column matrices, respectively, and X is an N X L rectangular matrix. Equation (31) can besolved for [HI byfirstpremultiplying both sides by [X] (i.e., transpose of [ X ]). Then

'

[XI = V I [HI = [XI =[ Y'I

Consequently the variance of the first element in 20' while it is 3u2 for the remaining elements. Letting k = 3 in (23),Chebyshev's inequality yields

< 3(0.0141)] 2 8/9 P[ I 1 1 - 0.9999 I < 3(0.0173)] > 8/9 P[ I - 0.0002 I < 3(0.0173)] 2 8/9. (29) P[ 110- 1.0001 I

{K } is

-

(32)

It follows that [HI = {[XI qx1 I- '[XI =[ Y'l .

(33)

It is interesting to note that [X] '[X] is the autocorrelation matrix of the available input time samples. Thus [X] [XI is a symmetricpositivesemi-definitematrix. In fact [X]=[X] is Toplitz a matrixandthereforeisinvertible bymeans of a recursiveprocedurerequiring L2 operations. ' 0

66 1

DETERMINATION SARKAR RESPONSE e t al. : IMPULSE

The inversion of [X]*[X] is numerically highly unstable where [PI and as has been pointed out by Ekstrom [ 191. He suggested are defined as a singular value decompositiontechnique for the inversion. In this procedure [X] T[X]is converted to a diagonal matrix bysimilarity a transformation.The eigenvalues arethen filtered (Le., eigenvalues belowacertain small amountare set' to zero). The whole process of premultiplying by [X] and then obtaining the singular value decomposition, is time consuming. In this paper an alternate approach is suggested and based upon the conjugate gradient method [ 201 . With the conjugategradient method it is not necessary to formthematrix [XI*[XI. Thematrixequation [XI [HI= [ Y'] is solved by the following procedure. By using any initial guess [ H I o for the solution, the matrices [R] and [PI aregenerated, where

[e]accountforthe

measurementerrorsand

*,

[ ~ l o = [ ~ l [ ~ l o - ~ ~ ' l

%I -

(35)

The (n

+ 2)th estimate is then given by

[HI n+

1

= [HI n

+ tn[PI

4N

(34)

and

[PI 0 =-[XI

=rJ

(36)

n

(43)

[e]

The elements of [PI and are assumed to be statistically independent zero mean random variables. By assuming the measurement errors are small, an estimate of the impulse response is obtained as

[a] = { [ X + P ] T [ X + P ] } - ' [ X + P ] T [ Y ' +

where

Q]

= i[XlT[X11-'wI -([XITIXl)-l[All

- [ Xf P] *[ Y' + Q] '+ higher order terms

(37)

(44)

where and [R] and [PI are obtained according to the relations

1 [PI n -

1

1 [PI n- 1

[I] = identity matrix, and

(38) (39)

[ A ] = [X]T[P]

+ [PI T[X] + [PI *[PI.

The expected value of the time domainimpulse response is given by E[G]

-{[x]q X 1 l - I

where The advantage of the conjugate gradient method is that it is an iterative scheme which usually yields excellent results within J iterations where J is the number of independent eigenvalues of [ X ] * [ X ] . It has been shown that the conjugate gradient method converges to a solution even when [ X ] * [ X ] is singular. For this case the iteration process is terminated once the solution stops converging (i.e., begins to oscillate). This method has been used successfully in the area of image processing for a similardeconvolution problem [ 2 11.However an analytical expression is available for the total number of iterations beyond which the iterative procedure can be terminated [251. The error analysis performed in the previous section is now repeated in order to find mean and variance of the solution obtained bytheleast squares method whenmeasurement errors are introduced.Corresponding to

[XI [HI = [ Y'I [I?], an estimateof 1 8 3 , satisfies the matrix equation [X+P][g] =[Y'+Q]

[up2] = EI[PI T[pll.

(45)

Also notethatwhen measurementerrors are associated with the input, there is a bias in the solution which is given by H-E'[E?] %{[X]T[X]}-2[~p2][X]T[Y'].

Similarly the variance of [I?] is defined to beacolumn vector whose elements are the variances of each of the elements in [I?].Note that

[ill

- E[R]

= {[XI '[XI

1- '[PI

- ([XI TIXI )--

*[ Y'I + [XI TIQl

'

{[XI T [ p l $- [PI *[XI )[XI

+ higher order terms = [Dl

+ higher order terms.

*r Y'I

662

IEEE TRANSACTIONS ON ANTENNASAND PROPAGATION, VOL. AP-30, NO. 4 , JULY 1982

Then

(47)

From (44) an estimateof the time domain impulse response cat be obtained at each instance through the column matrix [ H I . The expected value and the variance are obtained from (45) and(47),respectively, through the elements of the E [ E ] and u2 [HI.Using Chebyshev’s incolumnmatrices equality, as in the previous section,the probability of a particular ele_ment [ g i ] ,of thematrix [ E ] ,lyingbetween 2 E[E?;.I f k u [ H i l , is greater than (1- l / k 1, is., P{I[Ej] -E[fij]lSku[fij]}>l---

1 k2

(48)

Example 3: Consider, once again, the same problem stated in example 2. From the problem statement, (41) becomes

[;gg1.01

0 0.99 l.ol]

x[;j=[’:a:].

1.0

(49)

After premultiplying both sides of the equations byJX f PI and inverting [X f PI T[X f P I , theestimate [HI becomes

By comparing (50) to (25) we find that the least squares method yields a betterestimatethanthesynthetic division techniqueforthis problem. Also the variance ofthe least squares method is much smaller thanthe variance of the synthetic division method. However thesynthetic division scheme is computationally much simpler than the least squares method and yields excellentresultswhen u of the measurement error is very small.

VII. OPTIMUM INPUTS FOR WHITE NOISE In mostpractical cases theinput is knownwith great accuracy. The significant errorsarethen associated with measurement of the output. The question then arises, “Given p ( z ) 5 0-and q ( z ) # 0, which input x ( z ) minimizes the variances of h ( z ) in ( 2 2 ) and (47)?” It is apparent from (21), (22), and (47) that the variances of g ( z ) are reduced as the amplitude of x ( z ) is increased. In a practical situation, however, theamplitude of the input may be limited in order not to excite system nonlinearities. Also, the input amplitude may be constrained due to power limitations. Therefore the above questionis pursued subject to a constraint ontheinputamplitude.For convenience, the constraint is chosen to be on the mean square value of x ( z ) , 1.e., GxX(O) = autocorrelation of the inputx of lag zero

--x l

-

M

x,. 2

M f 1m = ~ Assuming p ( z ) = 0 and q ( z ) to be a white noise sequence, (47) becomes

(56 )

u2[B] = US*{[X] T[X]}-1.

Observe that [ X ] T [ X ] is a symmetric positive definite matrix with each element along the principal diagonal having the value (M f l)&JO). UsingLevin’s results [ 2 4 ] , the variance of (56) is minimized if and only if

This compares to 0.9802 1.0293

[XI TIX1 = (M + I)Q*X(O)

by the syntheticdivision method. In addition, theleast squares approach yields

where [ I ] is the identity matrix. Hence provided the input is deterministic, it follows that the solution forx ( z ) satisfies

E[Hl = {[XI T I X 1 1- l([II -{[XI TIX1 I-

--

[ u p 2 ])[XI

=I

0.9999

-

0.9999

- 111

Y’I

I

0

GXX(0)

f

$xx(Z)

= 0,

for 0 < I