c13-3 (Numerical Recipes in C)

54713.3 Optimal (Wiener) Filtering with the FFT13.3 Optimal (Wiener) Filtering with the FFTZ∞s(t) =−∞r(t − τ )u(τ ) dτor S(f) = R(f)U (f)(13.3.1)where S, R, U are the Fourier transforms of s, r, u, respectively. Second, themeasured signal c(t) may contain an additional component of noise n(t),c(t) = s(t) + n(t)(13.3.2)We already know how to deconvolve the effects of the response function r inthe absence of any noise (§13.1); we just divide C(f) by R(f) to get a deconvolvedsignal.

We now want to treat the analogous problem when noise is present. Ourtask is to find the optimal filter, φ(t) or Φ(f), which, when applied to the measuredsignal c(t) or C(f), and then deconvolved by r(t) or R(f), produces a signal ue(t)eor U(f) that is as close as possible to the uncorrupted signal u(t) or U (f). In otherwords we will estimate the true signal U bye (f) = C(f)Φ(f)UR(f)e to be close to U ?In what sense is Uleast-square senseZ∞−∞Z2|eu(t) − u(t)| dt =∞−∞(13.3.3)We ask that they be close in the2eU(f) − U (f) dfis minimized.(13.3.4)Substituting equations (13.3.3) and (13.3.2), the right-hand side of (13.3.4) becomes2 [S(f) + N (f)]Φ(f)S(f) −dfR(f)R(f) −∞Z ∞no−22222|S(f)| |1 − Φ(f)| + |N (f)| |Φ(f)|df|R(f)|=Z∞(13.3.5)−∞The signal S and the noise N are uncorrelated, so their cross product, whenintegrated over frequency f, gave zero.

(This is practically the definition of what wemean by noise!) Obviously (13.3.5) will be a minimum if and only if the integrandSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).There are a number of other tasks in numerical processing that are routinelyhandled with Fourier techniques.

One of these is filtering for the removal of noisefrom a “corrupted” signal. The particular situation we consider is this: There is someunderlying, uncorrupted signal u(t) that we want to measure. The measurementprocess is imperfect, however, and what comes out of our measurement device is acorrupted signal c(t).

The signal c(t) may be less than perfect in either or both oftwo respects. First, the apparatus may not have a perfect “delta-function” response,so that the true signal u(t) is convolved with (smeared out by) some known responsefunction r(t) to give a smeared signal s(t),548Chapter 13.Fourier and Spectral Applicationsis minimized with respect to Φ(f) at every value of f. Let us search for such asolution where Φ(f) is a real function.

Differentiating with respect to Φ, and settingthe result equal to zero givesΦ(f) =|S(f)|222(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 notcontain U , the “true” signal.

This makes for an important simplification: The optimalfilter can be determined independently of the determination of the deconvolutionfunction that relates S and U .To determine the optimal filter from equation (13.3.6) we need some way22of separately estimating |S| and |N | . There is no way to do this from themeasured signal C alone without some other information, or some assumption orguess. Luckily, the extra information is often easy to obtain. For example, wecan sample a long stretch of data c(t) and plot its power spectral density usingequations (12.0.14), (12.1.8), and (12.1.5). This quantity is proportional to the sum22|S| + |N | , so we have22|S(f)| + |N (f)| ≈ Pc(f) = |C(f)|20 ≤ f < fc(13.3.7)(More sophisticated methods of estimating the power spectral density will bediscussed in §13.4 and §13.7, but the estimation above is almost always good enoughfor 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 continuousnoise spectrum.

The noise spectrum may be flat, or tilted, or smoothly varying; itdoesn’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 regiondominated by the signal as well. Now draw a smooth curve through the signal plusnoise power. The difference between these two curves is your smooth “model” ofthe signal power. The quotient of your model of signal power to your model ofsignal plus noise power is the optimal filter Φ(f).

[Extend it to negative values of fby the formula Φ(−f) = Φ(f).] Notice that Φ(f) will be close to unity where thenoise is negligible, and close to zero where the noise is dominant. That is how itdoes its job! The intermediate dependence given by equation (13.3.6) just turns outto be the optimal way of going in between these two extremes.Because the optimal filter results from a minimization problem, the quality ofthe results obtained by optimal filtering differs from the true optimum by an amountthat is second order in the precision to which the optimal filter is determined. In otherwords, even a fairly crudely determined optimal filter (sloppy, say, at the 10 percentlevel) can give excellent results when it is applied to data.

That is why the separationof the measured signal C into signal and noise components S and N can usefully bedone “by eye” from a crude plot of power spectral density. All of this may give youthoughts about iterating the procedure we have just described. For example, afterdesigning a filter with response Φ(f) and using it to make a respectable guess at thee (f) = Φ(f)C(f)/R(f), you might turn about and regard Ue (f) as a freshsignal USample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).|S(f)| + |N (f)|54913.4 Power Spectrum Estimation Using the FFT N 2 (extrapolated) S 2 (deduced)fFigure 13.3.1. Optimal (Wiener) filtering.

The power spectrum of signal plus noise shows a signal peakadded to a noise tail. The tail is extrapolated back into the signal region as a “noise model.” Subtractinggives the “signal model.” The models need not be accurate for the method to be useful. A simplealgebraic combination of the models gives the optimal filter (see text).new signal which you could improve even further with the same filtering technique.Don’t waste your time on this line of thought. The scheme converges to a signal ofS(f) = 0. Converging iterative 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 datawhen you are constructing an optimal filter.

To apply the filter to your data, youcan use the methods described in §13.1. The specific routine convlv is not neededfor optimal filtering, since your filter is constructed in the frequency domain tobegin with. If you are also deconvolving your data with a known response function,however, you can modify convlv to multiply by your optimal filter just before ittakes the inverse Fourier transform.CITED REFERENCES AND FURTHER READING:Rabiner, L.R., and Gold, B.

1975, Theory and Application of Digital Signal Processing (EnglewoodCliffs, NJ: Prentice-Hall).Nussbaumer, H.J. 1982, Fast Fourier Transform and Convolution Algorithms (New York: SpringerVerlag).Elliott, D.F., and Rao, K.R. 1982, Fast Transforms: Algorithms, Analyses, Applications (NewYork: Academic Press).13.4 Power Spectrum Estimation Using the FFTIn the previous section we “informally” estimated the power spectral density of afunction c(t) by taking the modulus-squared of the discrete Fourier transform of someSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).log scale C 2 (measured).