13.4 Power Spectrum Estimation Using the FFT 553 leakage and data windowing which were brought up after equation(13.4.7)above. Data Windowing The purpose of data windowing is to modify equation(13.4.7),which expresses the relation between the spectral estimate P at a discrete frequency and the actual underlying continuous spectrum P(f)at nearby frequencies.In general,the spectral power in one"bin"k contains leakage from frequency components that are actually s bins away,where s is the independent variable in equation(13.4.7).There is,as we pointed out,quite substantial leakage even from moderately large values of s. g When we select a run of N sampled points for periodogram spectral estimation, nted for we are in effect multiplying an infinite run of sampled data cj by a window function in time,one that is zero except during the total sampling time NA,and is unity during that time.In other words,the data are windowed by a square window function.By the convolution theorem(12.0.9;but interchanging the roles of f and t),the Fourier transform of the product of the data with this square window function is equal to the convolution of the data's Fourier transform with the window's Fourier transform.In fact,we determined equation(13.4.7)as nothing more than the square of the discrete Fourier transform of the unity window function. University Press. 需 W(s)= 1 sin(ms) (13.4.8) 9 sin(πs/N) 心量 OF SCIENTIFIC( The reason for the leakage at large values of s,is that the square window function turns on and off so rapidly.Its Fourier transform has substantial components at high frequencies.To remedy this situation,we can multiply the input data ci,j=0,...,N-1 by a window function wj that changes more gradually from zero to a maximum and then back to zero as j ranges from 0 to N.In this case,the equations for the periodogram estimator(13.4.4-13.4.5)become 10.621 N-1 Dk三∑C四e2ikN 公 k=0,..,N-1 (13.4.9) Numerical Recipes 43106 j=0 P(0)=P(fo)-IDo (outside Software. PU)=DP+Dw-] North k=1,2, )PUa (13.4.10) 显 where Wss stands for "window squared and summed," N-1 (13.4.11) i=013.4 Power Spectrum Estimation Using the FFT 553 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). leakage and data windowing which were brought up after equation (13.4.7) above. Data Windowing The purpose of data windowing is to modify equation (13.4.7), which expresses the relation between the spectral estimate Pk at a discrete frequency and the actual underlying continuous spectrum P(f) at nearby frequencies. In general, the spectral power in one “bin” k contains leakage from frequency components that are actually s bins away, where s is the independent variable in equation (13.4.7). There is, as we pointed out, quite substantial leakage even from moderately large values of s. When we select a run of N sampled points for periodogram spectral estimation, we are in effect multiplying an infinite run of sampled data c j by a window function in time, one that is zero except during the total sampling time N∆, and is unity during that time. In other words, the data are windowed by a square window function. By the convolution theorem (12.0.9; but interchanging the roles of f and t), the Fourier transform of the product of the data with this square window function is equal to the convolution of the data’s Fourier transform with the window’s Fourier transform. In fact, we determined equation (13.4.7) as nothing more than the square of the discrete Fourier transform of the unity window function. W(s) = 1 N2  sin(πs) sin(πs/N) 2 = 1 N2      N
−1 k=0 e2πisk/N      2 (13.4.8) The reason for the leakage at large values of s, is that the square window function turns on and off so rapidly. Its Fourier transform has substantial components at high frequencies. To remedy this situation, we can multiply the input data cj , j = 0,...,N − 1 by a window function wj that changes more gradually from zero to a maximum and then back to zero as j ranges from 0 to N. In this case, the equations for the periodogram estimator (13.4.4–13.4.5) become Dk ≡ N
−1 j=0 cjwj e2πijk/N k = 0,...,N − 1 (13.4.9) P(0) = P(f0) = 1 Wss |D0| 2 P(fk) = 1 Wss |Dk| 2 + |DN−k| 2  k = 1, 2,..., N 2 − 1  P(fc) = P(fN/2) = 1 Wss  DN/2   2 (13.4.10) where Wss stands for “window squared and summed,” Wss ≡ N N
−1 j=0 w2 j (13.4.11)