c13-4 (Numerical Recipes in C), страница 3

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).(13.4.12)55513.4 Power Spectrum Estimation Using the FFT1square windowWelch window.8amplitude.4.2Hann window0050100150bin number200250Figure 13.4.1. Window functions commonly used in FFT power spectral estimation.

The data segment,here of length 256, is multiplied (bin by bin) by the window function before the FFT is computed. Thesquare window, which is equivalent to no windowing, is least recommended. The Welch and Bartlettwindows are good choices.much as a factor of two, however), but trade this off by widening the leakage tailby a significant factor (e.g., the reciprocal of 10 percent, a factor of ten).

If wedistinguish between the width of a window (number of samples for which it is atits maximum value) and its rise/fall time (number of samples during which it risesand falls); and if we distinguish between the FWHM (full width to half maximumvalue) of the leakage function’s main lobe and the leakage width (full width thatcontains half of the spectral power that is not contained in the main lobe); thenthese quantities are related roughly byN(window width)N(leakage width in bins) ≈(window rise/fall time)(FWHM in bins) ≈(13.4.16)(13.4.17)For the windows given above in (13.4.13)–(13.4.15), the effective windowwidths and the effective window rise/fall times are both of order 12 N .

Generallyspeaking, we feel that the advantages of windows whose rise and fall times areonly small fractions of the data length are minor or nonexistent, and we avoid usingthem. One sometimes hears it said that flat-topped windows “throw away less ofthe data,” but we will now show you a better way of dealing with that problem byuse of overlapping data segments.Let us now suppose that we have chosen a window function, and that we areready to segment the data into K segments of N = 2M points.

Each segment willbe FFT’d, and the resulting K periodograms will be averaged together to obtain aSample 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).Bartlett window.6556Chapter 13.Fourier and Spectral Applications1Hann.6BartlettWelch.4.2square0−8−6−4−202offset in units of frequency bins468Figure 13.4.2.Leakage functions for the window functions of Figure 13.4.1.

A signal whosefrequency is actually located at zero offset “leaks” into neighboring bins with the amplitude shown. Thepurpose of windowing is to reduce the leakage at large offsets, where square (no) windowing has largesidelobes. Offset can have a fractional value, since the actual signal frequency can be located betweentwo frequency bins of the FFT.PSD estimate at M + 1 frequency values from 0 to fc .

We must now distinguishbetween two possible situations. We might want to obtain the smallest variancefrom a fixed amount of computation, without regard to the number of data pointsused. This will generally be the goal when the data are being gathered in real time,with the data-reduction being computer-limited. Alternatively, we might want toobtain the smallest variance from a fixed number of available sampled data points.This will generally be the goal in cases where the data are already recorded andwe are analyzing it after the fact.In the first situation (smallest spectral variance per computer operation), it isbest to segment the data without any overlapping. The first 2M data points constitutesegment number 1; the next 2M data points constitute segment number 2; and soon, up to segment number K, for a total of 2KM sampled points.

The variance inthis case, relative to a single segment, is reduced by a factor K.In the second situation (smallest spectral variance per data point), it turns outto be optimal, or very nearly optimal, to overlap the segments by one half of theirlength. The first and second sets of M points are segment number 1; the secondand third sets of M points are segment number 2; and so on, up to segment numberK, which is made of the Kth and K + 1st sets of M points.

The total number ofsampled points is therefore (K +1)M , just over half as many as with nonoverlappingsegments. The reduction in the variance is not a full factor of K, since the segmentsare not statistically independent. It can be shown that the variance is instead reducedby a factor of about 9K/11 (see the paper by Welch in [3]). This is, however,Sample 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).amplitude of leakage.813.4 Power Spectrum Estimation Using the FFT557#include <math.h>#include <stdio.h>#include "nrutil.h"#define WINDOW(j,a,b) (1.0-fabs((((j)-1)-(a))*(b))) /* Bartlett *//* #define WINDOW(j,a,b) 1.0 */ /* Square *//* #define WINDOW(j,a,b) (1.0-SQR((((j)-1)-(a))*(b))) */ /* Welch */void spctrm(FILE *fp, float p[], int m, int k, int ovrlap)Reads data from input stream specified by file pointer fp and returns as p[j] the data’s power(mean square amplitude) at frequency (j-1)/(2*m) cycles per gridpoint, for j=1,2,...,m,based on (2*k+1)*m data points (if ovrlap is set true (1)) or 4*k*m data points (if ovrlapis set false (0)).

