13.4 Power Spectrum Estimation Using the FFT 557 significantly better than the reduction of about K/2 that would have resulted if the same number of data points were segmented without overlapping. We can now codify these ideas into a routine for spectral estimation.While we generally avoid input/output coding,we make an exception here to show how data are read sequentially in one pass through a data file (referenced through the parameter 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 the component P(fe)at the Nyquist frequency.It would also be straightforward to include that component. #include <math.h> #include <stdio.h> #include "nrutil.h" granted for 19881992 #define WINDOW(j,a,b)(1.0-fabs((((j)-1)-(a))*(b)))/*Bartlett * /*#define WINDOW(j,a,b)1.0 *//Square * 100 /*#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) to any from NUMERICAL RECIPES IN 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,.. based on (2*k+1)*m data points(if ovrlap is set true(1))or 4*k*m data points (if ovrlap (Nort is set false(0)).The number of segments of the data is 2*k in both cases: The routine calls four1 k times,each call with 2 partitions each of 2*m real data points. America server computer, to make one paper Cambridge University Press. THE void fouri(float data[],unsigned long nn,int isign); 是 ART int mm,m44,m43,m4,kk,joffn,joff,j2,j; float w,facp,facm,*w1,*w2,sumw=0.0,den=0.0; Programs mm=m+m; Useful factors m43=(m4=mm+mm)+3: st st m44=m43+1; wi-vector(1,m4); Copyright (C) w2-vector(1,m); facmem: facp=1.0/m; for (j=1;j<=mm;j++)sumw +SQR(WINDOW(j,facm,facp)); to directcustser OF SCIENTIFIC COMPUTING(ISBN 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]); v@cambri for (kk=1;kk<=k;kk++){ Loop over data set segments in groups of two. 1988-1992 by Numerical Recipes 10-621 for (joff =-1;joff<=0;joff++){ 43106 Get two complete segments into workspace. if (ovrlap){ for (i=1;i<=m;++)w1[ioff+i+i]=w2[i]; for (j=1;j<=m;j++)fscanf(fp,"%f",&w2[j]); (outside joffn=joff+mm; Software. for (j=1;j<=m;j++)v1[joffn+j+j]=v2[j]; North else for (j=joff+2;j<=m4;j+=2) Amer fscanf (fp,"%f",&w1[j]); for (j=1;j<=mm;j++){ Apply the window to the data j2=j+i: w=WINDOW(j,facm,facp); 1[j2]*=w; w1[j2-1]*= four1(w1,mm,1); Fourier transform the windowed data. p[1]+=(SQR(w1[1])+SQR(w1[2]); Sum results into previous segments.13.4 Power Spectrum Estimation Using the FFT 557 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). significantly better than the reduction of about K/2 that would have resulted if the same number of data points were segmented without overlapping. We can now codify these ideas into a routine for spectral estimation. While we generally avoid input/output coding, we make an exception here to show how data are read sequentially in one pass through a data file (referenced through the parameter 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 the component P(fc) at the Nyquist frequency. It would also be straightforward to include that component. #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 ovrlap is set false (0)). The number of segments of the data is 2*k in both cases: The routine calls four1 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