f13-5
PDF · 7 pages · 81.0 KB
Open PDF file
Excerpt from the book Numerical Recipes in Fortran 77 (Cambridge University Press), pp. 551 onward, not Phil's own writing. It argues for filtering in the Fourier domain, then covers causal linear filters, the FIR and IIR recursion formula, the filter response H(f), and designing FIR coefficients by truncating and cyclically shifting FFT-derived coefficients. It also includes the end of the preceding section's code and reference list.
AI-written summary; may contain errors.
Extracted text (machine-read; may contain errors)
13.5DigitalFilteringintheTimeDomain 551Sample 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).do19j=1,m
p(j)=p(j)/den Normalize the output.
enddo 19
return
END
CITED REFERENCES AND FURTHER READING:
Oppenheim, A.V., andSchafer, R.W. 1989, Discrete-Time Signal Processing (EnglewoodCliffs,
NJ: Prentice-Hall). [1]
Harris, F.J. 1978, Proceedings of the IEEE , vol. 66, pp. 51–83. [2]
Childers, D.G. (ed.) 1978, Modern Spectrum Analysis (New York: IEEE Press), paper by P.D.
Welch. [3]
Champeney,D.C.1973, FourierTransformsandTheirPhysicalApplications (NewYork:Academic
Press).
Elliott,D.F.,andRao,K.R.1982, FastTransforms:Algorithms,Analyses,Applications (NewYork:
Academic Press).
Bloomfield, P. 1976, Fourier Analysis of Time Series – An Introduction (New York: Wiley).
Rabiner,L.R.,andGold,B.1975, TheoryandApplicationofDigitalSignalProcessing (Englewood
Cliffs, NJ: Prentice-Hall).
13.5 Digital Filtering in the Time Domain
Suppose that you have a signal that you want to filter digitally. For example, perhaps
youwanttoapply high-pass orlow-pass filtering,toeliminatenoiseatloworhighfrequencies
respectively; or perhaps the interesting part of your signal lies only in a certain frequencyband, so that you need a bandpass filter. Or, if your measurements are contaminated by 60
Hz power-lineinterference, you may need a notch filter to remove only a narrow band around
that frequency. This section speaks particularly about the case in which you have chosento do such filtering in the time domain.
Before continuing, we hope you willreconsider thischoice. Remember how convenient
it is to filter in the Fourier domain. You just take your whole data record, FFT it, multiplythe FFT output by a filter function H(f), and then do an inverse FFT to get back a filtered
data set in time domain. Here is some additional background on the Fourier technique thatyou will want to take into account.
•Remember that you must define your filter function H(f)for both positive and
negative frequencies, and that the magnitude of the frequency extremes is alwaysthe Nyquist frequency 1/(2∆), where ∆is the sampling interval. The magnitude
of the smallest nonzero frequencies in the FFT is ±1/(N∆), where Nis the
number of (complex) points in the FFT. The positive and negative frequencies towhich this filter are applied are arranged in wrap-around order.
•If the measured data are real, and you want the filtered output also to be real, then
your arbitrary filterfunction should obey H(−f)=H(f)*. You can arrange this
most easily by picking an Hthat is real and even in f.
•If your chosen H(f)has sharp vertical edges in it, then the impulse response of
your filter (the output arising from a short impulse as input) will have damped“ringing” at frequencies corresponding to these edges. There is nothing wrongwith this, but if you don’t like it, then pick a smoother H(f). To get a first-hand
look attheimpulseresponse ofyour filter,justtaketheinverse FFTofyour H(f).
If you smooth all edges of the filter function over some number kof points, then
the impulse response function of your filter will have a span on the order of afraction 1/kof the whole data record.
552 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).•If your data set is too long to FFT all at once, then break it up into segments of
any convenient size, as long as they are much longer than the impulse responsefunction of the filter. Use zero-padding, if necessary.
•You should probably remove any trend from the data, by subtracting from it a
straightlinethroughthefirstandlastpoints(i.e.,makethefirstandlastpointsequalto zero). If you are segmenting the data, then you can pick overlapping segmentsand use only the middle section of each, comfortably distant from edge effects.
•A digital filter is said to be causalorphysically realizable if its output for a
particular time-step depends only on inputs at that particular time-step or earlier.It is said to be acausalif its output can depend on both earlier and later inputs.
Filteringin the Fourier domain is,in general, acausal, since the data are processed“in a batch,” without regard to time ordering. Don’t let this bother you! Acausalfilters can generally give superior performance (e.g., less dispersion of phases,sharper edges, less asymmetric impulse response functions). People use causalfilters not because they are better, but because some situations just don’t allow
access to out-of-time-order data. Time domain filters can, in principle, be either
causal or acausal, but they are most often used in applications where physicalrealizability is a constraint. For this reason we willrestrict ourselves to the causalcase in what follows.
Ifyouarestill favoringtime-domainfilteringafterallwehavesaid,itisprobablybecause
you have a real-time application, for which you must process a continuous data stream andwish to output filtered values at the same rate as you receive raw data. Otherwise, it maybe that the quantity of data to be processed is so large that you can afford only a very smallnumber of floating operations on each data point and cannot afford even a modest-sized FFT(withanumberoffloatingoperationsperdatapointseveraltimesthelogarithmofthenumberof points in the data set or segment).
Linear Filters
The most general linear filter takes a sequence xkof input points and produces a
sequence ynof output points by the formula
yn=M/summationdisplay
k=0ckxn−k+N/summationdisplay
j=1djyn−j (13.5.1 )
Here the M+1coefficients ckand the Ncoefficients djare fixed and define the filter
response. Thefilter(13.5.1)produceseachnewoutputvaluefromthecurrentand Mprevious
input values, and from its own Nprevious output values. If N=0, so that there is no
secondsumin(13.5.1),thenthefilteriscalled nonrecursive orfiniteimpulseresponse(FIR) .If
N/negationslash=0,thenitiscalled recursive orinfiniteimpulseresponse (IIR) .(Theterm“IIR”connotes
only that such filters are capableof having infinitely long impulse responses, not that their
impulse response is necessarily long in a particular application. Typically the response of anIIR filter will drop off exponentially at late times, rapidly becoming negligible.)
The relation between the c
k’s and dj’s and the filter response function H(f)is
H(f)=M/summationtext
k=0cke−2πik(f∆)
1−N/summationtext
j=1dje−2πij(f∆)(13.5.2 )
where ∆is,as usual, the sampling interval. The Nyquist interval corresponds to f∆between
−1/2and1/2. For FIR filters the denominator of (13.5.2) is just unity.
Equation (13.5.2) tells how to determine H(f)from the c’s and d’s. To design a filter,
though, we need a way of doing the inverse, getting a suitable set of c’s and d’s — as small
a set as possible, to minimize the computational burden — from a desired H(f). Entire
books are devoted to this issue. Like many other “inverse problems,” it has no all-purpose
13.5DigitalFilteringintheTimeDomain 553Sample 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).solution. One clearly has to make compromises, since H(f)is a full continuous function,
while the short list of c’s and d’s represents only a few adjustable parameters. The subject of
digital filter design concerns itself with the various ways of making these compromises. Wecannot hope to give any sort of complete treatment of the subject. We can, however, sketcha couple of basic techniques to get you started. For further details, you will have to consultsome specialized books (see references).
FIR(Nonrecursive) Filters
When the denominator in (13.5.2) is unity, the right-hand side is just a discrete Fourier
transform. Thetransformiseasilyinvertible,givingthedesiredsmallnumberof ckcoefficients
in terms of the same small number of values of H(fi)at some discrete frequencies fi. This
fact, however, is not very useful. The reason is that, for values of ckcomputed in this way,
H(f)will tend to oscillate wildly in between the discrete frequencies where it is pinned
down to specific values.
A better strategy, and one which is the basis of several formal methods in the literature,
is this: Start by pretending that you are willing to have a relatively large number of filtercoefficients, that is, a relatively large value of M. Then H(f)can be fixed to desired values
on a relatively fine mesh, and the Mcoefficients c
k,k=0,...,M −1can be found by
an FFT. Next, truncate (set to zero) most of the ck’s, leaving nonzero only the first, say,
K,(c0,c1,...,c K−1)and last K−1,(cM−K+1,...,c M−1). The last few ck’s are filter
coefficients at negative lag , because of the wrap-around property of the FFT. But we don’t
want coefficients at negative lag. Therefore we cyclically shift the array of ck’s, to bring
everything to positive lag. (This corresponds to introducing a time-delay into the filter.) Dothis by copying the c
k’s into a new array of length Min the following order:
(cM−K+1,...,c M−1,c0,c1,...,c K−1,0,0,..., 0) ( 13.5.3 )
To see if your truncation is acceptable, take the FFT of the array (13.5.3), giving an
approximation to your original H(f). You will generally want to compare the modulus
|H(f)|to your original function, since the time-delay will have introduced complex phases
into the filter response.
If the new filter function is acceptable, then you are done and have a set of 2K−1
filter coefficients. If it is not acceptable, then you can either (i) increase Kand try again,
or (ii) do something fancier to improve the acceptability for the same K. An example of
something fancier is to modify the magnitudes (but not the phases) of the unacceptable H(f)
to bring it more in line with your ideal, and then to FFT to get new ck’s. Once again set
to zero all but the first 2K−1values of these (no need to cyclically shift since you have
preserved the time-delaying phases), then inverse transform to get a new H(f), which will
often be more acceptable. You can iterate this procedure. Note, however, that the procedurewill not converge if your requirements for acceptability are more stringent than your 2K−1
coefficients can handle.
Thekey idea,inotherwords,istoiteratebetweenthespaceofcoefficientsandthespace
of functions H(f), untila Fourierconjugate pair thatsatisfiesthe imposed constraints in both
spacesis found. A more formal technique for this kind of iteration is the Remes Exchange
Algorithm which produces the best Chebyshev approximation to a given desired frequency
response with a fixed number of filter coefficients (cf. §5.13).
IIR(Recursive) Filters
Recursive filters,whose output atagiven time depends both on thecurrent and previous
inputsandonpreviousoutputs,cangenerallyhaveperformancethatissuperiortononrecursivefilters with the same total number of coefficients (or same number of floating operations perinput point). The reason is fairly clear by inspection of (13.5.2): A nonrecursive filter has afrequency response that is a polynomial in the variable 1/z, where
z≡e
2πi(f∆)(13.5.4 )
554 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).By contrast, a recursive filter’s frequency response is a rational function in1/z. The class of
rational functions is especially good at fitting functions with sharp edges or narrow features,and most desired filter functions are in this category.
Nonrecursive filters are always stable. If you turn off the sequence of incoming x
i’s,
then after no more than Msteps the sequence of yj’s produced by (13.5.1) will also turn off.
Recursive filters, feeding as they do on their own output, are not necessarily stable. If thecoefficients d
jare badly chosen, a recursive filter can have exponentially growing, so-called
homogeneous , modes, which become huge even after the input sequence has been turned off.
This is not good. The problem of designing recursive filters, therefore, is not just an inverseproblem; it is an inverse problem with an additional stability constraint.
How do you tell if the filter (13.5.1) is stable for a given set of c
kanddjcoefficients?
Stability depends only on the dj’s. The filter is stable if and only if all Ncomplex roots
of thecharacteristic polynomial equation
zN−N/summationdisplay
j=1djzN−j=0 ( 13.5.5 )
are inside the unit circle, i.e., satisfy
|z|≤1( 13.5.6 )
The various methods for constructing stable recursive filters again form a subject area
for which you will need more specialized books. One very useful technique, however, is thebilineartransformationmethod . Forthistopicwedefineanewvariable wthatreparametrizes
the frequency f,
w≡tan[π(f∆)] = i/parenleftbigg1−e
2πi(f∆)
1+e2πi(f∆)/parenrightbigg
=i/parenleftbigg1−z
1+z/parenrightbigg
(13.5.7 )
Don’tbefooledbythe i’sin(13.5.7). Thisequationmapsrealfrequencies fintorealvaluesof
w. Infact,itmapstheNyquistinterval −1
2<f∆<1
2ontothereal waxis−∞ <w< +∞.
The inverse equation to (13.5.7) is
z=e2πi(f∆)=1+iw
1−iw(13.5.8 )
In reparametrizing f,walso reparametrizes z, of course. Therefore, the condition for
stability (13.5.5)–(13.5.6) can be rephrased in terms of w: If the filter response H(f)is
written as a function of w, then the filter is stable if and only ifthe poles of the filterfunction
(zeros of its denominator) are all in the upper half complex plane,
Im(w)≥0( 13.5.9 )
Theidea of thebilinear transformation method isthat instead ofspecifying your desired
H(f),youspecifyonlyitsdesiredmodulussquare, |H(f)|2=H(f)H(f)* =H(f)H(−f).
Pick this to be approximated by some rational function in w2. Then find all the poles of this
functioninthe wcomplexplane. Everypoleinthelowerhalf-planewillhaveacorresponding
pole in the upper half-plane, by symmetry. The idea is to form a product only of the factorswith good poles, ones in the upper half-plane. This product is your stably realizable H(f).
Nowsubstituteequation(13.5.7)towritethefunctionasarationalfunctionin z,andcompare
with equation (13.5.2) to read off the c’s and d’s.
The procedure becomes clearer when we go through an example. Suppose we want to
design a simple bandpass filter, whose lower cutoff frequency corresponds to a value w=a,
and whose upper cutoff frequency corresponds to a value w=b, with aandbboth positive
numbers. A simple rational function that accomplishes this is
|H(f)|2=/parenleftbiggw2
w2+a2/parenrightbigg/parenleftbiggb2
w2+b2/parenrightbigg
(13.5.10 )
This function does not have a very sharp cutoff, but it is illustrative of the more general
case. To obtain sharper edges, one could take the function (13.5.10) to some positive integer
13.5DigitalFilteringintheTimeDomain 555Sample 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).power, or, equivalently, run the data sequentially through some number of copies of the filter
that we will obtain from (13.5.10).
The poles of (13.5.10) are evidently at w=±iaandw=±ib. Therefore the stably
realizable H(f)is
H(f)=/parenleftbiggw
w−ia/parenrightbigg/parenleftbiggib
w−ib/parenrightbigg
=/parenleftBig
1−z
1+z/parenrightBig
b
/bracketleftBig/parenleftBig
1−z
1+z/parenrightBig
−a/bracketrightBig/bracketleftBig /parenleftBig
1−z
1+z/parenrightBig
−b/bracketrightBig (13.5.11 )
We put the iin the numerator of the second factor in order to end up with real-valued
coefficients. Ifwe multiply out all the denominators, (13.5.11) can be rewrittenin the form
H(f)=−b
(1+a)(1+b)+b
(1+a)(1+b)z−2
1−(1+a)(1−b)+(1−a)(1+b)
(1+a)(1+b)z−1+(1−a)(1−b)
(1+a)(1+b)z−2(13.5.12 )
from which one reads off the filter coefficients for equation (13.5.1),
c0=−b
(1 + a)(1 + b)
c1=0
c2=b
(1 + a)(1 + b)
d1=(1 + a)(1−b)+( 1 −a)(1 + b)
(1 + a)(1 + b)
d2=−(1−a)(1−b)
(1 + a)(1 + b)(13.5.13 )
This completes the design of the bandpass filter.
Sometimesyoucanfigureouthowtoconstructdirectlyarationalfunctionin wforH(f),
ratherthanhavingtostartwithitsmodulussquare. Thefunctionthatyouconstructhastohaveitspolesonlyintheupperhalf-plane,forstability. Itshouldalsohavethepropertyofgoingintoitsowncomplexconjugateifyousubstitute −wforw,sothatthefiltercoefficientswillbereal.
For example, here is a function for a notch filter, designed to remove only a narrow
frequency band around some fiducial frequency w=w
0, where w0is a positive number,
H(f)=/parenleftbiggw−w0
w−w0−i/epsilon1w0/parenrightbigg/parenleftbiggw+w0
w+w0−i/epsilon1w0/parenrightbigg
=w2−w2
0
(w−i/epsilon1w0)2−w2
0(13.5.14 )
In(13.5.14)theparameter /epsilon1isasmallpositivenumberthatisthedesiredwidthofthenotch,asa
fractionof w0. Goingthroughthearithmeticofsubstituting zforwgivesthefiltercoefficients
c0=1+w2
0
(1 + /epsilon1w0)2+w2
0
c1=−21−w2
0
(1 + /epsilon1w0)2+w2
0
c2=1+w2
0
(1 + /epsilon1w0)2+w2
0
d1=21−/epsilon12w2
0−w2
0
(1 + /epsilon1w0)2+w2
0
d2=−(1−/epsilon1w0)2+w2
0
(1 + /epsilon1w0)2+w2
0(13.5.15 )
556 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).(a)
(b)
Figure 13.5.1. (a) A “chirp,”or signal whose frequency increases continuously with time. (b) Same
signal after it has passed through the notch filter (13.5.15). The parameter /epsilon1is here 0.2.
Figure 13.5.1 shows the results of using a filter of the form (13.5.15) on a “chirp”input
signal, one that glides upwards in frequency, crossing the notch frequency along the way.
While the bilinear transformation may seem very general, its applications are limited
by some features of the resulting filters. The method is good at getting the general shape
of the desired filter, and good where “flatness”is a desired goal. However, the nonlinear
mapping between wandfmakes it dif ficult to design to a desired shape for a cutoff, and
may move cutoff frequencies (de fined by a certain number of dB) from their desired places.
Consequently, practitionersoftheartofdigital filterdesignreservethebilineartransformation
for specific situations, and arm themselves with a variety of other tricks. We suggest that
you do likewise, as your projects demand.
CITED REFERENCES AND FURTHER READING:
Hamming, R.W. 1983, Digital Filters , 2nd ed. (Englewood Cliffs, NJ: Prentice-Hall).
Antoniou, A. 1979, Digital Filters: Analysis and Design (New York: McGraw-Hill).
Parks, T.W., and Burrus, C.S. 1987, Digital Filter Design (New York: Wiley).
Oppenheim, A.V., andSchafer, R.W. 1989, Discrete-Time Signal Processing (EnglewoodCliffs,
NJ: Prentice-Hall).
Rice, J.R. 1964, The Approximation of Functions (Reading, MA: Addison-Wesley); also 1969,
op. cit., Vol. 2.
Rabiner,L.R.,andGold,B.1975, TheoryandApplicationofDigitalSignalProcessing (Englewood
Cliffs, NJ: Prentice-Hall).
13.6LinearPredictionandLinearPredictiveCoding 557Sample 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).13.6 Linear Prediction and Linear Predictive
Coding
Webeginwithaverygeneralformulationthatwillallowustomakeconnections
to various special cases. Let {y/prime
α}be a set of measured values for some underlying
set of true values of a quantity y, denoted {yα}, related to these true values by
the addition of random noise,
y/prime
α=yα+nα (13.6.1 )
(compareequation 13.3.2, with a somewhat different notation). Our use of a Greek
subscript to index the members of the set is meant to indicate that the data points
are not necessarily equally spaced along a line, or even ordered: they might be“random”pointsinthree-dimensionalspace,forexample. Now,supposewewantto
constructthe “best”estimate ofthe truevalueof someparticularpoint y
⋆as alinear
combination of the known, noisy, values. Writing
y⋆=/summationdisplay
αd⋆αy/prime
α+x⋆ (13.6.2 )
we want to find coefficients d⋆αthat minimize, in some way, the discrepancy x⋆.
Thecoefficients d⋆αhavea“star”subscripttoindicatethattheydependonthechoice
of point y⋆. Later, we might want to let y⋆be one of the existing yα’s. In that case,
our problem becomes one of optimal filtering or estimation, closely related to the
discussion in §13.3. On the other hand, we might want y⋆to be a completely new
point. In that case, our problem will be one of linear prediction .
A natural way to minimize the discrepancy x⋆is in the statistical mean square
sense. Ifanglebracketsdenotestatisticalaverages,thenweseek d⋆α’sthatminimize
/angbracketleftbig
x2
⋆/angbracketrightbig
=/angbracketleftBigg/bracketleftbigg/summationdisplay
αd⋆α(yα+nα)−y⋆/bracketrightbigg2/angbracketrightBigg
=/summationdisplay
αβ(/angbracketleftyαyβ/angbracketright+/angbracketleftnαnβ/angbracketright)d⋆αd⋆β−2/summationdisplay
α/angbracketlefty⋆yα/angbracketrightd⋆α+/angbracketleftbig
y2
⋆/angbracketrightbig(13.6.3 )
Here we have used the fact that noise is uncorrelatedwith signal, e.g., /angbracketleftnαyβ/angbracketright=0.
The quantities /angbracketleftyαyβ/angbracketrightand/angbracketlefty⋆yα/angbracketrightdescribe the autocorrelation structure of the
underlying data. We have already seen an analogous expression, (13.2.2), for thecase of equally spaced data points on a line; we will meet correlation several times
againinitsstatisticalsenseinChapters14and15. Thequantities /angbracketleftn
αnβ/angbracketrightdescribethe
autocorrelationpropertiesof the noise. Often, forpoint-to-pointuncorrelatednoise,
we have /angbracketleftnαnβ/angbracketright=/angbracketleftbig
n2
α/angbracketrightbig
δαβ. It is convenient to think of the various correlation
quantities as comprising matrices and vectors,
φαβ≡/angbracketleftyαyβ/angbracketright φ⋆α≡/angbracketlefty⋆yα/angbracketright ηαβ≡/angbracketleftnαnβ/angbracketrightor/angbracketleftbig
n2
α/angbracketrightbig
δαβ (13.6.4 )
Setting the derivative of equation (13.6.3) with respect to the d⋆α’s equal to zero,
one readily obtains the set of linear equations,
/summationdisplay
β[φαβ+ηαβ]d⋆β=φ⋆α (13.6.5 )