The number of segments of the data is 2*k in both cases: The routine callsfour1 k times, each call with 2 partitions each of 2*m real data points.{void four1(float data[], unsigned long nn, int isign);int mm,m44,m43,m4,kk,joffn,joff,j2,j;float w,facp,facm,*w1,*w2,sumw=0.0,den=0.0;mm=m+m;Useful factors.m43=(m4=mm+mm)+3;m44=m43+1;w1=vector(1,m4);w2=vector(1,m);facm=m;facp=1.0/m;for (j=1;j<=mm;j++) sumw += SQR(WINDOW(j,facm,facp));Accumulate the squared sum of the weights.for (j=1;j<=m;j++) p[j]=0.0;Initialize the spectrum to zero.if (ovrlap)Initialize the “save” half-buffer.for (j=1;j<=m;j++) fscanf(fp,"%f",&w2[j]);for (kk=1;kk<=k;kk++) {Loop over data set segments in groups of two.for (joff = -1;joff<=0;joff++) {Get two complete segments into workspace.if (ovrlap) {for (j=1;j<=m;j++) w1[joff+j+j]=w2[j];for (j=1;j<=m;j++) fscanf(fp,"%f",&w2[j]);joffn=joff+mm;for (j=1;j<=m;j++) w1[joffn+j+j]=w2[j];} else {for (j=joff+2;j<=m4;j+=2)fscanf(fp,"%f",&w1[j]);}}for (j=1;j<=mm;j++) {Apply the window to the data.j2=j+j;w=WINDOW(j,facm,facp);w1[j2] *= w;w1[j2-1] *= w;}four1(w1,mm,1);Fourier transform the windowed data.p[1] += (SQR(w1[1])+SQR(w1[2]));Sum results into previous segments.Sample 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).significantly better than the reduction of about K/2 that would have resulted if thesame number of data points were segmented without overlapping.We can now codify these ideas into a routine for spectral estimation. Whilewe generally avoid input/output coding, we make an exception here to show howdata are read sequentially in one pass through a data file (referenced through theparameter FILE *fp).

Only a small fraction of the data is in memory at any one time.Note that spctrm returns the power at M , not M + 1, frequencies, omitting thecomponent P (fc ) at the Nyquist frequency. It would also be straightforward toinclude that component.558Chapter 13.Fourier and Spectral Applicationsfor (j=2;j<=m;j++) {j2=j+j;p[j] += (SQR(w1[j2])+SQR(w1[j2-1])+SQR(w1[m44-j2])+SQR(w1[m43-j2]));}den += sumw;Correct normalization.Normalize the output.}CITED REFERENCES AND FURTHER READING:Oppenheim, A.V., and Schafer, R.W. 1989, Discrete-Time Signal Processing (Englewood Cliffs,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, Fourier Transforms and Their Physical Applications (New York: Academic Press).Elliott, D.F., and Rao, K.R. 1982, Fast Transforms: Algorithms, Analyses, Applications (NewYork: Academic Press).Bloomfield, P. 1976, Fourier Analysis of Time Series – An Introduction (New York: Wiley).Rabiner, L.R., and Gold, B. 1975, Theory and Application of Digital Signal Processing (EnglewoodCliffs, NJ: Prentice-Hall).13.5 Digital Filtering in the Time DomainSuppose that you have a signal that you want to filter digitally. For example, perhapsyou want to apply high-pass or low-pass filtering, to eliminate noise at low or high frequenciesrespectively; 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 60Hz power-line interference, you may need a notch filter to remove only a narrow band aroundthat frequency. This section speaks particularly about the case in which you have chosen todo such filtering in the time domain.Before continuing, we hope you will reconsider this choice. Remember how convenientit 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 filtereddata 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 andnegative frequencies, and that the magnitude of the frequency extremes is alwaysthe Nyquist frequency 1/(2∆), where ∆ is the sampling interval.