13.3 Optimal(Wiener)Filtering with the FFT 547 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 from a"corrupted"signal.The particular situation we consider is this:There is some underlying,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, so that the true signal u(t)is convolved with(smeared out by)some known response function r(t)to give a smeared signal s(t), s(t)= r(t-T)u(T)dr or S(f)=R(f)U(f) (13.3.1) 蓉 where S,R,U are the Fourier transforms of s,r,u,respectively.Second,the measured signal c(t)may contain an additional component of noise n(t), RECIPES 9 c(t)=s(t)+n(t) (13.3.2) Press. We already know how to deconvolve the effects of the response function r in the absence of any noise($13.1);we just divide C(f)by R(f)to get a deconvolved signal.We now want to treat the analogous problem when noise is present.Our task is to find the optimal filter,o(t)or (f),which,when applied to the measured 色g台 signal c(t)or C(f),and then deconvolved by r(t)or R(f),produces a signal u(t) OF SCIENTIFIC( or U(f)that is as close as possible to the uncorrupted signal u(t)or U(f).In other words we will estimate the true signal U by 6 0U)=Cf)Φf) (13.3.3) R(f) In what sense is U to be close to U?We ask that they be close in the 10.621 least-square sense Numerica 431 )-vp a=-ud is minimized. (13.3.4) (outside Recipes Substituting equations(13.3.3)and (13.3.2),the right-hand side of (13.3.4)becomes North [S(f)+N(f)(f)s(f) R(f) R(f) d (13.3.5) IRf12{Sf)PI1-Φ(f)P+1NfI21()P}d西 The signal S and the noise N are uncorrelated,so their cross product,when integrated over frequency f,gave zero.(This is practically the definition of what we 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 a13.3 Optimal (Wiener) Filtering with the FFT 547 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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 from a “corrupted” signal. The particular situation we consider is this: There is some underlying, 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, so that the true signal u(t) is convolved with (smeared out by) some known response function r(t) to give a smeared signal s(t), 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, the measured 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 in the absence of any noise (§13.1); we just divide C(f) by R(f) to get a 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) or C(f), and then deconvolved by r(t) or R(f), produces a signal u(t) or U(f) that is as close as possible to the uncorrupted signal u(t) or U(f). In other words we will estimate the true signal U by U(f) = C(f)Φ(f) R(f) (13.3.3) In what sense is U to be close to U? We ask that they be close in the least-square sense ∞ −∞ |u(t) − u(t)| 2 dt = ∞ −∞    U(f) − U(f)    2 df is minimized. (13.3.4) Substituting equations (13.3.3) and (13.3.2), the right-hand side of (13.3.4) becomes ∞ −∞     [S(f) + N(f)]Φ(f) R(f) − S(f) R(f)     2 df = ∞ −∞ |R(f)| −2  |S(f)| 2 |1 − Φ(f)| 2 + |N(f)| 2 |Φ(f)| 2  df (13.3.5) The signal S and the noise N are uncorrelated, so their cross product, when integrated over frequency f, gave zero. (This is practically the definition of what we 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