Data adaptive signal estimation by singular value decomposition of a data matrix

PROCEEDINGS OF THE IEEE, VOL. 70, NO. 6, JUNE 1982

684

K. Iizuka and P. Kmtky, “Step frequency radar applied to a fiber ODtiC fault locator.” submitted to hoc.IEEE. Handbook of Chemistry and Phy&. Cleveland,OH:Chemical Rubber Pub. Co., 1961.

Data Adaptive Signal Estimation by Singular Value Decomposition of a Data Matrix DONALD W. TUFTS, RAhiDAS KUMARESAN, AND WARS KIRSTEINS

Akamct-A new method is pmented for estimating the s i g d component of a noisy record of data Only a little prior information about the &M is asmned. I S p e d i d l y , the appmximate value of rank of a matrbrwhichisformedfmmthesrmplesofthesignnlisrsannedtobe known or obtainable fmn singulnvdue decomposition (SVD).

I. INTRODUCTION In 1936, Eclrart and Young presented the derivation of a procedure for fmding the best lower rank approximation to a given matrix [ 11. We believe that this procedure can be applietl to a wide variety ofproblems in signal estimation,estimation of parametersof signals, and approximation theory. The starting point for our discussion is the singular valuedecomposition (SVD) of a rectangular matrixA which has real or complex entries. Eckart and Young [2] showed that the SVD of such a matrix A can be used to fmd an approximant to A oflower rank. The SVD of the matrix A can be specified by the following product of three matrices

131: A

=

U

I:

V*.

(1)

m X nm X m m X n n X n (C)

Fa.2.

Scope display of acoustic timedomain reflectometer (a) with a d m ITT graded-index optical fiber; (b) when the fiber glass is broken at 4.5 m and the broken part is slightly bent; (c) like (b) but thefiber was realigned to be smoothly straight.

The attenuation of the cable also varies depending on how the optical fiber is armored. Any restrictive forcecaused by the armor significantly influences the signal throughput. For instance,attenuationofthe armored Belden 220001 cable is about twice as high as those previously mentioned.

Thedimensions ofeach matrix are writtenbelow the matrix.The matrices U and V are unitary, and I: is a rectangular diagonal matrix of the same size as A with real, nonnegative diagonal entries. These diagonalentries, called the singularvalues of A , areconventionally ordered in decreasing orderwith the largest intheupperlefthand corner. These singuh values are the nonnegative square roots of the eigenvalues of A*A and AA*. The asterisk is used to denote the complex conjugate transpose of a matrix. The pseudoinverse of A , denoted by AI, is related to the S V D of A by the formula

Az = V d U *

(2)

CONCLUSION The ATDR was tried out as a fiber-optic fault locator. The reflected power from the break has very little to do with the condition of the fractured surface and certainly has no such limitation of a surface tilt angle of 6.5’. The disadvantage, however,is a limited range of detection due to high attenuation of the acoustic wave in the fiber. The ATDRwould be used as a short-rangefaultlocatorandcouldsupplement the other methods whichhavea limitation on the tilt angle of the fractured surface.

in which ZZ is obtained from I: by replacing each positive diagonal entry by itsreciprocal. The theorem ofEckart and Young [ 11 can be stated as follows: Eckmt-Young 7Reorem: Let A be an m X n matrix of rank r which has complex elements. Let Sp be the set of all m X n matrices of rank p < r . ThenforallmatricesBinSp

Acwow~e~om The author wishes to thank B. Maillard of the University of Toronto for helping him to conduct the experiments andreading themanuscript

and 2 is obtained from the matrix L of (1) by setting to zero all but its p largest singular The matrix norm of (3) is the Frobenius

REFERENCES Y. Ueno and M. Shimizu, “Optical fiber fault location method,” AppL O p t , v o l . 15, no. 6,pp. 1385-1388,1976. S. D. Personick, “Photo probe--An optical fiber timedomain reflectometer,” Ben Syst. Tech. /., vol. 56, no. 3, pp. 355-366. 1977. M. K. Barnoski, M. D. Rourke. S. M. Jemen. and R. T. Melville, “Optical time domain reflectometer,” AppL Opt, vol. 16, pp. 2375-2319.1978. K. Aoyama, K. Nakapwa, and T. Itoh, “Optical time domain reflectometer in a single mode fiber,” IEEE Tram Quantum Electron,vol. QE17, no: 6,pp. 862468,1981. R. I. MacDonald, “Frequency domain optical reflectometer,’’ AppL Opt.,vol. 20,no. 10,pp. 1840-1844,1981.

IIA -211< IIA -Ell

(3)

where

2=USv* norm. That is

values.

IIA -BIl2 = tr [(A - B ) * (A - E ) ] .

(5)

~ e n c ein words, 2 is the best least squares approximation of lower rank p to thegiven matrix. Because of the intimate connection between least squares approximationand maximum likelihood and maximum posterior probability signal estimation in Gaussiannoise, the applicability of the

=art-

Support for this work was Manuscript received December 18, 1981. providedby a grant from the Office of Naval Researchunder the Probability and Statistics Program. The authors are with the Department of ElectricalEngineering, Kelley Hall, University of Rhode Island, Kingston, RI 02881.

0018-9219/82/0600-0684$00.75 8 1982 IEEE

PROCEEDINGS OF THE IEEE, VOL. 7 0 , NO. 6 , JUNE 1982

685

Young theorem should not be surprising. To clarify this applicability we provide two examples. 11. EXAMPLES

Example I Considera two-dimensional signal which consists of rn uniformly spaced (in time) complex valuedsamples from each of n uniformly (or nonuniformly) spaced(along a line in space) sensor elements. A

monochromatic, but not necessarily plane, wave is incident on the line array. The signal is then given by the formula

S(k,I ) = c exp (twk) exp (ie ( I ) ) ,

i = &i

SGUAL

(6)

inwhichk=1,2,...,mandI=1,2,.--,n . The(mXn)matrixof signal samples has rank 1 because each row can be expressed as a complex scale factor times the f m t row. Therefore, if the A matrix consists white-noise samples, it of such a signal matrix plus a matrix of complex is reasonable to estimate the signal samples by least s q u a r e s approximation of the data matrix by a matrix of rank 1. If we had reason to believe thatthe signal consistedof thesuperpositionof three such and wavefronts shapes, then a components with distinct frequencies rank-3approximation would be appropriateforextractionof threecomponent signaL Example 2 Consider asnapshotof given by theformula

the

25 uniformlyspaced, real-valuedsamples

y(k) = c o s (2n(0.15) k) + c o s (2n(0.30) k) + w(k)

(7)

--

for k = 1 , 2 , . , 25 and in which w(k) is a sequence of zero-mean, mutually independent, Gaussian random variables, each of which has a variance o2 = 0 5 . (SNR is 0 dB for each sine wave.) Motivated by the results of forward/backward linear prediction [4], [5], we form the following data matrix:

? Y ( W ~ ( 1 7 )* ~ ( 1 9 ) ~ ( 1 8 )*

A=

* *

- *

but the fust p columns of U and V can be set equal to zero to form U and V p . One is then led to the following three additional forms of

(4:

h

-

A = UpC^V$= U p U $ A = A V p V $ .

~(1) ~(2)

~ b 5 ) ~ ( 2 4* *) * vi81 ~ ( 1 ) ~ ( 2 ) *.* ~(18) Y(2) Y ( 3 ) * * * Y(19)

-v&)

Fig. 1. The 25 samples of signal plus noise, signal, and signal estimate waveforms are shown. SNR = 0 dB. The 25 discrete samplevalues are joined smoothlyby interpolating in between the samples.

(8)

Because the columns of U and V are eigenvectors of A A * and A * A , respectively, the smaller size of one or the other of the latter square matrices can lead one to avoid the SVD by 1) computing the principal eJgenvectorsof A A * or A *A whichever is easier and 2) computing A from either of the last two expressions in(9). Conceptually, the last two expressions of (9; are useful in providing an adaptive Wiener ffltering interpretation of A . Note that UpUjand VpV$ are matrix filters which operate directlyon the data m a wA . ReFERENCEs

~ b )

