f13-3

13.3Optimal(Wiener)FilteringwiththeFFT 539Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).input data sets. If youwant all possible lags from Ndata points (nota usual thing), thenyouwillneedtopadthedatawithanequalnumberofzeros;thisis theextremecase. So here is the program: SUBROUTINE correl(data1,data2,n,ans) INTEGER n,NMAX REAL data1(n),data2(n)COMPLEX ans(n) PARAMETER (NMAX=4096) Maximum anticipated FFT size. C USES realft,twofft Computes the correlation of two real data sets data1(1:n) anddata2(1:n) (includ- ing any user-supplied zero padding). nMUST be an integer power of two. The answer is returned as the first npoints in ans stored in wrap-around order, i.e., correlations at increasingly negative lags are in ans(n) on down to ans(n/2+1) , while correlations at increasingly positive lags are in ans(1) ( z e r ol a g )o nu pt o ans(n/2) .N o t et h a t ans must be supplied in the calling program with length at least 2*n , since it is also used as working space. Sign convention of this routine: if data1 lagsdata2 , i.e., is shifted to the right of it, then ans will show a peak at positive lags. INTEGER i,no2 COMPLEX fft(NMAX) call twofft(data1,data2,fft,ans,n) Transform b oth data vectors at once. no2=n/2 Normalization for inverse FFT. do11i=1,no2+1 ans(i)=fft(i)*conjg(ans(i))/float(no2) Multiply to find FFT of their corre- lation. enddo 11 ans(1)=cmplx(real(ans(1)),real(ans(no2+1))) Pack first and last into one element. call realft(ans,n,-1) Inverse transform gives correlation. returnEND As in convlv, it would be better to substitute two calls to realftfor the one call to twofft,i fdata1anddata2have very different magnitudes, to minimize roundoff error. Thediscrete autocorrelation of a sampled function gjis just the discrete correlation of the function with itself. Obviously this is always symmetric withrespect to positive and negative lags. Feel free to use the above routine correl to obtain autocorrelations, simply calling it with the same datavector in both arguments. If the inefficiency bothers you, routine realftcan, of course, be used to transform the datavector instead. CITED REFERENCES AND FURTHER READING: Brigham, E.O. 1974, TheFast FourierTransform (EnglewoodCliffs, NJ: Prentice-Hall), §13–2. 13.3 Optimal (Wiener) Filtering with the FFT There are a number of other tasks in numerical processing that are routinely handled with Fourier techniques. One of these is filtering for the removal of noise froma“corrupted”signal. Theparticularsituationweconsideristhis: Thereissomeunderlying, uncorrupted signal u(t)that we want to measure. The measurement process is imperfect, however, and what comes out of our measurement device is a corrupted signal c(t). The signal c(t)may be less than perfect in either or both of two respects. First, the apparatus may not have a perfect “delta-function”response, 540 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).sothatthetruesignal u(t)isconvolvedwith(smearedoutby)someknownresponse function r(t)to give a smeared signal s(t), s(t)=/integraldisplay∞ −∞r(t−τ)u(τ)dτorS(f)=R(f)U(f)( 13.3.1 ) where S, R, Uare the Fourier transforms of s, r, u,respectively. Second, the measured signal c(t)may contain an additional componentof noise n(t), c(t)=s(t)+n(t)( 13.3.2 ) We already know how to deconvolve the effects of the response function rin the absenceofanynoise( §13.1);we just divide C(f)byR(f)to geta deconvolved signal. We now want to treat the analogous problem when noise is present. Our task is to find the optimal filter ,φ(t)orΦ(f), which, when applied to the measured signal c(t)orC(f), and then deconvolved by r(t)orR(f), produces a signal /tildewideu(t) or/tildewideU(f)that is as close as possible to the uncorruptedsignal u(t)orU(f). In other words we will estimate the true signal Uby /tildewideU(f)=C(f)Φ(f) R(f)(13.3.3 ) In what sense is /tildewideUto be close to U? We ask that they be close in the least-square sense /integraldisplay∞ −∞|/tildewideu(t)−u(t)|2dt=/integraldisplay∞ −∞/vextendsingle/vextendsingle/vextendsingle/tildewideU(f)−U(f)/vextendsingle/vextendsingle/vextendsingle2 dfis minimized. (13.3.4 ) Substitutingequations(13.3.3)and(13.3.2),theright-handsideof(13.3.4)becomes /integraldisplay∞ −∞/vextendsingle/vextendsingle/vextendsingle/vextendsingle[S(f)+N(f)]Φ(f) R(f)−S(f) R(f)/vextendsingle/vextendsingle/vextendsingle/vextendsingle2 df =/integraldisplay∞ −∞|R(f)|−2/braceleftBig |S(f)|2|1−Φ(f)|2+|N(f)|2|Φ(f)|2/bracerightBig df(13.3.5 ) The signal Sand the noise Nareuncorrelated , so their cross product, when integratedoverfrequency f,gavezero. (Thisis practicallythe definition ofwhatwe mean by noise!) Obviously (13.3.5)will be a minimum if and only if the integrand is minimized with respect to Φ(f)at every value of f. Let us search for such a solutionwhere Φ(f)isarealfunction. Differentiatingwithrespectto Φ,andsetting the result equal to zero gives Φ(f)=|S(f)|2 |S(f)|2+|N(f)|2(13.3.6 ) This is the formula for the optimal filter Φ(f). Notice that equation(13.3.6)involves S, the smeared signal, and N, the noise. The two of these add up to be C, the measured signal. Equation (13.3.6) does not 13.3Optimal(Wiener)FilteringwiththeFFT 541Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).contain U,the“true”signal. Thismakesforanimportantsimplification: Theoptimal filter can be determined independently of the determination of the deconvolutionfunction that relates SandU. To determine the optimal filter from equation (13.3.6) we need some way of separately estimating |S| 2and|N|2. There is no way to do this from the measured signal Calone without some other information, or some assumption or guess. Luckily,theextrainformationisofteneasytoobtain. Forexample,wecansamplea longstretchofdata c(t)andplotitspowerspectraldensityusingequations(12.0.14), (12.1.8),and(12.1.5). Thisquantityisproportionaltothesum |S|2+|N|2,sowehave |S(f)|2+|N(f)|2≈Pc(f)=|C(f)|20≤f<f c (13.3.7 ) (More sophisticated methods of estimating the power spectral density will be discussedin §13.4and §13.7,buttheestimationaboveisalmostalwaysgoodenough for the optimal filter problem.) The resulting plot (see Figure 13.3.1) will oftenimmediately show the spectral signature of a signal sticking up above a continuous noise spectrum. The noise spectrum may be flat, or tilted, or smoothly varying; it doesn’t matter, as long as we can guess a reasonable hypothesis as to what it is.Draw a smooth curve through the noise spectrum, extrapolating it into the region dominated by the signal as well. Now draw a smooth curve throughthe signal plus noisepower. Thedifferencebetweenthesetwocurvesisyoursmooth“model”ofthe signal power. The quotient of your model of signal power to your model of signal plusnoise poweris the optimalfilter Φ(f). [Extendit tonegativevaluesof fbythe formula Φ(−f)=Φ ( f).] Notice that Φ(f)will be close to unity where the noise is negligible, and close to zero where the noise is dominant. That is how it does its job! Theintermediatedependencegivenbyequation(13.3.6)justturnsouttobetheoptimal way of going in between these two extremes. Because the optimal filter results from a minimization problem, the quality of theresults obtainedbyoptimalfilteringdiffersfromthetrueoptimumbyanamount thatissecondorder intheprecisiontowhichtheoptimalfilterisdetermined. Inother words, evena fairly crudelydeterminedoptimalfilter (sloppy,say, at the 10 percentlevel)cangiveexcellentresultswhenitisappliedtodata. Thatiswhytheseparation ofthemeasuredsignal Cintosignalandnoisecomponents SandNcanusefullybe done“byeye” froma crudeplot ofpowerspectral density. All of this may giveyouthoughts about iterating the procedure we have just described. For example, after designinga filter withresponse Φ(f)andusingit tomakea respectableguessat the signal /tildewideU(f)=Φ ( f)C(f)/R(f), you might turn about and regard /tildewideU(f)as a fresh new signal which you couldimproveeven furtherwith the same filtering technique. Don’t waste your time on this line of thought. The scheme convergesto a signal of S(f)=0. Convergingiterative methods do exist; this just isn’t one of them. You can use the routine four1(§12.2) or realft(§12.3) to FFT your data when you are constructing an optimal filter. To apply the filter to your data, you can use the methods described in §13.1. The specific routine convlvis not needed for optimal filtering, since your filter is constructed in the frequency domain to beginwith. If youare also deconvolvingyourdatawith a knownresponsefunction, however, you can modify convlvto multiply by your optimal filter just before it takes the inverse Fourier transform. 542 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica). S 2 (deduced) N 2 (extrapolated) C 2 (measured)log scale f Figure 13.3.1. Optimal (Wiener) filtering. The power spectrum of signal plus noise shows a signal peak added to a noise tail. The tail is extrapolated back into the signal region as a “noise model. ”Subtracting givesthe“signalmodel. ”Themodelsneednotbeaccurate forthemethodtobeuseful. Asimplealgebraic combination of the models gives the optimal filter (see text). CITED REFERENCES AND FURTHER READING: Rabiner,L.R.,andGold,B.1975, TheoryandApplicationofDigitalSignalProcessing (Englewood Cliffs, NJ: Prentice-Hall). Nussbaumer,H.J.1982, FastFourierTransformandConvolutionAlgorithms (NewYork:Springer- Verlag). Elliott,D.F.,andRao,K.R.1982, FastTransforms:Algorithms,Analyses,Applications (NewYork: Academic Press). 13.4 Power Spectrum Estimation Usingthe FFT Intheprevioussectionwe “informally ”estimatedthepowerspectraldensityofa function c(t)bytakingthemodulus-squaredofthediscreteFouriertransformofsome finite,sampledstretchofit. Inthissectionwe ’lldoroughlythesamething,butwith considerablygreaterattentionto details. Ourattentionwill uncoversomesurprises. Thefirst detail is power spectrum (also called a power spectral density or PSD) normalization. In general there is somerelation of proportionalitybetween a measure of the squared amplitude of the function and a measure of the amplitude of the PSD. Unfortunately there are several different conventions for describingthe normalization in each domain, and many opportunities for getting wrong the relationshipbetween the two domains. Supposethat our function c(t)is sampled at Npoints to produce values c 0...c N−1, and that these points span a range of time T, that is T=(N−1)∆, where ∆is the sampling interval. Then here are several