The Use of Fast Fourier Transform for the Estimation
`of Power Spectra: A Method Based on Time Aver-(cid:173)
`aging Over Short, Modified Periodograms
`
`PETER D. WELCH
`
`Abstract-The use of the fast Fourier transform in power spec(cid:173)
`trwn analysis is described. Principal advantages of this method are a
`reduction in the number of computations and in required core
`storage, and convenient application in nonstationarity tests. The
`method involves sectioning the record and averaging modified
`periodograms of the sections.
`
`INTRODUCTION
`
`THIS PAPER outlines a method for the application
`
`of the fast Fourier transform algorithm to the
`estimation of power spectra, which involves sec(cid:173)
`tionmg the record, taking modified periodograms of
`these sections, and averaging these modified periodo(cid:173)
`grams. In many instances this method. involves fewer
`computations than other methods. Moreover, it in(cid:173)
`the transformation of sequences which are
`volves
`shorter than the whole record which is an advantage
`when computations are to be performed on a machine
`with limited core storage. Finally, it directly yields a
`potential resolution in the time dimension which is use(cid:173)
`ful for testing and measuring nonstationarity. As will be
`pointed out, it is closely related to the method of com(cid:173)
`plex demodulation described by Bingham, Godfrey, and
`Tukey. 1
`
`THE METHOD
`Let X(j), j=O, · · ·, N-1 be a sample from a sta(cid:173)
`tionary, second-order stochastic sequence. Assume for
`simplicity that E(X) =0. Let X(j) have spectral density
`P(j), Iii :::;;t. We take segments, possibly overlapping, of
`length L with the starting points of these segments D
`units apart. Let X1(j), j=O, · · · , L-1 be the first such
`segment. Then
`
`Similarly,
`X 2(j) = X(j + D)
`and finally
`
`j = 0, · · ·, L - 1.
`
`j = 0, · · ·, L - 1,
`
`XK(j)-= X(j+ (K - 1)D)
`
`j = 0, · · ·, L - 1.
`
`Manuscript received February 28, 1967.
`The author is with the IBM Research Center, Yorktown Heights,
`N. Y.
`1 C. Bingham, M. D. Godfrey, and J. W. Tukey, "Modern tech(cid:173)
`niques of power spectrum estimation," this issue, p. 56--66.
`
`We suppose we have K such segments; X 1(j), · · · ,
`XK(j), and that they cover the entire record, i.e., that
`(K-1)D+L =N. This segmenting is illustrated in Fig. 1.
`The method of estimation is as follows. For each
`segment of length L we calculate a modified periodo(cid:173)
`gram. That is, we select a data window W(j),j=O, · · ·,
`L-1, and
`sequences X1(j) W(j), · · · ,
`form
`the
`X K (j) W (j). We then take the finite Fourier transforms
`A1(n), · · · , AK(n) of these sequences. Here
`1 L-1
`Ak(n) = - L Xk(j)W(j)e-2kiin!L
`L
`i=O
`and i = ( -1) 1' 2• Finally, we obtain the K modified
`periodograms
`
`k = 1, 2, · · ·, K,
`
`where
`
`and
`
`n
`fn = (cid:173)
`L
`
`n = 0, · · ·, L/2
`
`1 L-1
`U = - L W 2(j).
`L j=O
`
`The spectral estimate is the average of these periodo(cid:173)
`grams, i.e.,
`
`1 K
`F(Jn) = - L h(f n) •
`K k=l
`
`Now one can show that
`
`E{ F(fn)} =
`
`f 1/2
`
`h(j)P(f -
`-1/2
`
`fn)dj
`
`where
`
`and
`
`L W(j)e2rifj
`1 I L-1
`h(f) = -
`LU j=O
`
`1
`
`2
`
`f 1/2
`
`h(f)df = 1.
`-1/2
`
`Reprinted from IEEE Trans. Audio and Electroacoust., vol. AU-15, pp. 70-73, June 1967.
`
`17
`
`1
`
`Micro Motion 1049
`
`

