THE UNLVERSI Y O MPffSAs ENGINtEERING LIPRIRY Technical Report No. 166 TIE: RECOVERY OF CLOCKED SIGINAL INPUTS TO LIITEAR DYNAITIIC SYST.EIS by F. N. Bailey and M. J. Darnborg Approved By: __ _ __ T. W. Butler, Jr. for COOLEY EIECTRONICS LABORATORY Department of Electrical Engineering The University of Michigan Ann Arbor, Michigan Contract No. DA-36-039 AMC-03733(E) U. S. Army Electronics Command Fort Monmouth, New Jersey Septemrber 1965

ACKNOWLEDGMENTS The authors would like to acknowledge the valuable assistance received from their colleagues at the Cooley Electronics Laboratory through discussions of our work and criticisms of the initial drafts of this material. We would especially like to thank Professors Arch W. Naylor and Keki B. Irani for their continued assistance and advice throughout this project. Special recognition also goes to Dr. Frederick M. Waltz for his many helpful ideas concerning circulant matrices. Finally, we would like to thank Miss Carol Wolanski for her intensive efforts in the typing and preparation of this manuscript. ii

E rrata p. 4 line 17: delete "if" before "'knowiedge". p. 5 line 7: change "near" to "linear". p. 15 line 13: delete "shifti;ng" before ".output". p. 16 Eq. 43: should be t-T f w[t - - (s- )] Pi( - ) d(s- T). -o p. 19 and 20 These two pages should be reversed. i. 24 Eq. 69: should be T 2 E S (v (t) 0(t) )2 d}t 0 p. 25 Eq. 73: should be S SllS S a,, N Ni N 2rl 2 p. 31 line 24: Beetwen "relativ'e" and "damaging" insert "to LT, it seems reasonable tha-t for some cases, depending,' uLpon the system, the". p. 49 Eq. A. 6: should be r 2.p L L Ui U'. L i _2 pj Le - J COO-LEY ELECTRONICS LABOiATOR-t!I T echnical Repor t No. 1'6 6137- 1 I-T

TABLE OF CONTENTS Page ACKNOWLEDGMENTS ii LIST OF ILLUSTRATIONS iv ABSTRACT v 1. INTRODUCTION I 1.1 Background 1 1.2 Representation of Clocked Sequences 2 1.3 State Space Representation of Linear Systems 5 1.4 Report Format 7 2. RECOVERY OF PERIODIC CLOCKED SIGNALS WITH KNOWN FILTERS 8 2.1 Introduction 8 2.2 Development of the Inverse Filter: Fundamentals 9 2.3 Development of the Inverse Filter: Transforms.12 2.4 Development of the Inverse Filter: Computation 15 2.5 Performance of the Inverse Filter in the Presence of Noise 21 2.6 Inclusion of A Priori Knowledge About the Coefficients 28 3. RECOVERY OF A PERIODIC CLOCKED SIGNALS WITH KNOWN FILTERS 31 3.1 Problem Representation and Notation 32 35.2 Solution Using Transforms 33 3-5.3 Performance in the Presence of Noise 38 3.4 Approximate Transform Using Sampling 39 4. PERIODIC CLOCKED SIGNALS INTO UNKNOWN FILTERS 43 5- CONCLUSION AND RECOMMENDATION 47 5.1 Conclusion 47 5.2 Recommendation 48 APPENDICES 49 REFERENCES 52 DISTRIBUTION LIST.553

LIST OF ILLUSTRATIONS Page Fig. 1. Block diagram of the system. 2 Fig. 2. A random binary sequence. 2 Fig. 3(a). A typical basis function, L = 10, T = 1. 3 Fig. 3(b). A typical clocked signal, L = 10, T = 1. 3 Fig. 4(a). A typical basis function, L = 8, T = 1. Fig. 4(b). A typical periodic clocked signal, L = 8, T = 1. Fig. 5. The input and output spaces for G and the transforms defined upon them. 14 Fig. 6. Pole-zero plot for the linear system of Eq. 61. 2L Fig. 7. Representation of the output space O, its orthogonal complement 0O and the various signal and noise vecors. P22 Fig. 8. The dependence of ~|kl, k = 1,2, upon the collinearity of z1 and z2. 25 Fig. 9. The linear system for the example of Section 2.5. 27 Fig. 10. Results of example showing DIIlkI vs. L for various filter parameters. 29 Fig. 11. Results of example showing saturation value Iekl s vs. filter parameters 30 Fig. 12. Comparison of the effect of noise in close subspaces (e.g. A' and B') to subspaces not close (e.g. A and B). 36 Fig. 13. Relation between a minimum phase, Wm(s); and all-pass, Wl(s); and nonminimum phase filter, W2(s). 44 Fig. 14. Minimum phase filter Wl(s) and possible non-minimum phase variations. 45 iv

ABSTRACT This report primarily deals with the recovery of a clocked input (see Section 1.2) to a known linear time-invariant system when the output, or at least a noisy version of the output, can be observed. Another topic, receiving minimal attention here, is the identification of the time-invariant linear system when the output can be observed and certain a priori statistical information about the clocked input is available (e.g., the power spectrum). Both of these problems are of interest in data transmission and in communication systems using clocked inputs. The input recovery problem is treated separately for periodic inputs and aperiodic inputs. In both cases, a transform (a type of "inverse filter" different from the standard inverse filter) is developed which operates on the output to recover the input. This transform can be realized in the form of a computer program. When the inputs are periodic, a special type of matrix called a circulant results which simplifies the computations necessary to recover the input for the entire period. When the inputs are aperiodic, the time axis is divided into intervals and the input is recovered one interval at a time. If the output is corrupted with noise, only an approximation of the actual input is recovered, in the noiseless case the actual input is recovered. A preliminary solution is suggested in the brief discussion concerning the system identification problem. This solution uses power spectrum information as well as the transform mentioned above for the periodic case. It has been verified in special situations. However, there is still some question about the general validity of this method.

1. INTRODUCTION 1.1 Background This report concerns an investigation of linear time-invariant filters having clocked signals as inputs (clocked signals will be defined explicitly in Section 1.2). This class of inputs is of considerable interest since clocked binary input signals are the basis for a number of schemes for developing communications systems with desirable properties.l Most of the results reported here apply to the more general case in which the inputs are not restricted to binary signals. The results are developed in terms of a solution to the following problem: Given the output (possibly corrupted by noise) of a linear time-invariant system2 forced by a clocked signal, and given certain a priori information, obtain a description of the linear system and/or the input sequence. Depending upon the assumed a priori information and the desired result, there is more than one problem to be found in this general statement. One such problem is to obtain a description of the linear system from observation of the output and from assumed statistical knowledge about the input (such as the power spectrum). Another problem is to obtain the input sequence from observation of the output assuming that the linear system and certain information about the form of the clocked input are already known. Besides being of interest in the study of communications systems, this problem of obtaining a description of the system or input sequence is of interest in a broader class of problems concerning data recovery. For instance, consider a device with a digital output feeding a data link or, alternatively, consider a digital signal being measured by a transducer. In these and similar problems, the recovery of the original signal from a distorted (filtered) version that is corrupted by noise involves the solution of this problem. The system and the signals involved will be represented by the block diagram of Fig. 1 1Reference 1 presents a study of communications systems using clocked binary signals, and is concerned with the properties of filtered binary signals and their relation to the binary sequences which generate them. 2The results will be developed only for those systems which can be modeled by linear, constant coefficient, ordinary differential equations.

v(t) YMt) + + c(t) Y w(t) (t) Fig. 1. Block diagram of the system. where c(t) is the clocked input sequence; w(t) is the impulse response of the linear timeinvariant system; y(t) is the system output; v(t) is additive noise corrupting the output; and, y(t) is the noisy form of the output. The input signal recovery problem receives major attention in this report. By using a small amount of prior information about the form of the clocked input (information which determines the particular linear space containing the inputs), we can circumvent the wellknown difficulties normally encountered in the design of an inverse filter. The "inverse filter" which is developed is not simply the ordinary inverse filter and actually, it does not operate in real time. However, it is a true inverse filter in that it operates upon the output signal to recover a delayed version of the input. The result basically depends on constructing efficient bases sets for both the input and output linear function spaces. 1.2 Representation of Clocked Sequences In this report, a clocked signal is a sequence of time functions defined over equal adjacent intervals of time called clock periods. The time functions in each interval are identical except for a constant multiplier. The random binary sequence shown in Fig. 2 is an example of this. c(t) 0 27 4T 6T 8T 07T 12T 14T 16T Fig. 2. A random binary sequence. Here the function in each clock period is identically equal to one, and the constant multipliers are random samples from a binary distribution of zeros and ones. In a more general case, the time function defined on each clock period has a basic waveform b(t) for O t < T

a segment of the clocked signal of length T = LT can be represented as L c(t) = ai qi(t) 0 < t < T (1) i=l where L is the number of clock periods in the segment, T is the length of the clock period, and where b(t-[i-l]r) (i-l)T < t < iT, i = Z,2,...,L qi(t) = 0 elsewhere; (2) and ai c R = the set of real numbers. Note that the functions in the set (qi(t)) are linearly independent and thus form a basis for the L-dimensional space of all such clocked signals. This general case is illustrated in Figs. 3(a) and (b). q6(t) 0 1 2 3 4 5 6 7 8 9 10 Fig. 3(a). A typical basis function, L = 10, T = 1. c(t) 0 31 4 _ 7 10 t Fig. 5(b). A typical clocked signal, L = 10, - = 1.