*

~b5)>

If there were nonoise, the rank of theabove matrix would be 4. Hence, assuming that we do not know the frequencies of the sinusoidal components, we estimate the signal portion of the matrixA of (8) by fmding the least squares approximationto A which has rank 4. To estimate the original 25 samples of signal we arithmetically average the multiple occu~ancesof each such samplewhich arise because of the redundancies in the matrix A of (8). An example of the waveforms of signal plus noise, signal, and estimated signal are shown in Fig. 1. We have performed experiments on the estimation of one real sinusoidal component using data like that of (7) with one sinusoid and different SNR values. The result is that there is about 6 d B improve ment in SNR of the estimated signal with respect to that obtained by using an adaptive linear prediction filter with an impulse response of length 8 samples, the coefficients of which are obtained by minimizing the sum of squared prediction erzors aver (25 8) = 17 samples of the 25 samplesof observed data. The Eckart-Young theoremcanbeused toimprove theill-coditioned nature of signal parameter estimation via linear prediction by replacing the raw data matrixbya “cleaned-up” data matrixof lower rank [61-[101* To applythe Mart-Young theorem it is not necessary that the signal be deterministic nor that the signal matrix be of less than full rank. What is important is that the signal matrix be of approximately lower rank in the sense that, with high probability, the signahnly matrix can be well approximated by a matrix of lower rank. This will happen for a random signal when its estimated correlationmatrix has a few eigenvalues which are si@fhntly g g e r than the other [8], [ 111. W~threspect to the computationof A of (4), we hrst note that there are many unnecesary, implied computations in (4). For example, all

-

[ 11 C. Eckart and G. Young, “The approximation of one matrix by another of lower rank,” Psychometrfku, vol. 1, pp. 211-218, 1936. [ 21 -, “A principal axis transformationfornon-Hermetian matrices,”BuU. Amer. Math. Soc.,vol. 45,pp. 118-121, 1939. [ 31 C. L. L a w n and R. J. Hanson, Solving Least Square8Roblem.s. Englewood Cliffs, NJ: Rentice-Hall, 1974. [ 4 ] T. J. Ulrych and R. W. Clayton, “Time series modelling and maximum entropy,” Physics Earth Planetary Interiom, vol. 12, pp. 188-200, Aug. 1976. [ 5 ) A. H. Nuttall, “Spectral analysis of a univariate process with bad

datapoints via maximum entropyand linearpredictive techniques,” in Specmrl Estimution (NUSC scientific and engineering studies). New London, CT: NUSC, Mar. 1976. [ 61 D. W. Tufts and R. Kumaresan, “Frequency estimation of multiple sinusoids: Making linearpredictionperform likemaximum likelihood,” submitted for publication to IEEE Trans. Acoust., Speech, Sllglllll Roces. [ 7 ] -, “Singularvalue decomposition andimproved frequency estimation using linear prediction,” to be published in IEEE Trcvu Acoust., Speech,SignaJ R o c e s . , vol. ASSP-30, Aug. 1982. [ 8 ] R.Kumaresanand D. W. Tufts, “Singularvalue decomposition and spectral analysis,” in 1st ASSP Workshop on Spectml Estimation, Hamilton, Ont., (Aug. 17 and 18, 1981) pp. 6.4.1-

Roc.

6.4.12. [ 9 ] -, “Estimation of parameters of exponentially damped signals in noise and pole-zero modeling,” submitted for publication to the IEEE Tmnr Acoust., Speech,Sllglllll Rocem. [ 101 -, “ImprovedF’isarenko-type methods in spectral analysis,’’ submitted for publication to the IEEE Tmnr Acoust., Speech, SignalProcess. [ 11 ] -, “Data adaptiveprincipalcomponent signal processing,” in IEEE Con$ on Denldon and Control (Albuquerque, NM), Dec. 1980.

Roc.