`
`WELCH: FAST FOURIER TRANSFORM FOR POWER SPECTRA
`
`J _______ X_(_j) ______ J
`0
`N-l
`I x,(j)
`0
`
`L-l
`I
`X2(j)
`D+L+l
`
`D
`
`~-------·-------·----------'
`
`Fig. 1.
`
`Illustration of record segmentation.
`
`Hence, we have a spectral estimator F(j) with a resul(cid:173)
`tant spectral window whose area is unity and whose
`width is of the order of 1 / L.
`
`CHOICE OF DATA WINDOWS
`We suggest two reasonable choices for the data win(cid:173)
`dow WU); one of them has the shape 1-t2 : -1~t~1
`and gives a spectral window which, when the two are
`normalized to have the same half-power width, is very
`close in shape to the hanning or cosine arch spectral
`window; the other data window has the shape 1- It I :
`-1 ~ t ~ 1 and gives the Parzen spectral window. The
`actual functions for a particular segment length L are
`
`For ~(j) the half-power width is
`
`(1.28)
`!l2f:::<--•
`L+l
`
`THE VARIANCES OF THE ESTIMATES
`As developed above our estimator is given by
`
`A
`
`1 K
`P(fn) = - L h(fn)'
`K k=l
`
`(n = 0, 1, · · · , L/2).
`
`Now, if we let
`
`d(j) = Covariance {!,.(Jn), l1;+J(fn)}
`
`then it is easily shown that
`
`1 {
`K-1 K - j }
`d(O) + 2 L --d(j) .
`Var {P(fn)} = -
`K
`K
`
`;-1
`
`Further, if
`
`p(j) = Correlation { Ik(fn), h+J(fn)}
`
`d(j)
`d(O)
`
`then,
`
`j
`}
`K-1 K -
`d(O) {
`1+2 L --p(j)
`Var {P(fn)} = -
`K
`K
`=Var{!,.,(!,.)} {l + 2 E K - j p(j)}.
`
`K
`
`;-1
`
`K
`
`;-1
`
`and
`
`L-1
`j - - -
`2
`L+l
`2
`L-1
`j - - -
`2
`L+l
`2
`
`j = 0, 1, · · · , L - 1.
`
`Assume now that X(j) is a sample from a Gaussian
`process and assume that P(f) is flat over the passband
`of our estimator. Then we can show2 that
`
`The resultant spectral windows corresponding to these
`data windows are given approximately by
`[sin { (L + l)rf}
`1 {
`2
`h1(f) "'" LU r 2(L + l)f2
`(L + l)rf
`- cos{ (L + l)rf} ]}
`{ (L + l)rf/2}} 2
`h __ 1 {(L + 1) sin2
`{ (L + l)rf/2} 2
`2(f) - LU
`2
`In the preceding approximations Lis a scale parameter.
`In changing L we change the shape of h1(f) and ~(j) only
`in stretching or shrinking the horizontal dimension. For
`h1(j) the half-power width is
`
`2
`
`(1.16)
`!:J.if"'"--.
`L+l
`
`Further, under the above assumptions and assuming
`that h(f-fn) =0 for f <0 and/>! we can show3 that
`
`]2
`]2 [ L-1
`[ L-1
`.E W(k)W(k + jD)
`.E W 2(k)
`
`p(j) =
`
`.
`
`Hence, we have the following result which enables us to
`estimate the variances of F(f,.) when f,. is not close to 0
`or!.
`Result: If X(j) is a sample from a Gaussian process,
`and P(D is flat over the passband of the estimator, and
`h(f-f,.) =0 for f <0 andf>!, then
`
`j
`}
`p2(j, ) {
`K-1 K -
`Var {P(fn)} = _n_ 1+2 L --p(j)
`K
`K
`
`;=1
`
`1 P. D. Welch, "A direct digital method of power spectrum estima(cid:173)
`tion," IBM J. Res. and Dev., vol. 5, pp. 141-156, April 1961.
`• In Welch2 we obtained the variance spectrum of I ,.(f,.) considered
`as a function of time. The above result is obtained by taking the
`Fourier transform of this spectrum.
`
`18
`
`2
`
`

`
`IEEE TRANSACTIONS ON AUDIO AND ELECTROACOUSTICS, JUNE 1967
`
`where
`
`]2
`]2/ [ L-1
`[ L-1
`p(j) = Eo W(k) W(k + jD) E W2(k)
`
`For estimating the spectrum of P(j,.) at 0 and ! the
`variance is twice as great, as given by the following
`result:
`Result: If X(j) is a sample from a Gaussian process
`and P(f) is flat over the passband of the estimator, then
`Var I P(o or 1/2)}
`j
`}
`1 K. -
`2P2(0 or 1/2) {
`1 + 2 :E - - p(j)
`K
`K
`.
`
`K-
`
`-1
`
`=
`
`where p(j) is as defined above.
`In the above results note that p(j) ~O and that
`p(j) = 0 if D ~L. Hence, if we average over K segments
`the best we can do is obtain a reduction of the variance
`by a factor 1/K. Further, this 1/K reduction can be
`achieved (under these conditions) if we have nonover(cid:173)
`lapping segments. Hence, if the total number of points
`N can be made sufficiently large the computationally
`most efficient procedure for achieving any desired vari(cid:173)
`ance is to have nonoverlapping segments, i.e., to let
`D =L. In this case we have
`
`}
`
`P 2(j,.)L
`
`·
`
`P 2(f,.)
`I A
`Var P(f,.) = - - =
`K
`N
`Further, under these conditions E { F(f,.)} = E { h(f,.)}
`=P(j,.) and, hence,
`
`£21 P(f,.)}
`- - - - -=K
`Var I P(j,.)}
`and the equivalent degrees of freedom of the approxi(cid:173)
`mating chi-square distribution is given by
`E.D.F. { P(f,.)} = 2K.
`If the total number of points N cannot be made arbi(cid:173)
`trarily large, and we wish to get a near maximum reduc(cid:173)
`tion in the variance out of a fixed number of points then
`a reasonable procedure is to overlap the segments by
`one half their length, i.e., to let D=L/2. In this case, if
`we use W1(j) as the data window we get p(l) = 1/9 and
`p(j) =0 for j> 1. Letting F1(f,.) be the estimate, we have
`
`=----
`9K
`
`The factor 11/9, compared with the factor 1.0 for non(cid:173)
`overlapped segments, inflates the variance. However, an
`overall reduction in variance for fixed record length is
`achieved because of the difference in the value of K. For
`nonoverlapped segments we have K = N / L; for the
`overlapping discussed here
`
`2N
`2N
`N
`K= - - 1 = - -1 =-·
`L/2
`L
`L
`
`Therefore, for fixed N and L the overall reduction in
`variance achieved by this overlapping is by a factor of
`11/18. Now again E { F1(j,.)} =P(j,.) and, hence,
`
`E2{ F1(/n)}
`Var { F1(/,.)}
`
`18N
`9K
`=-:::;,--·
`11
`11L
`
`(1.16)
`'11(/) ""' - - ""' -
`L + 1
`6L
`
`7
`
`.
`
`Finally,
`
`Thus,
`
`and the equivalent degrees of freedom of the approxi(cid:173)
`mating chi-square distribution is
`E.D.F. { F1(j,.)} = 2.8N '1f.
`Similarly, if we use W2(j) as our data window we get
`p(1) = 1/16 and p(j) =0, j> 1. Letting F2(j,.) be the
`estimate in this case we get, by following the above
`steps and using the result tl2f= (1.28)/(L+t), that the
`equivalent degrees of freedom is again approximately
`
`Thus, both W1(j) and W2(j) yield roughly the same
`variance when adjusted to have windows of equal half
`power width. Finally, we should point out that the
`above variances need to be doubled and the equivalent
`degrees of freedom halved for the points! n = 0 and !.
`
`DETAILS IN THE APPLICATION OF THE FAST
`FOURIER TRANSFORM ALGORITHM
`Our estimator F(j,.) is given by
`1 K
`L
`
`F(f,.) = K Ei h(f,.) = UK k~ I Ak(n) 12,
`
`K
`
`where L is the length of the segments, and K is the
`number of segments into which the record is broken, and
`
`1 L-1
`U = - L W 2(j).
`L j=O
`
`We will first discuss how the complex algorithm can be
`used to obtain the summation Lk-1KI Ak(n) I 2 two terms
`at a time with K/2[or (K+t)/2, if K is odd] rather
`than K transforms. Suppose K is even and let
`
`Y1(j) = X 1(j)W(j) + iX2(j)W(j) l
`
`·
`YK12(j) = XK-1(j)W(j) + iXK(j)W(j)
`
`j = 0, · · · , L - 1.
`
`19
`
`3
`
`

`
`WELCH: FAST FOURIER TRANSFORM FOR POWER SPECTRA
`
`Let Bk(n) be the transform of Yk(j). Then, by the linear(cid:173)
`ity property of the finite Fourier transform
`
`Further,
`Bk(N - n) = A21:-1(1Y - n) + iA2k(N -· n)
`__.___..,
`= A21·-1(n) + iA2k(n).
`____,,.,
`~ow,
`I Bk(n) 12 = (Az1·-1(n) + iA2k(n))(A2k-1(n) -
`i.421.(n))
`____,,
`iA2k(n))(A2k-1(n) + iA2k(n)).
`2 = (A2k-1(n) -
`I B1:(N - n) 1
`These equations yield, with some algebra,
`2 + I Bk(N - n) 12
`= 2( I A2k-1(n) 12 + I Ak(n) 12
`
`I Bk(n) 1
`
`).
`
`Hence, finally,
`L
`F(fn) = - - L (I Bk(n) 12 + I Bk(N - n) 12).
`2UK k~t
`
`K/2
`
`If K is odd this procedure can be extended in an obvious
`fashion by defining Y<K+1i12U) =XK(j) and summing
`from 1 to (K+1)/2.
`A second observation on the actual application of the
`algorithm concerns the bit-inverting. If the algorithm is
`applied as described here, and one is especially con(cid:173)
`cerned with computation time, then the bit-inverting
`could be postponed until after the summation. Thus,
`instead of bit-inverting K/2 times, one would only have
`to bit-invert once.
`
`spectral estimation. Complex demodulation is discussed
`in Tukey, 4 Godfrey,5 and Bingham, Godfrey, and
`Tukey. 1 The functions Ak(n)e-2rikDJL considered as
`functions of k are complex demodulates sampled at the
`sampling period D. In this case the demodulating func(cid:173)
`tion is e-2rifni. A phase coherency from sample to sample
`is retained in the complex demodulates. This phase is
`lost in estimating the spectrum and, hence, as a method
`of estimating spectra, complex demodulation is identical
`to the method of this section. However, additional in(cid:173)
`formation can be obtained from the time variation of the
`phase of the demodulates.
`
`THE SPACING OF THE SPECTRAL ESTIMATES
`This method yields estimates spaced 1/L units apart.
`If more finely spaced estimates are desired zeros can be
`added to the sequences Xk(j) W(j) before taking the
`transforms. If L' zeros are added giving time sequences
`L+L'=M long and we let Ak'(n) be the finite Fourier
`transforms of these extended sequences, i.e.,
`
`1 L-1
`Ak'(n) = - L Xk(j)W(j)e-2rijn/M
`M ;~o
`
`then the modified periodogram is given by
`
`y2
`h(Jn) = -
`LU
`
`I A/(n) 12
`
`where
`
`n
`fn = M
`
`n = 0, 1, · · · , M/2.
`
`COMPUTATION TIME
`The time required to perform a finite Fourier trans(cid:173)
`form on a sequence of length L
`is approximately
`k' L log2 L where k' is a constant which depends upon the
`program and type of computer. Hence, if we overlap
`segments. by an amount L/2 we require an amount of
`computing time (performing two transforms simul(cid:173)
`taneously) approximately equal to
`
`( ~ ) ( L~2 ) k' L log2 L = k' N log2 L,
`
`plus the amount of time required to premultiply by the
`data window and average. If we only consider the time
`required for the Fourier transformation this compares
`with approximately k' N(log2 N)/2 for the smoothing of
`the periodogram. Hence, if L < (N) 1i 2 it requires less
`computing time than the smoothing of the periodogram.
`
`Everything proceeds exactly as earlier except that we
`have estimates spaced at intervals of 1/ M rather than
`1/L.
`
`ESTIMATION OF CROSS SPECTRA
`Let X(j), j=O, · · · ,N-1, and Y(j), j=O,
`N -1, be samples from two second-order stochastic
`sequences. This method can be extended in a straight(cid:173)
`forward manner to the estimation of the cross spectrum,
`P,,u(f). In exactly the same fashion each sample is
`divided in K segments of length L. Call these segments
`X 1(j), · · · , XK(j) and Yi(j), · · · , Y K(j). Modified
`cross periodograms are calculated for each pair of seg(cid:173)
`ments Xk(j), Yk(j), and the average of these modified
`cross periodograms constitutes the estimate Fx11(fn).
`The spectral window is the same as is obtained using
`this method for the estimation of the spectrum.
`
`RELATION OF THIS METHOD TO
`COMPLEX DEMODULATION
`It is appropriate to mention here the process of com(cid:173)
`plex demodulation and its relation to this method of
`
`'J. W. Tukey, "Discussion, emphasizing the connection between
`analysis of variance and spectrum analysis," Technometrics, vol. 3,
`pp. 191-219, May 1961.
`'M. D. Godfrey, "An exploratory study of the bispectrum of
`economic time series," Applied Statistics, vol. 14, pp. 48-69, January
`1965.
`
`20
`
`4