A special case of interest occurs when the sequence is periodic, that is, when the sequence is repeated after L clock periods. An example is a periodic version of the binary sequence shown in Fig. 2. Such a periodic sequence will consist of a sequence of ones and zeros which repeats every T = LT seconds. These signals, sometimes called "pseudo-random sequences," are readily generated by shift-register generators with proper feedback (Ref. 2). The name "pseudo-random" comes from the fact that pieces (O < t < T) of such sequences, while completely deterministic, have the same appearance as the truly random sequence shown in Fig. 2. Since these sequences are commonly employed in digital communication systems, they are particularly applicable to this report. A general periodic sequence of period LT, where the basic waveform (the time function defined on each clock period for 0 < t < T) is b(t), can be represented as L c(t aj pi(t) = o < t< (3) i=l where b(t-[i-l]v-nT) nT + (i-l) T < t < nT + iT pi(t) = 0 elsewhere (4) for n = 0, +1, +2,... and i = 1,2,...,L and where ai c R Like the nonperiodic case described by Eqs. 1 and 2, the periodic functions in the set (pi(t)} are linearly independent and may be considered as a basis for the L-dimensional space of all such periodic clocked signals. This general case is illustrated in Figs. 4(a) and (b). The representations of clocked signals given by Eqs. 1 and 2 or Eqs. 3 and 4 will be used throughout the remainder of this report. If it is assumed that if knowledge of the system implies knowledge of the set of functions (pi(t)} or (qi(t)) (for the periodic or aperiodic cases, respectively), then the input is uniquely specified by the set of coefficients (ail. In many cases, the only quantities of interest in the signal recovery problem are the ai's. We can consider this problem solved when these are found.

p3(t) le T - - I _, I i I T 8 I I i -6 -4 -2 0 2 4 6 8 10 12 14 16 18 20 t Fig. 4(a). A typical basis function, L = 8, T = 1. c(t) -6 -4 -2 2 4 6 1 12 14 m 1 20 t Fig. 4(b). A typical periodic clocked signal, L = 8, T = 1. 1.3 State Space Representation of Linear Systems The basic problem stated in Section 1.1 involves the study of clocked signals in linear, time-invariant systems. Although such systems can be represented either in the time domain or in the complex frequency domain, the time domain -will generally be used since it is a more natural setting for this particular problem. A basic tool in time domain treatment of dynamic system problems is the state space representation (Ref. 3). It has been shown (Ref. 3) that any physically realizable, singleinput, single-output device that can be modeled by near, constant coefficient, ordinary differential equations has a mathematical model of the form. = Mx + b c(t) (5) y(t) = h'x (6) where x is the state vector (a column vector), x is the time derivative of x, M is a constant matrix, b and h are constant column vectors, h' is the transpose of h, c(t) is the scalar input and y(t) is the scalar output. At each instant, the state vector x may be considered to represent information stored in the system (charges on capacitors, currents in inductors, etc.) because of past inputs and initial conditions. Future values of the state vector will

be determined by the present state and by the future inputs. This is clear from the wellknown solution of Eq. (5) (Ref. 4) t x(t;Xoto) = (t-to)Xo0 + O (t-s) b c(s) ds (7) to where x(t;xo,to) is the value of the state vector at time t when the initial conditions were x at time too The matrix X(t-to) is obtainable from the matrix M by various methods (Ref. 4), and is sometimes called the transition matrix (or fundamental matrix), since it shows the transition of xo = x(to) into x(t) for t > to when c O. Equation (7) shows that the present state x(t;xo,to) depends upon the initial conditions xo and the inputs in the interval (to,t). The output y(t) is now obtained as t y(t) = h'x = h' 0(t-to)xo + I h'm(t-s) b c(s) ds (8) to If the linear system has all of its poles in the left-half-plane (which implies that the characteristic roots of M are in the left-half-plane),3 the effects of initial conditions become less important as t-to becomes larger. That is, for fixed t lim (t-t) = (9) and so in the limit as t + - oo 0o t y(t) = [h' (t-s)b] c(s) ds. (10) -00 It is also known (Ref. 5) that y(t) can be represented as the convolution of the impulse response w(t) of the system with the input c(t), so that t y(t) = w(t-s) c(s) ds (11) -00 Hence, the quantity in the brackets in Eq. 10 is, in fact, the impulse response of the system. 3This follows from the fact that the characteristic roots of the system matrix M are the poles of the system transfer function. The existence of (10) for c(s) bounded and of the form of Eq. 1 is assured since the transition matrix ~(t-s) is exponentially bounded.

1.4 Report Format We will proceed with the solution of the problem in these terms: In Section 2 of this report, we consider the recovery of a clocked periodic input signal for a known linear system. A special "inverse filter" is constructed to operate on noiseless or noisy observation of the output to reconstruct the input c(t). The recovery of arbitrary (nonperiodic) input signals for a known linear system is discussed in Section 3. A method, similar to that of Section 2, is presented to handle this case. In Section 4, we consider the problem of determining an unknown linear system from statistical information about its input and from observation of its output. The conclusions reached as a result of this research and our recommendations for further study are discussed in Section 5.

2. RECOVERY OF PERIODIC CLOCKED SIGNALS WITH KNOWN FILTERS 2.1 Introduction The simplest case of the problem outlined in Section 1.1 occurs when the filter is known exactly and the input signal is periodic. This case will be treated in considerable detail since the approach developed in this section provides a background for understanding the treatment of the more complex problems considered in Sections 3 and 4. First, consider the noise-free case where y = y in Fig. 1. From this block diagram representation of the problem it is clear that a possible solution would be an inverse filter of some type to recover c(t) when given y(t). Not only are they often physically unrealizable but they are usually undesirable in a practical sense because of problems with bandwidth and noise. The source of this difficulty is suggested by the following development. Consider the situation where c(t) is any periodic function with finite period T which can be expanded in the Fourier series j 2~k t c t c(t) = X Ck e (12) k= -oo where T 2ik 1 Kj t(Is Ck = T c(t) e T dt (13) o for k=O,+l, +2,.... Equation 13 is the discrete frequency Fourier transform of c(t) and Eq. 12 is the inverse transform. The coefficient Yk in the expansion of the output y(t) is found by noting that Y =H(j T ) Ck (14) where H(jw) is the Fourier transform of the impulse response of the linear filter. The time function y(t) can now be obtained from the inverse transform Co. 2rk y(t) Yke T (15) k= -oo

In the problem of interest here, the output y(t) is given and c(t) is sought. To obtain c(t) in the form of Eq. 13, we require the inverse operation H-i (j) on Yk such that H-1 2k -l 2k 2k C (16) — T k T T Ck = Ck (16) It is usually true from practical consideration that H(jw) + O as w + o so that (see Eq. 16), the inverse operation H-l(jw) is such that H-l(j) + oo as c + + oo. Thus, H-l(jc) has an arbitrarily large bandwidth which results in the large amplification of any noise at high frequencies. This undesirable feature of the proposed inverse filter is sufficient to make it unacceptable as a solution. 2tk Objections can also be raised against using the functions ej T t for the expansion of c(t) as in Eq. 13 since then an infinite number of the coefficients Ck must be recovered to specify a general periodic c(t). This problem can be eliminated by using a priori knowledge of c(t) to define a finite dimensional subspace containing all possible inputs c(t). A set of basis functions for this finite dimensional subspace can be used to expand c(t) using a finite number of coefficients. Finally, a transform analogous to the Fourier transform above, can be defined in terms of this finite set of basis functions. This approach is used below in the case in which it is known a priori that c(t) is a periodic clocked signal with a known basic waveform b(t) (see Section 1.2) and clock period T. The transform is first developed for the noise-free case and then applied to find an approximate solution to the same problem in a noisy environment. The nature of these results is then examined. At this point, it is desirable to establish the following convention: the transform operation to be defined below is not carried out in real time so that c(t) is not actually recovered but rather a delayed version c(t - 8) for some positive i. However, throughout this discussion this recovered input will simply be denoted as c(t) with the delay suppressed. 2.2 Development of the Inverse Filter: Fundamentals It was shown in Section 1.2 that a periodic clocked signal c(t) with clock period T, basic waveform b(t), and period T = LT could be represented as L c(t) = ai Pi(t) (7) i=l where the pits are defined in Section 1.2. When the ai's are arbitrary real numbers these functions form a finite dimensional linear space, the input space, which will be denoted

by Ip. If the linear filter in Fig. 1 has impulse response w(t), then t y(t) = w(t-u) c(t) du. (18) -00 Substitution of Eq. 17 for Eq. 18 and interchange of summation and integration gives L t y(t) = ai /w(t-u) pi(u) du (19) i=lor L y(t) ai Zj(t) (20) i=l where t zi(t) = w(t-u)Pi(u) du (21) -00 Since pi(t) is periodic, it follows that zi(t) is also periodic. Note that the functions zi(t), i=1,2,...,L, span the output space containing the functions y(t) and, if they are also linearly independent, form a basis for the output space. It is shown below that this linear independence of zi(t), i=1,2,...,L, is guaranteed if the system transfer function W(s) exists and is nonzero everywhere on the jw axis. To see the linear independence of the functions zi(t), i=1,2,...,L, under the assumption that W(s) / 0 on the jw axis, assume that the contrary is true. That is, assume that there exists some set of coefficients (ai), not all zero, such that L y(t) = i zi(t) 0 (22) i=l for t E [O,T]. Recall that y(t) can be written as y(t) = | w(t-u) c(u) du (23) -00 where L c(t) = ai Pi(t) (24) i=l 5Since w(t-u) = 0 for u > t, the upper limit on the integration can be taken as 0 rather than t as in Eq. 18. 10

which can not be identically zero because of the linear independence of the pi(t)'s. Expand c(t) in a Fourier series as co jnot c(t) = A e (25) n= -co where T A= - c(t) e dt (26) n T 0o 2it do T and express Eq. 22 using Eq. 23 as co00 00 y(t) = w(t-u) A e du -o on=-oo n=-oo - n= -oo -O = A W(jrwo) te 0 (27) 6or=t. s that for t c [OT]. Since c(t) 0 in [,T], it follows that the coefficients A in Eq. 25 are not all zero. Now Eq. 27 implies th at An W(jnOo) = 0 for every n and since not all the An can be zero, it follows that for some n, W(jnwo) = 0. This result is contrary to the original assumption and thus, the linear independence of the zi(t)'s has been shown. Note that in this argument, for the fixed period LT, it was only necessary that W(s) i O at the discrete points jmno along the ju axis and not on the entire jc axis as was originally assumed. The stronger assumption is still necessary, however, since the critical points jnu0 change with the period T = LT and it is desired that the result hold for an arbitrary period. Since the same coefficients, the ai's, weigh the basis functions in both the input and output spaces, the problem of recovering c(t) from y(t)will be solved if a method can be Interchanging the integration and summation in Eq. 27 can be justified using the Lebesque Dominated Convergence Theorem. (Ref. 6). 11

devised for obtaining the ai's from y(t), that is, a method of finding the "amount" of each zi(t) basis vector in y(t). Such a method can be devised by finding a set of functions Qi(t), i=1,2,..., L, such that (Qi(t), zj(t)) = ij (28) where i6.j is the Kronecker delta function and the inner product is that given by Eq. 31 of the next section. When a set of ai(t)'s is found, each desired coefficient ai would simply be given by ai = (ei(t), y(t)) i=l,2,...,L. (29) One method of finding these &i(t)'s is provided by the transform T defined in the next section. The process of Eq. 29 can be instrumented in a number of ways. One way is to use digital computation which involves quantizing. This approach is like the one considered in Section 3.4 for the aperiodic case. Another method using continuous time filters is to build a filter with impulse response hi(t) = Ri(T-t). Then, by writing out the expression for the inner product of Eq. 29, the coefficient ai is given by T ai= i(t) y(t) dt 0o T = /hi(T-t) y(t) dt (30) 0o which is simply the output of the filter at t = T under the input y(t) starting at t = O with zero initial conditions. 2.3 Development of the Inverse Filter: Transforms Let the input space as restricted in Section 2.2 be denoted by Ip. This space is the linear function space spanned by the set P = (Pi(t) I i=l,...,L). On this space define the inner product (x(t), y(t)) = V x(t) y(t) dt 12

for all x, y c I and the norm x(t)ll = (x(t),x(t))2 (32) for all x c I. In Eq. 31 the bar above y(t) denoting its complex conjugate is included for mathematical completeness even though the inner product will be restricted to real functions here. The elements of P are clearly orthogonal with respect to this inner product. Any c(t) C Ip has a representation L c(t) = E ai pi(t) = a' p(t) (53) i=l where the column vectors a and p(t) have the form a = col[al,...,aL] and p(t) = col[pl(t),..., PL(t)]. It is possible to define a generalized transformation, TI, on Ip so that (cypl) a = TI(c) = [(Pi,,pj)]L- (34) where [(Pi,pj)] is a matrix whose ij element is (Pi,pj). Because of the orthogonality of the members of P, this matrix is diagonal (with all elements on the diagonal nonzero) and, hence invertible. Consequently, the transform TI maps the input space Ip onto an input transform domain space. The inverse transform of the input is given by Eq. 33, or formally as -1 c(t) = TI (a) = at p(t) (35) A mapping similar to TI can be carried out on the output space (denoted by Op) which, recall, is spanned by the set of functions Z = [zi(t) I i=1,2,...,L). This transform on Op. denoted by To, is defined as a = T (y(t)) = [(zizj)] (36) and the corresponding inverse transform is y(t) = To (a) = a' z(t) (37) Here z(t) and [(zi.zj)3 are defined in a fashion analogous to p(t) and [(Pi,p;)]. The inner product on O0 is the same as that on Ip. The functions ai(t), i=1,2,...,L, of Eq. 29 are 13

seen to be imbedded in the transform To. In fact, denoting the ith row of [(zi,z)] -l by g, i(t) is given by z1 i(t) = [. (38) The transformations TI and T,as well as their inverses and the spaces involved, are shown in Fig. 5. Note that here the operation of the system in transforming inputs c(t) into outputs y(t) is represented as a transformation G from the input space into the output space. R R Ti TI T G Fig. 5. The input and output spaces for G and the transforms defined upon them. The existence of To (hence, the existence of the inverse of [(Zizj)]) follows from the fact that the transformation G is invertible when operating upon the class of clocked inputs in the space Ip. That is, G has only a trivial null space in Ip. Indeed, if G has a nontrivial null space in Ip, there would be some nonzero vector of coefficients a such that y(t) = _ z(t) (39) for t c [O,T]. However, this is forbidden by the linear independence of the functions The vector a, obtained using the transforms TI and To, is a vector in a finite dimensional space as compared to the infinite dimensional spaces employed in Eq. 13. Also, in the transform domain (the domain of the coefficients ai), the input and output are related by the

identity operator. This is contrasted with the Fourier transform case where Yk = H(j T )Ck. Finally, observing y(t) for an entire period of the input (i.e., for all t C [t0, to + T] for some t ) we can take the transform defined by Eq. 36 to get the vector a. The desired signal, the input c(t), is then specified by this a vector and can be recovered by using the inverse transform of Eq. 35. Note that this total operation is actually a type of "inverse filter" in the sense that it operates on the output y(t) to produce a delayed version of the input c(t). However, this operation cannot be carried out in real time. At this point, the problem can be solved in the single step a = To(y(t)). If the matrix [(ziz;)l of Eq. 36 is easily inverted, then the inverse filter can be realized in a very simple and yet exact fashion —no approximations having been made. The feasibility of this matrix inversion step is considered in the next section. 2.4 Development of the inverse Filter: Computation As described in Section 2.3, the recovery of c(t) from the noise-free shifting output y(t) can be accomplished by merely taking the transformation in the output space. This operation is (y,zl) a = To(Y) = [(i'zij )] (40) -- 0 _ (YzL) _ From a computational standpoint, the only difficulty in this procedure is the inversion of the matrix [(zi,zj)]. The fact that this matrix is invertible has been discussed in Section 2.3. This matrix is an L x L matrix where L is the period of c(t) in terms of clock periods. Thus, L is usually large and a straightforward attempt to invert [(zi,zj)] is futile for long periods. Fortunately, a closer examination of this matrix will reveal that this difficulty can be circumvented. In developing an inverse for [(zi,zj)], we use the fact that this matrix has the following special structure: (1) it is symmetric, and (2) it is a circulant (Ref. 7). The symmetry is obvious from the fact that (zi,zj) = (zj,zi). A matrix is a circulant if each row contains the elements of the previous row, shifted one element to the right.7 This shifting is performed in such a way that the last element of any given row becomes the first element of the following row. This circulant property for a given L x L matrix [(zi,zj)] with elements fij, can be characterized by the equations 7When the matrix is formed by shifting the rows one element to the left we get what F. M. Waltz calls "left circulants" (Ref. 8). Many of the properties displayed here are not properties of left circulants. 15

fi. = f i,j = 1,2, L-1 fiL = fijl 1 i = 1,2,...,L-1 (41) The fact that [(zizj)] is actually a circulant is easily verified by showing that its elements satisfy Eq. 41. First note that the basis functions, Pi(t), i=l,...,L have the property Pi+l(t) = Pi(t-T) i = 1,2,...,L-1 Pl(t) = PL(t-T) (42) Then from Eq. 21 t Zi+l(t) = w(t-s)pi+l(s) ds -00 t = J w(t-s)pi(s —r) ds -00 t-T =- w([t-T-(S-T)]pI(s-T) d(s-T) — 00o = Zi(t-T) (43) Using this result in the inner product on 0 gives T (Zi+lzj+l) = Zi(t-T)zj(t-T) dt 0 T-T T zi(s) z (s) ds (45) o (zi, Zj) (46) The equality between Eq. 44 and Eq. 45 follows from the periodicity of the functions involved. A comparison of Eq. 46 and Eq. 41 shows that [(zi,zj)] is indeed a circulant. From the definition of (zi,zi) it is clear that [(zizi)] is a real symmetric (hermitian) matrix and, thus, there is a unitary matrix U such that [(zizj)] = u r u* (47) 16

where r = diag [7o,...-yL_1] Yk = an eigenvalue of [(zi,zj)], k=O,1,...,L-l U* = complex conjugate transpose of U. That the eigenvalues are nonzero follows from the existence of To established in Section 2.3. Moreover, since U is unitary, U* = U and the inverse of [(zi,zj)] is simply written as -l-l1 [(Zizj)1 = U FP U* (48) where F- is easily obtained as r = diag [ -... ]. Hence, if the eigenvalues of 70 7L-1 - [(Ziz;)] and if the unitary matrix U can be found, the required matrix [(Zizj)]- is easily obtained from Eq. 48. Except for the special properties of circulants which facilitate the calculation of their eigenvalues and eigenvectors, Eq. 48 would not be a simplification of the inversion problem. Let dj, j=O,1,...,L-l, be the elements of the first row of the matrix [(zi, zj)] so that the dj are related to the fij of Eq. 41 by the expression d. = flj+l j=0,1,2,...,L-1 (49) The symmetric circulant structure of [(zi,zj)] then implies that the dj are interrelated by the expression j=l,2,..., 2 L odd dj = dLj L j=l,2,...,- - 1 L even (50) Letting rk, k=0,1,2,...,L-l, denote the L roots of the equation r = 1, that is 2kkj rk = e r = cos 2kk + j sin 2k k=O,l,....,L-l (51) L L 8The notation diag[yo,...,y*_L- denotes the diagonal matrix with the elements Yon AL 1 along the main diagonal. 17

it has been shown (Ref. 7) that the eigenvalues, Yk, of a circulant are given by 2 L-l k do + dl rk + d2 r k +... + dL l rk k=0,1,...,L-l (52) That yk' k=O,l,...,L-l, is real can be emphasized by rewriting Eq. 52 as 2Trk 4_k 7k = do + 2dl cos L + 2d2 cos (3) 2d cos i(L-l)k if L odd L-l L ir(L-2)k + (-l)k 2dL cos + ( -1) k dL d - L L if L even L 2 2 which is possible by the symmetry of the matrix and the fact that rp and rL-p are complex conjugates. It is a property of real symmetric circulants that the eigenvalues occur in L pairs of equal eigenvalues such that if L is even we get two unpaired and - - 1 sets of L-l 9 paired eigenvalues and if L is odd we get one unpaired and 2 sets of paired eigenvalues. The corresponding eigenvector is given by 1 21 Ui= r2 i=O,l,..,L-l (54) i rL-1 1 1 2 where (L) is the constant necessary to make the Ui of unit norm. This orthonormality of Ui for i = 0,l,...,L- is shown in Appendix A. The eigenvector Ui forms the i+l column of the matrix U (which is unitary because of the orthonormality of the vectors Ui) and, hence, the inversion can be performed by Eq. 48. It is a general property of symmetric circulants that the inverse is also a symmetric circulant. Thus, to specify [(zi,zj)], only the first row (denote it by g ) is needed. Note from the definition of rk(see Eq. 51) the properties 9The authors would like to acknowledge their indebtedness to Dr. F. M. Waltz of the Cooley Electronics Laboratory for much helpful information concerning symmetric circulants and the example appearing at the end of this section. 18

If we define the vector r_ whose components are defined by i = =, i = 0,1,...,L-, then gl can be expressed by gi 1 ( U ) -1 Thus, Eq. 58 gives the explicit calculations for the first row of [(zi,zj)] and its circulant structure gives the rest. It can be shown that the general i, j element of this inverse is given by the equation L-1 i-1 -k -1 1 rk rj (60) [(zizj)l iS = Z k=O It is estimated that by using the above method for finding [(zi,zj)] at most 4L computer storage locations are necessary where L is the number of clock pulses per period of c(t). Comparing this with L locations necessary for a general inversion, we see that the above method will permit much larger sequences to be handled. Time estimates for carrying out this inverse can be determined from the fact that there are approximately 3L2/4 multiplications in this process of inverting an L x L circulant. On an IBM 7090 digital computer it takes at most 30 microseconds to perform a multiplication so that the total time involved should be on the order of (0.9L /4) x 10-4 seconds. This means that if L is one thousand, the inversion time would be less than 25 seconds. Before the inversion technique of Eq. 48 was developed, a computer program was written (using a general inverse) to check the feasibility of solving the signal recovery problem by means of Eq. 36. Several examples were checked and excellent results were consistently obtained. One such example used a linear system with impulse response given by w(t) = 0.385 e 2t -e 3t[2.39 cos 5t - 1.60 sin 5t] (61) The pole-zero plot of this system is given in Fig. 6 below. A seven pulse (L=7) periodic binary input was used and the vector a defining this sequence was a = col. [1, 0, 1, 1, 1, 0, 0]. (62) This vector was recovered by the transform with only a small amount of error which can be attributed to computer round off error. A generalization in the form of the basic function pi(t), i = 1,2,..., L, is possible. Previously they have been assumed to be mutually orthogonal (a result of their definition, Eq. 4). However, for the existence of the transform To it is only necessary that these basis 19

p r = r L-1 L-p rp = r k pk rpL+k = rk are seen to hold where p = 0,1,2,...,L. Using these properties and the rows of U given by Eq. 54, we see that U has the form 1 1 1 1... 1 1 r1 r2 r3 rL-l 1 r2 r4 r6 rL-2 1 r3 r6 r9 rL-3 U (56) NrL 1 rL-1 rL_2 rL-3... r -1 From the fact that U is unitary and symmetric, it follows that U = U* = U where U is the complex conjugate of U. Thus 1 1 1 1 Y70 Yo 70 70 1 rl r2 rL-1 Y1 Y1 71 rl -l 1 1 r2 r4 rL-2 L = N Y2 Y2 Y22 (57) 1 rL-1 rL-2 l 7L-l 7L-l 7L-l 7L-1 where rk is the complex conjugate of rk. Finally then g, the first row of [(zi,zj)]- -l UT U*, is given by L-l L-l i L-l - L- i 1= 2 L-8) 0i=O i=o i=O i=O i2 20

I - -- 2 2 |i _ _3 s-Plane 1-II O-.. -5 -4 1-3 -2 1-l I-2 ~I _i_-4 Fig. 6. Pole-zero plot for linear system of Eq. 2.50. functions be linearly independent since that (plus an assumption on the linear system, see Section 2.2) is sufficient to assure the linear independence of the functions zi(t), i=1,2,...,L. If, in addition to linear independence, the pi(t)zs retain the shift properties of Eq. 42, the matrix [(zizj)] will have the highly desirable circulant structure. 2.5 Performance of the Inverse Filter in the Presence of Noise When v = O in Fig. 1, so that a noise-free version of y(t) is observable, the mapping To(y) results in a complete recovery of the coefficients a of the input signal and, thus, an exact recovery of a delayed version of the input c(t). On the other hand, when only y = y+v is available,the resulting a = To(y) will have an error associated with it. In general, the noise v will not be in Op (the space spanned by the zi(t), i = 1,2,...,L) so that y will not lie in Op. When the transform To is applied directly to y, the vector of inner products (the rightmost vector of Eq. 36) effectively carries out an orthogonal projection of y onto Op. Any component of y in op, the orthogonal complement to Op is eliminated by this projection.10 It is sufficient for the existence of these orthogonal spaces that Op be contained in a complete normed linear inner product space (i.e., a Hilbert space) containing all functions Y. 21

Thus, if y = y + vl + v2 where y + vl C Op and v2 E Op we get'z = y + v1 (63) where a = To(y). The portion of the noise which is in Op i.e., vl, cannot be removed using linear operators since there is no way to distinguish it from the signals. The above idea of the orthogonal projection of y onto Op is displayed in Fig. 7 in which y is the true signal, y is the noisy signal and the noise v is broken up into components v1 E Op and v2 E Op 4L. 0 p /~/ / Zp O Fig. 7. Representation of the output space Op, its orthogonal complement 0p, and the various signal and noise vectors. Here a'z is the signal recovered by using a = To(-) and is seen to be in error by an amount equal to the noise component in Op. It can be shown that To has the effect of choosing a to minimize IIY - a zll, a result which coincides with the intuitive idea of orthogonal projections. It is desirable to have a measure of the amount of noise corrupting the recovery of a signal coefficient ak. This can be studied in terms of a signal to noise ratio in the kth "channel" of the transform To (i.e., the part of To which recovers the coefficient ak). Recall that the functions ~i, i=1,2,...,L, (see Eq. 28) are imbedded in the transform To. That is, the vector a is given by To as (z,Iy) (21,Y) a = To(y)=[(zizj)] ( (Z2'Y) (,y) 64) (zL'Y) 2Ly)2 22

where the inner product is T (~1,y) = ~1S(t) y(t) dt (65) and T = LT. The ~i(t)'s span the same L-dimensional subspace as the zi(t)'s and satisfy the condition of Eq. 28. Let the noise v(t) be zero-mean white noise with a power spectrum of r volts2/c.p.s. (that is, the covariance of v(t) is R(T) = r2 6(T) volts2; 6(T) is the Dirac delta function). It is now possible to find the amount of energy for the interval [O,T] in the space spanned by the function ~k(t). Note that because of the lack of orthogonality of the ~k(t), the sum over k of these component energies will exceed the total signal energy. The component of v(t) colinear with ~k(t) is (v(t), ~k(t)) vk(t) = lI k(t) (66) That is, the noise v(t) is split into components vk(t) colinear with Sk(t) and vk (t) orthogonal to ~k(t) and only the colinear component affects the kth channel. The expected energy of this noise signal in the interval [O,T] is T NIK = E vk(t)dt E [T T(s) kd s ] k(t)dt} 0 Sk(T),T (Ev(v)v(s)) k(s)ds d T T 1 1 (T L VT TVSkdk(s) dT 1 2 f ~k(T) R(T-S)1k(s)di dT (67) ak o Since the convariance R(T-s) is given by r 2(T-S), this becomes 2 2 = r volts - seconds. (68) 25

Note that this is the expected noise energy in the kth channel at the input to the transform To. The noise in the kth channel at the output of the transform To is given by the quantity (v(t), k(t)) so that the corresponding expected noise energy in the interval [O,T] is Nok = E (t)k(t)) dt =- / ~k(T)E (v(T)v(s))ek(s)dT ds dt O O O T T r ~k(T)dT dt o o T r r l i dt 0o =r l kIT. (69) L In the same manner, the signal y(t) = E ak zk(t) is broken into a component yk(t) colinear with ~k(t) and y I(t) orthogonal to ~k(t). The colinear portion is (~k(t),y(t)) yk(t) II 12 k(t) ak 2 k(t) volts (70) which is the only part of interest to the kth channel since the orthogonal portion is eliminated by the operation of Eq. 65. The corresponding energy in the interval [O,T] at the input to To in the kth channel is T SIk =/ Yk t)dt 2 S = -yak / (t)dt Ik o ak IIQk12 T 24

At the output, the signal energy in the kth channel for the interval [O,T] is T Sok = f(~k(t), y(t))2 dt = - a2 dt o 0 = a T (72) Thus, the signal to noise ratio in the kth cahnnel at both the input and the output of the transform To is 2 sk SIk Sok ak (73) Nk NIk Nok k (73) We see that as the norm of ~k increases, the signal to noise ratio becomes poorer. To gain an understanding of how ||ekiI can change, consider the two-dimensional case of Fig. 8. Note that Z and z2 span the same space, but in one case they are orthogonal, while in the other, z2 has a component colinear with z1. To satisfy Eq. 28, ~1 must be orthogonal to z2 and its inner product with zl must be unity. Thus as zl and z2 become more colinear Il|kfl,k = 1,2, must increase to satisfy Eq. 28. Note that when zl and z2 are orthogonal I1|k| will be minimum and, in fact, ikll = 1 i zkl1 z1 zA z2 2 2 Fig. 8. The dependence of IJ|kl, k = 1,2, upon the collinearity of Z1 and z2 To be able to evaluate the signal to noise ratio of Eq. 73 it is necessary to find |l[kl|2 This can be expressed in terms of the eigenvalues of the matrix [(zi,zj)]. As before, denote the first row of the matrix [(zi,zj)]- by the vector g so that (see Eq. 64) fl(t) is given as l(t) -= gj z(t) (74) 25

Then, using the inner product of Eq. 31, IIIK112 can be expressed as 2 = (gt z(t), gl z(t)) T = [ g [z(t)z'(t)dt]1g 0 go [(Z Jzj)]-l (75) where gl is the complex conjugate of the vector gl. When we replace these quantities by their equivalent representations of Eqs. 47 and 59, we have iI2lL = [k11t U] [U r U] KL U I] 1 I r a (76) From the definitions of F and j it then follows that 2 L-l =1 E _1 (77) i=o Since the index k simply represents a shift in time, it follows that L-1 2= 1111 =1 2 l 1 (78) i=o for k = 1,2,... L. It should be noted that the quantity on the right side of Eq. 78 is the same as the diagonal elements of the matrix [(zi,zj))] (i.e., the first element of', see Eq. 58). Thus, once [(zi,zj)]-l is found we also have II kt Example: To become acquainted with the results that can be expected, a simple example was studied and |kIr k was calculated for a variety of sequence length and filter parameters. For simplicity, the filter was chosen to be a single pole RC network (see Fig. 9). The clocked periodic inputs (period Lv) were chosen to be all clocked signals whose basis functions,

pi(t), i = 1,2,...,L, consist of square-topped pulses with pulse width equal to the clock period T. Note that ||k|| will not depend on the particular input chosen but rather only on the subspace involved which, in turn, is defined by the basis functions. R c(t) C z(t) 0'I Fig. 9. The linear system for the example of Section 2.5. For this filter and this set of basis functions, the first row of the circulant [(zizj)] is found to be given by (l-b)(l-b ) do = (ZZ1) RC (1-bL) 1'i- (zl z _) RC(l-b)2 L-i i-2 ( [b + b ] (79) 2( 1-bL) i = 2,3,...,L -T/RC where b = e. The clock frequency was chosen as 1000 c.p.s. (so that T = 0.001 seconds) and sequence lengths of from 2 to 500 pulses per period were used. The resulting values of le|kfl for RC time constants of from 0.002 to 0.024 seconds are shown in Table 1 below. The graphs in Figs. 10 and 11 illustrate these results. Note from Fig. 10 that H k e quickly levels off, i.e., saturates, as L increases and in fact for this example, DOk|I is constant within 0.01% when L is larger than light. Thus, the signal to noise ratio N of Eq. 73 does not change with increasing L once IIIk|| has reached its saturation value lsklls. From Fig. 11 we see that lI|klls increases in a straight line relation with RC over the range tested. It is likely, but not an established fact, that this straight line can be extended to give II~klls outside of the range tested. In any event, an upper bound on the saturation value'IJklls for any combination of T and RC in this example can be obtained from the expression -1/2 kls limit kll < [T - 2RC (-b) ] (80) 27

11 ikll RC L= seconds 2 4 6 8 10 100 500 0.002 158.43 138.38 136.98 136.88 136.87 136.87 136.87 0.004 311.61 270.46 267.51 267.30 267.29 267.29 267.29 0.008 620.56 537.73 531.76 531.33 531.29 531.29 0.012 930.07 805.70 796.71 796.07 796.02 796.02 0.016 1239.93 1073.95 1061.97 1061.10 1061.05 1061.03 0.020 1549.55 1342.05 1327.07 1326.02 1325.93 1325.92 0.024 1859.64 1610.72 1592.68 1591.39 1591.34 1591.26 Table 1. Results of example giving I|kl for various filter parameters and various length sequences. The experimentally determined values of Ik lls were consistently about 60% of the value of the right hand side of this expression. For the case of RC = 0.002 seconds, the term |1~kl2 in the expression for - (see Eq. 75) Nk results in a contribution of approximately -43db in this signal to noise ratio. Though this is poor, it may not be so bad as to negate the usefulness of the transform To. Note that Sk - also depends upon ak and r which will vary with the particular situation. Nk We found that by using the circulant structure of the matrix [(zizi)] a very efficient computer program could be written to find the eigenvalues. When L = 500, the largest sequence for which the example was run, a computing time of omly about 15 seconds was required on an IBM 7090 digital computer to find the eigenvalues of the 500 x 500 matrix [(zi, zj)]. 2.6. Inclusion of A Priori Knowledge About the Coefficients Unitl now the output of the transform To has been considered to be the solution to the problem. When additional a priori information is available with regard to the values of the aI's (say, that the ai's are binary in which case the inputs c(t) no longer form a linear space), this information can be included as a nonlinear filtering operation used in conjunction with the linear filtering which has been discussed. One possible method for performing this nonlinear operation would be the following two step process. First, apply the transform To to the observed signal to obtain the unconstrained a vector. Second, assuming the elements a. to be restricted to some proper subset (say S) of R, choose a as that vector with elements in S which minimizes some distance function such as [ Z (ai-ai ) ]. This constrained estimation problem has been discussed by Rice (Ref. 9). Other types of nonlinear operations could be devised depending upon the a priori information available concerning the ai's. 28

2000 1800 RC = 24 m. s. 1600 1400 RC = 20 m. s. II kL 1200 - RC = 16 m.s. 1000 RC = 12 m.s. 800 600 \ 0 0 —- ~ _ RC = 8 m.s. 400 |- RC = 4m.s. 200 RC = 2m.s. O L I I I I I I I I I I I I 0 4 8 12 16 20 24 28 L 29

1600 1400 1200 1000 800 600 - 400 200 0 I I I I I I 0 0. 004 0. 008 0. 012 0. 016 0. 020 0. 024 Fig. LL. Results of example showing saturation value 12kfls vs. filter parameters. 30

3. RECOVERY OF APERIODIC CLOCKED SIGNALS WITH KNOWN FILTERS When the input clocked signal c(t) is not periodic the signal recovery technique described in Section 2 breaks down. Recall that the transform To(see Eq. 36) for the periodic case requires the output for a full period of the input signal before the transform can be carried out. Clearly this is impossible for an aperiodic input. There are also situations where c(t) is periodic but the period of c(t) is so long that the techniques described in Section 2 are not practical. One is forced to work with pieces of y(t) that do not contain a complete period of c(t). In any case, it is desirable to be able to handle the problem stated in Section 1 for clocked input signals without useful periodic structure. This section considers a method for the recovery of these aperiodic signals from both noise-free and noisy observations of the output when the filter is known (see Fig. 1). A transform somewhat like To is developed for this purpose. Another approach to this problem which was studied, but which will not be developed here, is to use Wiener-Kalman filtering theory (Ref. 10) to estimate the initial conditions at some time to. Once the initial conditions are known and their effect removed, transform very similar to To of Section 2 can be used to recover c(t) for t > to. This method was found to be inferior to the result of the method this section discussed because it is a twostep process: first, estimating, and second, applying a transform. The method in this section on the other hand, is just a one step process of applying a transform. No noticable simplification in applying the transform is achieved by the estimation step. At this point, it should be suggested that an approximate solution to this problem for some situations could be obtained by simply applying the transform To of Section 2 to a segment of the output of length LT seconds. This solution ignores the initial conditions and assumes that the output is periodic with period LT. Assuming that the impulse response will decay rapidly enough relative damaging effect of incorrect initial conditions will only be important in the first part of the LT time interval. Then if the basis functions qi(t) are of a pulse type (such as in Fig. 3a), the approximation of the coefficients, the ails, will be better toward the end of the LT time interval than at its beginning. Except for the first few clock periods, this approximation would be expected to be rather good. 31

3.1 Problem Representation and Notation Consider a section of a clocked signal c(t) that is L clock pulses long starting at time to. Any such section of c(t) can be written as L c(t) = E ai qi(t) to < t < to + LT (81) i-l Now, referring to Eq. 8, define two row vectors _'(t-to) and _'(t-to) as at(t-to) = h'T(t-to) (82) and _'(t-to) = [Pl(t-to),...,IL(t-to)] where t %j(t-to) = 0 h' D(t-s)b qj(s)ds j = 1,2,...,L. to Then, using Eq. 81 in Eq. 8, interchanging the order of integration and summation, and by using the definitions of Eq. 82, y(t) can be represented by L y(t) = a'(t-to)o + X ajPj(t-tO) (83) j=l for t < t < to + LT. Again defining the vector a = col. [al, a2...,aL], Eq. 83 can finally be written as y(t) = a'(t-to)xO + P'(t-to)a to < t < to + LT (84) Note that if the system is observable (Ref. 11), the elements of the vector a(t-to) are linearly indpendent. This can be seen from Eq. 84 by letting a be zero (i.e., the input c(t) is zero over this time interval). The system is observable, if and only if, there is no nonzero xo such that c'(t-to)o - 0 in [to, to + LT). This is precisely the condition required to ensure that the elements of a(t-to) are linearly independent. Observability will not be assumed for the system since, as will be discussed later, linear independence of the element of a(t-to) is not necessary for the recovery of a. We will now establish the linear independence of the elements of the vector p(t-to) by using a result of Appendix B. Referring to Eq. 84, let the initial conditions -o be zero and consider only the forced response. 32

y(t) = P'(t-to)a (85) The elements of P(t-to) are linearly independent if there is no nonzero a such that y(t) 0 in [to,to + LT). We show in Appendix B that a necessary and sufficient condition for an input to give a zero output over this interval is that the input c(t) must be given by the control law h'M x(t) c(t) = - b (86) where x(t) is the state of the system and the other quantities are defined in Appendix B. From the control law of Eq. 86, we see the additional condition that Ao = x(to) = 0 implies that c(to) = O. Substituting this control law into the system representation (Eq. 1 of Appendix B), we see that the system behaves like the unforced linear system [ hb Because of the zero initial conditions, the state x(t) will remain at zero and the control c(t) remains at zero. Thus, y(t) - 0 in [to,to + LT), when xo = 0, implies that c(t) - 0 in this interval which, in turn, implies a = 0. This establishes the linear independence of the elements of P'(t-t ). 0 It is important that the subspace spanned by the elements of a(t-to) be disjoint from the subspace spanned by the element of _(t-to). In general it turns out that these two subspaces are not disjoint. A more complete discussion of this appears in the next section. With a known filter and a given set of basis functions (q.(t)) for the sequence c(t), the vectors a_'(t-to) and W'(t-to) in Eq. 84 can be evaluated. Therefore, since we observe y(t), the only unknowns in Eq. 84 are the vectors 2o and a. The solution of Eq. 84 for xo and a by means of a transform method similar to that of Section 2.3 is discussed below. 3.2 Solution Using Transforms When y(t) is available for observation without additive noise, the solution of Eq. 84 for xo and a is straightforward using a transform similar to Eq. 36. Rewrite Eq. 84 to give y(t) as the scalar product of two n + L dimensional vectors, that is y(t) =_'(t-to)w to < to + Lv (88) 33

where y(t-to) = col. [a(t-t0),9(t-to)] w = col. [xo,a]. Note that the function Yi(t-to), i = 1,2,..., n + L (where Yi(t-to) is the ith element of the vector, y (t-to)) span the output space for t E [to, to + LT). It is now possible to define a transform ToL on this output space to recover the vector w. This transform is formed by taking the inner product of both sides of Eq. 88 with 7Yi(t-to) for all i=1,2,..., n + L, and then multiplying both sides of this set of equations by the matrix [(yi.,Y)] (defined in a manner analogous to [(zizj)] of Section 2.3). Thus, ToL is given by (Y,71) (Y72 W = ToL(y(t))= [(yiyj)] (89) (Y,7n+L) where now the inner product (f,g) is given by t + LT 0o (f,g) = f(t)g(t)dt (90) to Note that the transform ToL operating on y(t) gives the vector of coefficients w (consisting of a and xo). This is precisely the vector necessary to expand y(t) for t c [to, to + LT) in terms of' (t-to) and' (t-to) as in Eq. 84. Only the vector a is required for recovering the clocked function c(t), and x is really unnecessary. It is seen from Eq. 99 that by considering only the last L rows of [(yi,yj)] we could obtain a and not x -1-0 An inverse transform TIL can be defined as was done in Section 2 which will operate on w to reproduce the input c(t) in the interval [to, to + LT). This is simply the result of Eq. 81, that is -1 c(t) TIL (W) = a' q(t) (91) where q(t) = col.[ql(t),...,qL(t)]. A few words must be said about the existence of the transformation ToL. Let X denote the space of all initial conditions x, and let D denote the space of all clocked inputs 34

for t C [to,to + LT) written in the form of Eq. 81. It is then possible to consider the linear system as a mapping GL of the input product space Ip = X D onto the space 0 of all P p outputs in the interval [to, to + LT). That is, the output y(t) for t c [to, to + LT) can be given by y(t) = GL(2,o, c(t)) = 7'(t-to)w. (92) Thus, the question of the existence of the operator ToL (i.e., the existence of the inverse [7i,7j)l1) reduces itself to the equivalent question of whether or not GL is invertible. Note that an equivalent condition to GL being invertible is that the elements of 7(t-to) are linearly independent. Appendix B shows that, in general, GL is not invertible (i.e., that it has a nontrivial null space). However, GL can be made invertible over a suitably restricted class of allowable inputs. To illustrate what this restriction must be let N denote the null space of GL and let S be the subspace of Ip consisting of the class of allowable inputs. If S and N are disjoint then the operator GL giving the output y(t) as GL(f(to,t)) = y(t) (where f(tO,t) denotes the ordered pair (xo, c(t)) C S) is invertible so that the operator ToL does exist. If, however, S and N are not disjoint, there exist elements fl(tot) and f2(tot) in S so that fl(to,t) - f2(to,t) E N and we can never distinguish between fl(tot) and f2(to,t) since they are mapped into the same point in Op by GL. Thus, the condition for invertibility of GL (hence, the condition for the existence of ToL) is that the subspace of allowable inputs S cX D be disjoint from the null space N X D of the transformation GL. In the particular situation considered here, a nontrivial null space N can arise only if the subspace A spanned by the elements of a(t-t0) is not disjoint from the subspace B spanned by the elements of _(t-t0). That is, assuming observability, this is the only way in which the elements of y(t-to) can be linearly dependent. This follows from the previously established fact that both the elements of o(t-to) and the elements of P(t-to) are linearly independent. It should be mentioned that if the system is not observable (i.e., the elements of a(t-t.) are not linearly independent);but if the subspaces A and B are disjoint, the input vector a can still be recovered, even though the initial condition vector -o cannot. In this case, the matrix [(yi,yj)l will be singular but the pseudo inverse (Ref. 12) can still 11Note that two disjoint linear subspaces have the zero element in common, that is, S and N are disjoint, if and only if, S n N = (0).

be used. This will provide a vector w whose x part will be the one of minimum norm of all those vectors yielding the same a'(t-to)x product. Because of the disjointness of A and B and the linear independence of the elements of P(t-to), the resulting a vector will be correct. Assuming that A and B are disjoint, a comment should be made concerning the degree of "closeness" of A and B. In the ideal case in which there is no noise present, this is of no concern; but with any noise present at all (whether from a noisy signal or round off error), this can be important. Note that the signal y(t) is represented by two vectors: x determining an element of A which is the projection of y(t) onto A along B, and the vector a determining the projection of y(t) onto B along A. When operating upon the noisy signal y(t), the projections do not give vectors xo and a representing the time signal y(t). That these resulting errors become worse, relative to the time vectors as the subspace A and B get close, is demonstrated in Fig. 12. Here we consider two dimensional Euclidean space (the time B B' A' /I —- V V V' / 0 t ~~~x a ~~x axa =A X a X aY 0 0 Fig. 12. Comparison of the effect of noise in close subspaces (e.g., A' and B') to subspaces not close (e.g., A and B). dependence being suppressed) where the subspaces A and A' are spanned by a and a' respectively, and the subspaces B and B' are spanned by P and I', respectively. Note that a and P are orthogonal while a' and P' have colinear components, thus, A' and B' are considered closer than A and B (a more general measure of closeness will be defined later). Let y be the signal and v be the corrupting additive noise. In the situation considered in this figure, the error in the initial condition xo is not affected by the closeness of A and B. However, considering the vector a, we see that a' a - >-_. Thus, the noise resulted in a larger relative error when the two subspaces were a' a close. 36

It is now necessary to consider a more general measure of closeness. The subspaces A and B can be considered to be close if the quantity supremum (hl,h2) is almost unity, where hleAh2eB the inner product is given by Eq. 90 and where hi and h2 have unit norm. The subspaces A and B are not close if this quantity is small compared to unity. Note that this measure is a generalization of the idea that as A and B became more colinear, they get close. From the above discussion, we see that in a noisy situation it is desirable that the subspaces A and B are not close. Note that the elements of A consist of decaying exponentials starting at time to. However, most of the clocked signals of which one thinks are highly discontinuous signals like sequences of impulses and sequences of finite length square pulses which do not give elements of B that "look like" exponentials starting at time to. Thus, with these inputs it is unlikely that supremum (hl,h2) will be close to unity. hleAh2eB It is useful to consider some properties of the null space N of the transformation GL. Two properties of N, established in Appendix B, are: (1) if the state of the system, x(t), has dimension n, then N has dimension n - 1, and (2) a necessary condition for an input to be in N is that c(t) have an arbitrary number of derivatives in the interval [to, to + LT). From the second property, it would appear that any discontinuous c(t) would give an input not in N. However, since discontinuous functions can be approximated closely by continuous functions, it is not safe to say in any practical situation that if a c(t) "looks" discontinuous the input is not in N. It is again more meaningful to talk about the closeness of the subspaces A and B. The following question concerning the null space remains open: Can N contain any inputs which are clocked signals of the form of Eq. 81, that is, can the c(t) given by Eq. 85 ever be a clocked signal? In general, the matrix [(Yi,yj)] of Eq. 89 lacks the symmetric circulant structure discussed in Section 2.4 for the analogous matrix [(zi,zj)]. Thus, it is not possible to use the simplified inversion procedure of Section 2.4. Rather, a standard procedure must be used. This presents no real problem, however, since we have the freedom of choosing an L small enough to keep the computation of this inverse within reason. Note that this transform procedure can also be considered as an inverse filter in the generalized sense discussed before. That is, TOL operating on y(t) gives the vector w which 1 defines the input c(t)(by means of TIL). In applying this transform technique for the case of aperiodic inputs, the user has a certain amount of flexibility in his approach. When considering the first Lr time interval of interest, the transform ToL must be applied directly to find both the initial conditions xO and the input vector a. However, when considering the second LT time interval and all

the succeeding intervals, either the transform ToL can again be used or the knowledge gained from the past interval can be used to construct the initial conditions for the next interval. To be more specific about this alternative procedure, let to denote the beginning of the first LT time interval so that ti = to + iLT, i = 0,1,2,..., is the beginning of the i + 1 st time interval. Now, to recover c(t) for the i + 1 st time interval we use the initial conditions, 2o i, and the a vector for the ith interval (which were found previously) to construct the initial condition -0 oi+l' for the i + 1 st interval. The effect of these initial conditions, x.i+l, can then be subtracted from y(t) to give (see Eq. 83) L r(t) A y(t)- (t-ti)o i+l = ajj(t-ti) (93) j=l The problem now reduces to recovering the ai's from r(t) where the Pits are given (see Eq. 82). Note the functions % (t-ti),j = 1,2,...L, span the space of all r(t) for ti < t < ti+l so using this finite spanning set a transform operating on r(t) can be defined in the same manner as To of Eq. 36 to give the vector a. In defining this transform replace zj(t) of Eq. 36 by j(t-ti) and the inner product to be used is like that of Eq. 31 but with the integration taken over the interval [titi+l]. Note that, this method of updating the initial condition has the undesirable effect of allowing errors to accumulate and it is expected that applying the transform TOL directly to each time interval will give more dependable results, especially if there is an appreciable amount of noise present. 3.3 Performance in the Presence of Noise Suppose that only a noisy version of the output y(t) = y(t) + v(t) (see Fig. 1) is available for observation. In general, the noise v(t) will have a component which lies in output space spanned by the functions yi(t-to0),i = 1,2,...,n + L, so that applying TOL directly to y(t) will result in a vector w which is in error. Following the same argument as that of Section 2.5, if the transform TOL is applied to the noisy signal y(t), all noise orthogonal to the output space is removed. Thus, again, an estimate of the input is obtained by projecting y(t) orthogonally onto the output space and finding the vector input w which would result in this projected output. The underlying ideas of Section 2.5 concerning signal-to-noise ratios could be carried over to study an analogous concept here. However, compared to Section 2.5, the calculations involving I[kfl would be complicated by the lack of circulant structure in the matrix

The comments of Section 2.6 concerning procedure when given a priori information about the vector a apply directly to this case. 3.4 Approximate Transform Using Sampling To simplify the computation involved in applying the transform ToL, a method which uses sampling to approximate TOL can often be used. This is carried out in the following manner. Sample y(t) at the set of times (ti) where ti C (to,t0 + LT) for i = 12,... m > n + L. The m equations formed in this way can be put in the matrix form y Hw (94) where w is defined as in Eq. 88 and where: y = col.[y(tl)Y(t2),...,y(tm)] Y'(t1 - to) 7'(t2 - to)' (tM - t ) Recall that y'(t-to)is a row vector with n + L elements. It is necessary to establish that, whenever the transform ToL exists, the set of sampling times (ti) can be chosen in such a way that the matrix H is of maximum rank. Recall that existence of TOL implies that the elements of the vector valued function 7(t-to) are linearly independent, that is, whenever a' Y(t-to) O (95) for every t E [to,to + LT) then a = 0. It is sufficient to let m = n + L and show that there does exist a set of n + L sampling times (ti) such that the n + L vectors in the set (Y (ti-to) i = 1,2,...,n+L), are linearly independent. Note that the vectors in this set are linearly independent, if and only if, they span the n + L dimensional Euclidean space En+L. Assume no such set of sampling times can be found. Then the values of the vector valued n+L functions 7(t-to) lie only in a proper subspace, say Q, of E for all t E [to,to + LT). n+L It is then possible to choose a vector r which lies in E and is orthogonal to Q. That is r'y(t - to) 0 (96) for every t C [to, to + Lv). However, this is contrary to the linear independence of the 39

elements of Y(t-to0) so that there does indeed exist a set (t,) of n + L sampling times such that H is of maximum rank. It should be pointed out that in practice such a set (ti) will not be difficult to find and, in fact, a great many such sets exist. The particular method of choosing the set (ti) will depend upon the situation but, in many cases, a suitable set can be found by picking at least one ti in each of the L clock periods in the interval [to,to + LT) until n + L times ti have been chosen. Note that for vectors formed by sampling, such as y above, the operation analogous to the inner product of Eq. 90 is the scalar product <f, g = f' g (97) where f = col.[f(tl),f(t2)...f(tm)] g = col.[g(tl),g(t2),...g(tm)] It is evident from this definition of <f, g> that replacing the inner product (f,g) in Eq. 89 by <f, g>, the analogous transform to ToL of Eq. 89 is ToLS given by w = TLS () = (H'H)H'y = Hy (98) Here H is the standard pseudo inverse (Ref. 12). It is known that H as given above exists whenever H is of maximum rank so that the operation of Eq. 98 can be carried out whenever TOL exists. If the noise free signal y(t) is observed, it is sufficient to let m = n + L. The matrix H is then square and, since the set of sampling times (ti9 can be chosen such that H is of maximum rank, it is invertible. We then obtain w directly as w=H-1 y (99) However, if only 5(t)(the noisy signal) can be observed, it is desirable to introduce redundancy into the sample vector y by letting m > n + L. Hence, a signal y(t), describable by the n + L dimensional vector w, is represented along with corrupting noise in an m > n + L dimensional space. This provides added information to aid in making a better estimate of the true y(t), i.e., the true w. When this is done, w can be estimated by means of the pseudoinverse H+ to get w_ = H+y (100) 40

which has the effect of projecting the m dimensional representation of y(t) onto an n + L dimensional space. This is an orthogonal projection so that all the noise outside the n + L dimensional space spanned by the vectors H w (for all w) is removed, the only component left being that in the space spanned by the H w's. Note that because of the orthogonal nature of the projection, V = H+ y is the vector which minimizes I|| - H w|l over all w. If digital computation is to be used, this sampling method will most likely require considerably fewer computations than the full transform ToL. This is due to the ease of calculating the sequence inner product <f,g> compared to the integral inner product (f,g) of Eq. 90. However, there are other methods of forming (f,g), such as by analog computation. This same type of sampling technique could also be used in the periodic case considered in Section 2. Note that this sampling method is actually an approximation to the transform To or TOL (of Section 2.3 and 3.2, respectively) and, in noisy situations, this sampling technique will give poorer estimates than the original transform To or ToL. To see that the result using sampling is poorer when noise is present, compare the two methods in the following way (done here for the periodic cases To and ToS). Recall that using the transform To the coefficient ai is obtained by the expression (see Section 2.2) a. = (y(t), ~i(t)) (101) In a similar manner, the sampled transform ToS gives ai by the expression ai = (S y(t), Li) (102) where S is a sampling operator and the Li are vectors such that (S z (t), Li) = bij (the Kronecker delta) for ij = 1,2,..., L. Using the adjoint operator of S, Eq. 102 can be expressed as ai = (y(t),S* Li) = (y(t), ~i(t) + ~iL(t)) (103) where ~i(t) is the component of S* Li in the subspace Op spanned by the functions zi(t), i = 1,2,...,L, and QiL(t) is the component in the orthogonal complement of this subspace. Note that the ~i(t) of Eq. 103 are identical to the ~i(t) of Eq. 101. Also note 12The norm Ilxll of a vector x is defined as L1 xm 1/2 41

that, in general, IiL(t) is not zero for if it were then S* Li E Op However this would, in turn, imply that there was no function z (t) in 0O (the orthogonal complement of Op) such that (S z (t), Li) i 0 (104) and, indeed, such a z (t) can be found. This is done simply by picking its values at the sampling times such that Eq. 104 holds (an easy way to do this is to let z (t) = zi(t) at these sampling instants) and then constructing z (t) between these sampling instants such that z (t) Op. Now consider the result when these two operators are applied to a noisy signal y(t) = y(t) + v(t) where y(t) is in 0 and where v(t) has a component vl(t) in 0 and another component v2(t) in Op. Using To we get an approximation ai of the true a. as ai = (y(t), ~i(t)) (y(t),~i(t)) + (vi(t),~i(t)) ai + (vl(t),~i(t)) (105) A The resulting approximation ai when using ToS is i = ((t),(t) + +il(t)) = (y(t),~i(t)) + (vl(t),~i(t)) + (v2(t), i(t)) = ai + (vl(t),Ii(t)) + (v2(t),iL(t)) (106) Thus, we see that the result, ai, using TOS has an additional term corrupting the true ai so that ai is a better approximation. 42

4. PERIODIC CLOCKED SIGNALS INTO UNKNOWN FILTERS Another aspect of the basic problem outlined in Section 1.1 involves the determination of a description of the linear system itself from a priori information about the input and observation of the output. To date this problem has not received a great deal of attention. However, preliminary results presented below suggest one possible approach to this problem when the input clocked signal is known to be periodic. The procedure developed is based on determination of the squared modulus of the transfer function of the unknown system from power spectra of the input and output signals. The inverse filter developed in Section 2 is then used to find the actual transfer function corresponding to the squared modulus function. Referring to Fig. 1, let Oc(w) be the power spectrum of the periodic input signal, c(t) (with period T), and let by(X) be the power spectrum of the noise-free output, y(t). That is, 0c(M) and by(X) are line spectra whose amplitudes at the discrete frequencies n Co = n n = 0, ~1, 2,..., represent the power in the signals c(t) and y(t), respectively, at those frequencies. The square of the modulus of the system transfer function W(s) is then given by (Ref.13)13 2 (my(m) W(jm)l- = (107) For a binary clocked periodic signal the spectrum envelope is independent of the particular input signal and determined entirely by the basic waveform, b(t), and the clock period v, (see Eqs. 2 and 4). For example, if c(t) is a periodic binary sequence made up of positive and negative square pulses (i.e., if b(t) is a unit pulse) with clock period T, then it is easily shown (Ref. 14) that the envelope of Oc(m) is Envelope sin( (108) where K is a constant depending on pulse amplitude. Moreover, if T is not known a priori, it can be obtained from observation of y(t) using some of the techniques suggested in Ref. 14. Note that Eq. 107 holds regardless of whether the input is perodic or not, however, the following discussion concerning Oc(m) depends upon the periodicity of c(t).

The output power spectrum (or its envelope) can, at least in principle, be obtained from the observation of y(t). Methods for doing this are discussed in Ref. 15. If only y = y + v is observable and v is white noise, then (m) is simply y(X) plus a constant In practice, techniques for the recovery of by(X) from observations on y will depend on the signal-to-noise ratios encountered. Assuming that by and (c have been obtained, the squared modulus of the system transfer function W(s) can be obtained from Eq. 107 (up to an unknown multiplicative constant). If it is assumed that the transfer function is a realizable physical network,then the poles of W(s) corresponding to the known IW(ji)12 are uniquely defined in the left half plane. However, there still remains an unresolvable ambiguity in the location of zeros of W(s). This ambiguous W(s) can be represented as a minimum phase version Wm(s) plus an unknown all-pass section (see Fig. 13). This is seen to be simply an ambiguity in the sign of the real part of the zeros of W(s). To completely determine the transfer function W(s), it is necessary to remove the ambiguity in zero locations, that is, to remove the "all-pass ambiguity" noted above. This all-pass ambiguity can possibly be removed in the following way. If c(t) is a periodic clocked signal with period T = LT and known basic waveform, b(t), then for each possible combination of Wm(s) with all-pass sections, as suggested in Fig. 13, there will be only one inverse filter in the sense of Section 2.3. If, in addition we know a priori thee o d x do do o +d -c -b -a -b b -c -a b o -d x -d o -d o Wm(S) Wl(s) W2(s) Fig. 13. Relation between a minimum phase, Wm(s); an all-pass, Wl(s); and a nonminimum phase filter, W2(s). the set containing the coefficients, ai i=1,2,...,L (recall c(t) = a'p(t), Eq. 33), one might hope to select the (hopefully unique) inverse filter which produces an output of ai's in the given set. The combination of Wm(s) and all-pass sections corresponding to this inverse filter would then be selected as the proper transfer function W(s). This approach would be of value if the following conjecture could be proved for the set of coefficients in a particular situation: The set of coefficients produced by the inverse filter 44

a = To(y(t)) (see Eq. 36) will lie in the correct coefficient set only if the inverse filter To corresponds to the actual filter W(s) used in producing y(t). To date, attempts to prove, disprove, or obtain more general or specific true versions of this conjecture for particular coefficient sets have failed. However, in a sample problem involving a binary clocked signal and a relatively simple filter, this conjecture has always been found to hold true. The filter and input signals used were those described earlier in the example of Section 2. In all cases, regardless of whether the minimum-phase filter or one of its nonminimum-phase variations was used as the true filter, the above method did pick out the correct filter and recover the sequence. Figure 14, shows the pole zero plot for the minimum-phase filter Wl(s) and its possible nonminimum phase variations W2(s) and W3(s), Table 2 gives the computer output for the three inverse filters corresponding to Wl(s), W2(s), and W3(s) where W2(s) was used as the true filter. From this, it is seen that the output of W2j1(s) (the inverse filter in the sense of Section 2.3 for W2(s)) is much more nearly binary than that of W-1 (s) or W3 (s). Note that this use of the information about the coefficients is really a form of nonlinear filtering. X 5 X 5 X 5 0 2 2o 0 2 -4 -3 - -3 -3 1 -3 -2 -1 4 o -2 o- 2P -2 X -5 X -5 X -5 W1(s) W 2(s) W3(s) Fig. 14. Minimum phase filter W1(s) and possible nonminimum phase variations.

Inverse Filter Outputs True 1 True 1 1.279709 x 10ol 9.999999 x 10o- 1.436202 x 100 0 -3.062445 x 10'o 1.564622 x 10-7 -1.122676 x o-1 1 -7.172019 x 10-2 9.999998 x 10o-l 1.356072 x 10~ 1 2.088381 x o10-2 1.000000 x 10~ 1.288630 x 10~ 1 1.914498 x 1o-l 9.999998 x 10-1 1.493220 x 100 0 -8.180243 x 10- 2 8.940697 x 10-8 2.083545 x 1o-2 0 -3.517730 x 1o-l -2.980232 x 10-8 2.621020 x 10-1 Table 2. Comparisin of inverse filter outputs with the true sequence where W2 (s) corresponds to the true filter. 46

5. CONCLUSION AND RECOMMENDATIONS 5.1 Conclusion It was shown in Section 2 of this report that when a known linear system is driven, by a periodic clocked signal (as given by Eq. 3), this input signal can be recovered by observation of the output. A transform, or "inverse filter," To(y(t)) (see Eq. 36) was developed to perform this task. It was shown that To exists provided the system transfer function has no zeros on the jw axis. This transform can be constructed in the form of a computer program and the length of period (in terms of clock pulses per period) that can be handled depends upon the size of available computing equipment. To apply To, it is necessary to compute the inverse of a particular matrix, a problem greatly simplified by the symmetric circulant structure of this matrix. It was found that using this structure, a 500 x 500 matrix, could be inverted in less than 25 seconds. Much of the effectiveness of this transform method comes from the fact that the transform is not carried out in real time. The method is restricted to cases in which the output (or a noisy version of the output) can be observed over an entire period of the input signal. If the linear system output can be observed without appreciable additive noise, the transform can be applied directly provided the sequence length is short enough for the available computing equipment. When there is additive noise in the observation, the transform can still be applied to give an estimate of the true input. The problem of recovering a clocked input signal to a known linear system when the signal is aperiodic is considered in-Section 3. The method developed for accomplishing this is an extension to aperiodic signals of the transform technique which was developed in Section 2 for periodic signals. This new transform, ToL(y(t)) (see Eq. 89), recovers a segment of the clocked input signal14 from observation of a segment of the output. The existence of TOL depends upon the inputs satisfying certain necessary and sufficient conditions. As with the transform of Section 2, ToL does not operate in real time and can be constructed in the form of a computer program. When the system output is observed without noise, the transform is applied directly to recover the input. When noise is present, the result obtained by applying 1For a description of such an input segment see Eq. 1. 47

ToL removes all noise orthogonal to the linear space containing the output to give an estimate of the true input. This transform (and that of Section 2) can be approximated by using a sampling technique. In Section 4 of this report, it was noted that the square of the modulus of the transfer function of an unknown linear system driven by a clocked signal can be determined, if the power spectrum of the input sequence is known or can be found from knowledge of the input, and if the power spectrum of the linear system output can be measured. Assuming the system to be realizable, the square of the modulus of its transfer function determines the poles and zeros of the system except for a sign ambiguity in the real part of the zeros. There is evidence to indicate that in some cases the correct zeros can be found by an application of the inverse filter of Section 2. 5.2 Recommendations Recall that, in the periodic case, having a large number of clock periods per signal period does not present any actual problem. That is, if L is so large that using To is unreasonable, the period can be broken into segments and TOL applied to each segment. Thus, in practical situations, signal recovery can be carried out using the transforms To and TOL (given by Eqs. 36 and 89, respectively) providing noise doesn't render the estimates of the input unacceptable. The authors suggest a more thorough study of the effects of noise to determine relationships between the noise level and the quality of the estimates resulting from To and ToL. The method discussed in Section 4 for the determination of an unknown filter from certain information about the input and observation of the output requires further investigation in two areas. The feasibility of finding the power spectrum of the system output should be studied as well as particular methods for doing this. It is also necessary to verify the validity of the suggested method of finding the correct zeros of the system. Finally, another suggestion for future work is the additional development of nonlinear filtering methods to include a priori information about the allowable values of the coefficients ai. In the present paper, this idea was simply introduced and only a minimum of investigation has been conducted on this topic. 48

APPENDIX A PROOF OF THE ORTHONORMALITY OF THE VECTORS Ui To show the eigenvectors of Eq. 54 to be orthonormal, it is sufficient to show 1 for i = k Ui* Uk = O for i f k i, k = 1,2,...,L (A.1) where Ui* denotes the complex conjugate transpose of the vector Ui. From Eq. 51 it is seen that 2ikj L rk = e k = 0,1,2,...,L-1 (A.2) (where j = rT) and putting this in Eq. 54 gives 2F_(k-i)j 4t(k-i)j 2(L-1)i(k-i)j U* UkL 1 + e + e +... + e L (A3) U i kL L L for i,k = 1,2,...,L. For the case when i = k, each term in the brackets on the right side of Eq. A.3 is unity, and then, Uk* Uk = 1 k = 1,2,...,L. (A.4) For i f k, let p = k - i and Eq. A.3 becomes the geometric progression F /2tpj / 2tpj \2 2 pj U* Uk = + + ( +... + L (A.5) The sum of this finite term geometric progression is known so that Eq. A.5 becomes i k =L II 2rp;3) L j (A.6) Since p is a nonzero integer, the numerator of Eq. A.6 becomes e2~pj - 1 = 0. (A.7) Equations A.4 and A.7 show Eq. A.1 to hold and hence, the eigenvectors are orthonormal. 49

APPENDIX B DEVELOPMENT OF THE NONTRIVIAL NULL SPACE OF THE MAPPING GL Consider the linear time-independent system described by the equation x(t) = Mx(t) + b c(t) y(t) = h' x(t) (B.lj as a mapping GL from the input space onto the output space. For the aperiodic case of Section 3, to which this applies, the input space is the direct sum of the space of initial conditions and the space of all clocked inputs c(t) in the interval [o,LT) which are represented by Eq. 1. The output space consists of all system outputs on the interval [O,LT). The existence of the nontrivial null space N of GL is to be shown first and then some of the properties of N will be discussed in greater detail. For a system modeled by Eq. B.1, the output for t E [0, LT) is given by t y(t) = h' [eMtx + eM(t- b c(T)dT] (B.2) - GL(o, c(t)) Recall that w is the vector consisting of the initial condition xo and the input coefficient vector a. To show that GL is not invertible (i.e., that it has a nontrivial null space) it is sufficient to show the existence of a nonzero input w such that y(t) 0- for t E [O,LT). Necessary and sufficient conditions for y(t) 0- in [ O,LT) are that y(t) 0- in [O,LT) and that y(O) h'x = 0 (B.3) From Eq. B.1 these conditions are seen to imply that y(t) = h'[Mx(t) + b c(t)] 0- (B.4) for t a [O,L-) which, in turn, gives the control law required to establish this identity. That is, if Eq. B.3 is satisfied and if the input c(t) is given as c(t) = h M x(t) (B.5) h'b

then y(t) 0 O for all t c [O,LT). Note that with this form of input the system behaves like the unforced linear system X' b ] M ~(t) = M- L-] x(t) (B.6) h' b If the initial condition xo = x(O) = 0, the state of the unforced linear system remains at zero, i.e., x(t) = 0 for t > 0. Under this condition c(t) = 0 for t > 0 so that w = 0 and this is simply the trivial null space of GL. However, if the initial condition O = x(O) i 0 if follows that x(t) f O for t > O so that, assuming observability [Ref. 11], c(t) 8 O for t > O. In this case, assuming c(t) can be expressed as a clocked signal (such as in Eq. 1), w i 0 and the existence of a nontrivial null space of GL has been established. This is also seen to be the only null space of GL since the conditions of Eqs. B.3 and B.4 are both necessary and sufficient. As a side comment, it should be noted that if the system is completely observable and if c(t) 0 O on [0,LT) then whenever x 0, there does exist some t 1(OLv) such that y(tl) i O. Now, consider some of the properties of the null space N of the operator GL. We want to show that if the system state vector x has dimension n, the dimension of N is n - 1. Note that for every initial condition x satisfying Eq. B.3 we get an element of N from -o Eq. B.5. Then, since the subspace of n-dimensional Euclidean space which is orthogonal to h (i.e., which satisfies Eq. B.3) has dimension n - 1, N must also have dimension n - 1. Another property of N is: a necessary condition for an input to be in N is that c(t) be in Cc (i.e., c(t) has an infinite number of derivatives). That c(t) must be continuous in the interval (O,LT) follows from Eq. B.5 and the fact that x(t) must be continuous (since it is the solution of the differential equation representing the system). The continuity of c(t) in the interval (0,LT) follows by differentiating Eq. B.5, that is ~(dt) = -h' M _(t) (B.7) h' b Then the continuity of x(t) as seen from Eq. B.1 establishes the continuity of c(t) and, hence, the existence of c(t) in (0,LT). This argument can be carried on for an arbitrary number of steps, thus c(t) E Co. From this necessary condition, we see that when GL operates on a class of inputs not in C, the operator is invertible, i.e., the null space N is disjoint from the input space. 51

REFERENCES 1. Birdsall, T. G. et al. An introduction to pseudo-random systems. Vols. I & II. Ann Arbor: Cooley Electronics Laboratory Technical Reports No. 104-I & 104-II, The University of Michigan, August 1960. (Secret) 2. Birdsall, T. G. & Ristenbatt, M. P. Introduction to linear shift-register generated sequence. Ann Arbor: Cooley Electronics Laboratory Technical Report No. 90, The University of Michigan, October 1958. 3. Zadeh, L. A. & Desoer, C. A. Linear system theory: The state space approach. New York: McGraw-Hill, 1963. 4. Coddington, E. A. & Levinson, N. Theory of ordinary differential equations. New York: McGraw-Hill, 1955. 5. Davenport, W. B. & Root, W. L. Random signals and noise. New York: McGraw-Hill, 1958. 6. Natanson, I. P. Theory of functions of a real variable. New York: Frederick Ungar, 1961, p. 161. 7. Bellman, R. Introduction on matrix analysis. New York: McGraw-Hill, 1960. 8. Waltz, F. M. Some properties of circulants. Ann Arbor: Cooley Electronics Laboratory Technical Memorandum No. 96, The University of Michigan, September 1965. 9. Rice, J. R. "Approximation with convex constraints," J. Soc. Indust. Appl. Math., Vol. II, No. 1, March 1963. 10. Kalman, R. E. & Bucy, R. S. "New results in linear filtering and prediction theory," Trans. ASME, J. of Basic Engr., March 1961, p. 95. 11. Kreindler E. & Sarachik, P. E. "On the concepts of controllability and observability of linear systems," IEEE Trans. on Automatic Control, Vol. AC-9, No. 2, April 1964, 129. 12. Penrose, P. "A generalized inverse for matrices," Cambridge Philos. Soc. Proc., Vol. 51, 1955. 13. Laning, J. H. & Battin, R. H. Random processes in automatic control. New York: McGraw-Hill, 1956. 14. Ristenbatt, M. P. "Investigation of narrowband waveforms generated by clocked pulses. Ann Arbor: Cooley Electronics Laboratory Technical Report No. 112, The University of Michigan, October 1960. 15. Lee, Y. W. Statistical theory of communication. New York: John Wiley and Sons, 1960.

(U) DISTRIBUTION LIST Copy Copy No. No. 1-2 Commanding Officer 11 Headquarters U. S. Army Electronics Command USAF U.S. Army Electronics Laboratories Washington 25, D. C. Fort Monmouth, New Jersey ATTN: AFRDR ATTN: Senior Scientist Electronic Warfare Division 12 AFAL (AVWW/ECM Technology) 3 Commanding General Wright-Patterson Air Force Base U. S. Army Electronic Proving Ground Ohio 45433 Fort Huachuca, Arizona ATTN: Director 13 Commander Electronic Warfare Division Aeronautical Systems Division Wright-Patterson Air Force Base 4 Commanding General Ohio U.S. Army Materiel Command ATTN: ASAPRD Bldg. T- 7 Washington 25, D. C. 14 Commander ATTN: AMCRD-DE-E Aeronautical Systems Division Wright-Patterson Air Force Base 5 Commanding General Ohio U. S. Army Materiel Command ATTN: ASNP Bldg. T- 7 Washington 25, D. C. 15 ESD (ESTI) ATTN: AMCRD-RP-E L. G. Hanscom Field Bedford, Massachusetts 6 Commanding Officer Signal Corps Electronics 16 Commander Research Unit Rome Air Development Center 9 560th USASRU Griffiss Air Force Base P.O. Box 205 New York Mountain View, California ATTN: RAYLD 7 U.S. Atomic Energy Commission 17 Commander 1901 Constitution Avenue, N. W. Air Proving Ground Center Washington 25, D. C. Eglin Air Force Base ATTN: Chief Librarian Florida ATTN: ADJ/Technical Report Branch 8 Director Central Intelligence Agency 18 Chief of Naval Operations 2430 E Street, N.W. EW Systems Branch Washington 25, D. C. OP-35, Department of the Navy Washington 25, D. C. 9 U.S. Army Research Liaison Officer MIT-Lincoln Laboratory 19 Chief, Bureau of Ships Lexington 73, Massachusetts Code 691C Department of the Navy 10 Commander Washington 25, D. C. Air Force Systems Command Andrews Air Force Base 53

(U) DISTRIBUTION LIST (Cont.) Copy Copy No. No. 20 Commander 32 President Bu Naval Weapons Code RRRE-20 U. S. Army Airborne and Department of the Navy Electronics Board Washington 25, D. C. Fort Bragg, North Carolina 21 Commander 33 U. S. Army Anti-Aircraft Artillery Naval Ordnance Test Station and Guided Missile School Inyokern Fort Bliss, Texas China Lake, California ATTN: ORL ATTN: Test Director - Code 20 34 Commander 22 Commander USAF Security Service Naval Air Missile Test Center San Antonio, Texas Point Mugu, California ATTN: CLR 23 Director 35 Chief of Naval Research Naval Research Laboratory Department of the Navy Countermeasures Branch, Code 5430 Washington 25, D. C. Washington 25, D. C. ATTN: Code 427 24 Director 36 Commanding Officer Naval Research Laboratory 52d U. S. Army Security Agency Washington 25, D. C. Special Operations Command ATTN: Code 2021 Fort Huachuca, Arizona 25 Director 37 President Air University Library U.S. Army Security Agency Board Maxwell Air Force Base Arlington Hall Station Alabama Arlington 12, Virginia ATTN: CR-4987 38 The Research Analysis Corporation 26 Commanding Officer-Director McLean, Virginia 22101 U. S. Navy Electronics Laboratory ATTN: Document Control Officer San Diego 52, California 39-48 Headquarters 27 Commanding Officer Defense Documentation Center U. S. Naval Ordnance Laboratory Cameron Station Silver Spring 19, Maryland Alexandria, Virginia 28-30 Chief 49 Commanding Officer U.S. Army Security Agency U.S. Army Electronics Research and Arlington Hall Station Development Laboratory Arlington 12, Virginia 22212 Fort Monmouth, New Jersey ATTN: 2 Cpys - IADEV ATTN: U.S. Marine Corps Liaison 1 Copy - EW Div. IATOP Office, Code: SIGRA/SL-LNR 31 President 50 Director U.S. Army Defense Board Fort Monmouth Office Headquarters Communications-Electronics Combat Fort Bliss, Texas Developments Agency Building 410 Fort Monmouth, New Jersey

(U) DISTRIBUTION LIST (Cont.) Copy Copy No. No. 51-65 Commanding Officer 68 U.S.A. F. Project Rand U. S. Army Electronics Command The Rand Corporation U. S. Army Electronics Laboratories 1700 Main Street Fort Monmouth, New Jersey Santa Monica, California ATTN: AMSEL-RD- DR AMSEL-RD-NSR 69 Stanford Electronics Laboratories AMSEL-RD-SM Stanford University AMSEL-RD-SA Stanford, California AMSE L-RD-SEA AMSEL-RD-SEJ 70 Director AMSEL-RD-SES National Security Agency AMSEL-RD-SEE Fort George G. Meade, Maryland AMSEL-RD-ADO ATTN: RADE- 1 AMSEL-RD-SR AMSEL-RD-SE 71 Bureau of Naval Weapons AMSEL-RD-ADT Representative AMSEL-RD-GFR Lockheed Missiles and Space AMSEL-RD- PRM Company AMSEL-RD-RHA P.O. Box 504 Sunnyvale, California 66 Commanding Officer U. S. Army Signal Missile Support 72 Dr. T. W. Butler, Jr., Director Agency Cooley Electronics Laboratory White Sands Missile Range The University of Michigan White Sands, New Mexico Ann Arbor, Michigan ATTN: SIGWS- MEW 73-83 Cooley Electronics Laboratory 67 Commanding Officer The University of Michigan U. S. Naval Air Development Center Ann Arbor, Michigan Johnsville, Pennsylvania ATTN: Naval Air Development Center Library Above distribution is effected by Electronic Warfare Division, Surveillance Department, USAEL, Evans Area, Belmar, New Jersey. For further information contact Mr. I. O. Myers, Senior Scientist, Telephone 59-61252.

UNCLASSIFIED Security Classification DOCUMENT CONTROL DATA - R&D (Security claesification of title, body of abstract and indexing annotation must be entered when the overall report is classified) " ORIGINATIN G ACTIVITY (Corporate author) 2a. REPORT SECURITY C LASSIFICATION The University of Michigan UNCLASSIFIED Cooley Electronics Laboratory 2b. GROUP Ann Arbor, Michigan REPORT TITLE The Recovery of Clocked Signal Inputs to Linear Dynamic Systems DESCRIPTIVE NOTES (Type of report and inclusive dates) Technical Report No. 166. AUTHOR(S) (Last name, first name, initial) Bailey, F. N. and Damborg, M. J. REPO RT DATE 7.. TOTAL NO. OF PAGES 7b. NO. OF REF. September 1965 1_ a. CONTRACT OR GRANT NO. 9a. ORIGINATOR'S REPORT NUMBER(S) DA 36-039-AMlC-03733(E) 6137-l1-T b. PROJECT NO. 1PO 21101 Ao42 02 c. 9b. OTHER REPORT NO(S) (Any other numbers that may be assigned this report) d. A VaiL ABILITY/LTATI ed reuesters mayON NOTICpies of thiES Availability: Qualified requesters may obtain copies of this report from DDC. Foreign announcement and dissemination by DDC is limited. I. SUPPL EMENTARY NOTES 12. SPONSORING MILITARY ACTIVITY U. S. Army Electronics Command.... __~~~Attn: AMSEL-WL-S Ft. Monmouth, N. Jo I. ABSTRACT The primary topic treated in this report is that of the recovery of a clocked input to a known linear time-invariant system when the output, or at least a noisy version of the output, can be observed. A secondary topic, receiving minimal attention here, is that of identifying the time-invariant linear system when the output can be observed and certain a priori statistical information about the clocked input is available (e.g. the power spectrum). Both these problems are of interest in data transmission and communication systems using clocked inputs. The input recovery problem is treated separately for the two cases of periodic inputs and aperiodic inputs. In both cases a transform (a type of "inverse filter" different from the standard inverse filter) is developed which operates on the output to recover the input. This transform can be realized in the form of a computer program. D 1JAN 64 1473 UNCLASSIFIED Security Classification

UNC LASS IF IED Security Classification 14. LINK A LINK B LINK C KEY WORDS ROLE WT ROLE WT ROLE WT C ocked sequences Linear system INSTRUCTIONS 1. ORIGINATING ACTIVITY: Enter the name and address imposed by security classification, using standard statements of the contractor, subcontractor, grantee, Department of De- such as: fense activity or other organization (corporate author) issuing (1) "Qualified requesters may obtain copies of this the report. report from DDC." 2a. REPORT SECURITY CLASSIFICATION: Enter the over- (2) "Foreign announcement and dissemination of this all security classification of the report. Indicate whether report by DDC is not authorized" "Restricted Data" is included. Marking is to be in accord-s ance with appropriate security regulations. (3) "U. S. Government agencies may obtain copies of this report directly from DDC. Other qualified DDC 2b. GROUP: Automatic downgrading is specified in DoD Di- users shall request through rective 5200. 10 and Armed Forces Industrial Manual. Enter the group number. Also, when applicable, show that optional.. markings have been used for Group 3 and Group 4 as author- (4) "U. S. military agencies may obtain copies of this ized. report directly from DDC Other qualified users 3. REPORT TITLE: Enter the complete report title in all shall request through capital letters, Titles in all cases should be unclassified., If a meaningful title cannot be selected without classifica- tion, show title classification in all capitals in parenthesis (5) "All distribution of this report is controlled. Qualimmediately following the title. ified DDC users shall request through 4. DESCRIPTIVE NOTES: If appropriate, enter the type of - b report, e.g., interim, progress, summary, annual, or final. If the report has been furnished to the Office of Technical Give the inclusive dates when a specific reporting period is Services, Department of Commerce, for sale to the public, indi. covered. cate this fact and enter the price, if known. 5. AUTHOR(S): Enter the name(s) of author(s) as shown on 11. SUPPLEMENTARY NOTES: Use for additional explanaor in the report. Entel last name, first name, middle initial. tory notes. If military, show rank and branch of service. The name of the principal au..thor is an absolute minimum requirement. 12, SPONSORING MILITARY ACTIVITY: Enter the name of the departmental project office or laboratory sponsoring (pay6. REPORT DAT: Enter the date of the report as day, ing for) the research and development. Include address. month, year; or month, year. If more than one date appears on the repor month, use date of publication 13. ABSTRACT: Enter an abstract giving a brief and factual summary of the document indicative of the report, even though 7a. TOTAL NUMBER OF PAGES:.The total page count it may also appear elsewhere in the body of the technical reshould follow normal pagination procedures, i.e., enter the port. If additional space is required, a continuation sheet shae number of pages containing information. be attached. 7b. NUMBER OF REFERENCES: Enter the total number of It is highly desirable that the abstract of classified report references cited in the report. be unclassified. Each paragraph of the abstract shall end witt 8a. CONTRACT OR GRANT NUMBER: If appropriate, enter an indication of the military security classification of the inthe applicable number of the contract or grant under which formation in the paragraph, represented as (TS), (S), (C), or (U. the report was written. There is no limitation on the length of the abstract. How8b, 8c, & 8d. PROJECT NUMBER: Enter the appropriate ever, the suggested length is from 150 to 225 words. military department identification, such as project number, subproject number system numbers, task number, etc. 14. KEY WORDS: Key words are technically meaningful term, subprject-umber sysemnuor short phrases that characterize a report and may be used as 9a. ORIGINATOR'S REPORT NUMBER(S): Enter the offi- index entries for cataloging the report. Key words must be cial report number by which the document will be identified selected so that no security classification is required. Identiand controlled by the originating activity. This number must fiers, such as equipment model designation, trade name, milita be unique to this report. project code name, geographic location, may be used as key 9b. OTHER REPORT NUMBER(S): If the report has been words but will be followed by an indication of technical conassigned any other report numbers (either by the originator text. The assignment of links, rules, and weights is optional. or by the sponsor), also enter this number(s), 10. AVAILABILITY/LIMITATION NOTICES: Enter any limitations on further dissemination of the report, other than those UNCLASSIFIED Security Classification

UNIVERSITY O.F MICHIGAN 3 9015 02493 89801111111ll111111I111111l 3 9015 02493 8980