/*

fftwbw.c 

A Fast Fourier Routine for computing power, cospectra and quadrature spectra
for meteorological variables.

This program adapts the subroutine twofft.c  from Press et al.s Numerical
Recipes in C to compute coefficients to compute power spectra
and real and imaginary components of the cross spectrum.  The twofft subroutine
was modified and renamed fftxy. This was needed as there was need to extract sub components of
FFT coefficients to compute co and quadrature spectra.


This program also adapts features from Fortran FFT program of Carter and
Ferrier to average spectral densities across spectral bands

This program uses many Numerical C utilities such as vector, matrix, free_vector and
free-matrix so it needs to be linked with the  nrutil.c.  We also need the nr.h and
nrutil.h include files. 


There are some nuances with using NR FFT.  First, if one is to use the four1 subroutine one will need
to input a sequence of real and imaginary components.  data(re1, im1, re2, im2...ren, im2) yielding an
array that is 2n.

The program twofft packs the two real arrays into one complex one, then runs four1.  The data are set up as

data(s1_1, s2_1, s1_2, s2_2, s1_3, s2_3...s1_n, s2_n).

The program outputs to FFT arrays.  FFT1 has the real and imaginary components of the first variable

Sre(s1)_1, Sim(s1)_1, Sre(s1)_2, Sim(s1)_2, ...Sre(s1)_n, Sim(s1)_n. FFT2 has the re and im components
of the second variable.


A few tricks with FFt

S(k)=y(k)y*(k), where y*(k) is the complex conjugate.  With a fft y(N-k) equal y*(k)


1-5-2000

Dennis Baldocchi
ESPM
Ecosystem Science Division
UC Berkeley
345 Hilgard Hall
Berkeley, CA 94720

baldocchi@nature.berkeley.edu

510-642-2874


DEBUG NOTES

1-5-2000
Found errors in the cospectra code.
Also need to normalize the power spectra by the sum of the window squared.  Otherwise
the integral of the spectra and the variance were not equal, according to tests
with sine waves

12-7-99

Had terrible debug problems even though the code was set up ok at home on NT.  It came
down to needing to up the STACK option on the LINKER.  It is now set at 10000000

In the meantime I ended up re-installing WATCOM 4 times as I lost the editor and was chasing my
tail as declared values would instantly change numbers as memory was overwritten. What a pain.

Finally with the correct STACK I was able to also add the SMOOTH routine to AVESPEC.  This is a
.1, .2, .4, .2, .1 filter on the spectral estimates.  It was gleaned from the old Carter and Ferrier
spectra.

12-4-99

Debugged version. Am able to read the file.raw files and process data.  It was hanging up on some arrays

*/



#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <io.h>
#include <dos.h>
#include <float.h>
#include <string.h>
#include <malloc.h>

/* include files for Numerical Recipes also remember to compile and link NRUTIL.C */

#include "nr.h"
#include "nrutil.h"

#define nch 6  /* this version is for OSUFL, other studies may have 7 channels */
#define nisamp 16384

#define NRANSI  

/* #define WINDOW(j,a,b) (1.0-fabs((((j)-1)-(a))*(b)))        Bartlett Window */

 #define WINDOW(j) 0.5*(1.0-cos((float)j*2*3.14159/(float)nsamp)) /* Hann cosine window */


/* declare subroutines */

void lremv (float xx[], float nnn, float iswch, float dc, float slope);  
void INDATA(float xa[], float ya[], int i1, int i2, float *wx);
void fftxy(float x1[], float y1[], float a[], float b[],unsigned long ni,
           float m[], float n[]);                                                                                                      

void avespec(float p[], float px[], float f[], float x, float y, float z, int *i);


/* declare global variables */

                                                                                                                      

const float nsamp=16384; /* number of samples closest to a factor of 2 */
const int isamp=16384;
const float   pi=3.14159;

                                                                                                                                
/*                                                                                                                                  
File pointers for i/o
*/                                                                                                                                  
                                                                                                                                    
FILE *fptr1, *fptr2, *fptr10;


    main ()
        {

/* pointer declaration for arrays that are passed */
float *powerx, *powery, *cospec_re, *cospec_im, *data1, *data2;
float *fftx, *ffty, *nfreq;
float *newpowerx, *newpowery, *newcovar_re, *newcovar_im;

float dx, sx, dy, sy, meanx, meany, tempa, tempb, tempc, tempd;
float real_x, imagin_x, realconj_x, imconj_x, real_y, imagin_y, realconj_y, imconj_y;
float samprate, delt,facp,facm,w, varx, vary,covar, df, flow, fhigh, ubar;
float integSx, integSy, integSxy;
float sumwindow;

long int isamp2;  /* 2 times nsamp */
long int halfsamp,isampp1,m44,m43,m4,j2,j,j1,jj;
        
int is3, ich1, ich2, inum, iloop;   
int handle1, handle2, handle10;   

char FNAME[31];

char outname[35], BERNSUB[9], varname[2];

samprate=10.;      /* sampling rate */
delt=1./samprate;  /* time step of samples */
isamp2=isamp*2;    /* twice number of samples */
halfsamp=isamp/2;  /* half number of samples */

df=1./(delt*nsamp);  /* frequency step */

flow=0.0;            /* lowest frequency */
fhigh=samprate/2.;   /* highest frequency, delimited by Shannon limit, one-half sample rate */


/*

Read the data 5 times for wu, wv, wt, wq and wc covariances
or reset iloop if one only wants to read it once and a certain value

*/

for(iloop=5;iloop <=5; iloop++)
{


switch(iloop)
{
   case 1: sprintf(varname,"%s","u");break;   
      case 2: sprintf(varname,"%s","v");break; 
        case 3: sprintf(varname,"%s","t");break; 
           case 4: sprintf(varname,"%s","q");break; 
                case 5: sprintf(varname,"%s","c");break; 
                      
default:
sprintf(varname,"%s","w");
}

       
        sprintf(BERNSUB,"%s","wbw_97");

/*
     OUTPUT DATA FROM TOWER OR FLOOR SYSTEMS

     FILE, outname
*/
      sprintf(outname,"%s%s%s%s%s","d:\\",BERNSUB,"\\wbwsp97",varname,".csv");



/* open the file of names  note they should have names as d:\sub\xdayhour.yr or *.raw.  The file names
do not need \\ to denote the difference between directories as is needed in direct C lines. */


 if(( fptr10=fopen("d:\\wbw_97\\filewbw97.raw","r")) != NULL);

      handle10=fileno(fptr10);         

/* output file */

if((fptr2=fopen(outname,"w"))!= NULL)

 handle2=fileno(fptr2);      

                                    
/* while file.raw is open */

     while(! feof(fptr10))
{

/* read name of input *.raw file */

fscanf(fptr10,"%s",FNAME);

/* open the *.raw file */

if ((fptr1=fopen(FNAME,"rb")) !=NULL)
printf(" %s \n",FNAME);
else
/* need to stop program */
exit(0);


handle1=fileno(fptr1);

fprintf (fptr2, "%s\n",FNAME);
printf("vector p \n");


/* declare the matrices for the power and cross spectra coefficients.
each array has a real and imaginary component, in the 1 to n and n+1 to 2n ranges */
     
        powerx=vector(1,isamp2);
        powery=vector(1,isamp2);
        cospec_re=vector(1,isamp2);
        cospec_im=vector(1,isamp2);

        fftx=vector(1,isamp2);
        ffty=vector(1,isamp2);

        data1=vector(1,nsamp);
        data2=vector(1,nsamp);

        newpowerx=vector(1,halfsamp);
        newpowery=vector(1,halfsamp);
        newcovar_re=vector(1,halfsamp);
        newcovar_im=vector(1,halfsamp);
        nfreq=vector(1,halfsamp);
        

/* initialize the vectors */

        meanx=0.0;
        meany=0.0;

        for (j=1;j<=isamp2;j++)
        {
        powerx[j]=0.0;
        powery[j]=0.0;
        cospec_re[j]=0.0;
        cospec_im[j]=0.0;
        }

        for (j=1; j <=halfsamp; j++)
        {
            newpowerx[j]=0.0;
            newpowery[j]=0.0;
            newcovar_re[j]=0.0;
            newcovar_im[j]=0.0;
            nfreq[j]=0.0;
        }

/*
read the data from a binary file, 2 byte integers. The number of channels depends on the
experiment. also input the index of the channels on which the FFT will be performed

note:
ch0: w
ch1: u
ch2: v
ch3: T
ch4: q
ch5: CO2
ch6: ???

*/


ich1=0;       /* The first variable of the cospectrum, in this case w */
ich2=iloop;  /* incremented by iloop, for u, v, t, q, c */

/* Velocity coordinates are rotated in INDATA */


     INDATA(data1,data2,ich1,ich2, &ubar);  


/* compute means */


for (j=1;j<=nsamp;j++)
{


/*

testing with sin waves

data1[j]=sin(2*pi*j/180);
data2[j]=sin(2*pi*j/360);
*/

meanx +=data1[j];
meany +=data2[j];
}


meanx /=(float)nsamp;
meany /=(float)nsamp;

/* subroutine to remove trends */


 /*
        remove the linear trend and compute the variance
        if is3 = 0, do not remove dc component or slope
        = 1, remove the dc component
        > 1, remove the dc component and slope
        */
        
        is3 = 1;   
        dx = 0.0;
        sx = 0.0;
        dy = 0.0;
        sy = 0.0;

      lremv (data1, nsamp, is3, dx, sx);
      lremv (data2, nsamp, is3, dx, sx);


/*
compute variances and co-variances on the deviations from the mean. This
is performed after the lremv subroutine
*/

varx=0.0;
vary=0.0;
covar=0.0;
sumwindow=0.0;

for (j=1;j<=nsamp;j++)
{
varx += data1[j]*data1[j];
vary += data2[j]*data2[j];
covar += data1[j]*data2[j];
}

varx /=(nsamp-1.);
vary /=(nsamp-1.);
covar /=(nsamp-1.);

fprintf (fptr2, "meanx, meany,varx,vary,covar\n");
fprintf (fptr2, "%9.5f,%9.5f,%9.5f,%9.5f, %9.5f\n",meanx, meany,varx,vary,covar);


facm=halfsamp;
facp=1.0/facm;

       
/* apply a cosine window to data matrices */


        isampp1=isamp+1;

        for (j=1;j<=halfsamp;j++)
        {
         w=WINDOW(j);
         data1[j] *=w;     
         data1[isampp1-1] *=w;
        }

        for (j=1;j<=halfsamp;j++)
        {
         w=WINDOW(j);
         data2[j] *=w;
         data2[isampp1-1] *=w;
        }

         for (j=1;j<=isamp;j++)
        {w=WINDOW(j);
        sumwindow +=w;
        }

printf("call fftxy \n");

                fftxy(data1,data2,fftx,ffty,nsamp,cospec_re,cospec_im);

                m43=(m4=isamp2+isamp2)+3;
                m44=m43+1;
                
            
                /*

                compute spectral density for the x variable as the product of
                the fourier transform and its complex conjugate
                fftx(k) and fft*(k).  Remember the fft from N/2 to N is the complex
                conjugate so, fft(N-k) equals fft*(k).  We have to be careful
                as 1,3,5,7... are real components and 2,4,6,.. are imaginary

                We are interested in S(0) to S(N/2)

                */
                  
                powerx[1]=SQR(fftx[1]+fftx[2]);
                powery[1]=SQR(ffty[1]+ ffty[2]);
              
                
                for (j=2;j<=halfsamp;j++)
                {
                jj=j+j;
                
                /* the real component, 0.5(a(k)+a(N-k),associated with indices 1,3,5,7... */

                real_x= fftx[jj-1];  

                /* the imaginary component, i (b(k)-b(N-k))), associated with indices 2,4,6,..*/

                imagin_x=fftx[jj];   

                realconj_x=fftx[isamp2-jj-1]; /* the real complex conjugate, N-1, N-3, N-5.. */

                imconj_x=fftx[isamp2-jj];    /* the imaginary component of the complex conjugate, N, N-2, N-4.. */
/*
The power spectral density, Sx(k)=dt/N Fx(k) F^*x(k),
where F^*x(k) is the complex conjugate, such that a+ib has a-ib as its complex conjugate
its magnitude is Re^2 + Im^2. Sx is technically divided by 4 but this is factored into
the fftx and ffty coefficients of the twofft program.
*/


                powerx[j] += SQR(real_x)+SQR(imagin_x);

                /* the real component, 0.5(a(k)+a(N-k),associated with indices 1,3,5,7... */
                real_y= ffty[jj-1];

                /* the imaginary component, i (b(k)-b(N-k))), associated with indices 2,4,6,..*/
                 
                imagin_y=ffty[jj];   

                realconj_y=ffty[isamp2-jj-1]; /* the real complex conjugate, N-1, N-3, N-5.. */

                imconj_y=ffty[isamp2-jj];    /* the imaginary component of the complex conjugate, N, N-2, N-4.. */


/*
The power spectral density, Sy(k)=dt/N Fy(k) F*y(k),
where F*y(k) is the complex conjugate, such that a+ib has a-ib as its complex conjugate
its magnitude is Re^2 + Im^2. Sy is technically divided by 4 but this is factored into
the fftx and ffty coefficients of the twofft program.
*/
                
                powery[j] += SQR(real_y)+SQR(imagin_y); 



/* wrong cospec  need to work out algebra */

/*
The cospectral density, Syx(k)=dt/N Fx(k) F*y(k),
where F*y(k) is the complex conjugate, such that a+ib has a-ib as its complex conjugate
its magnitude is Re portion dt/2N.  The real part is dt/2N(a(k)b(N-k) +a(N-k)b(k))
*/

               
                                                      
                }

        for (j=1;j<=halfsamp;j++)
        {
        /* powerx[j] *= delt/nsamp;   note the factor 0.25 is already incorporated in the fft coeficients */
        /* powery[j] *= delt/nsamp;   note the factor 0.25 is already incorporated in the fft coeficients */

        powerx[j] *= delt/nsamp;
        powery[j] *= delt/nsamp;
                
        cospec_re[j] *= delt*0.5/nsamp;
        cospec_im[j] *= delt*0.25/nsamp; 
        }
  


/* integrate spectra */

integSx=0.;
integSy=0.;
integSxy=0.;

   for (j=1;j<=halfsamp;j++)
   {
integSx += powerx[j];
integSy += powery[j];
integSxy += cospec_re[j];
   }

   integSx *= (df)*2.; 
   integSy *= (df)*2.;
   integSxy *= (df)*2.;

fprintf (fptr2, "integSx,integSy,integSxy \n");

fprintf (fptr2, "%9.5f, %9.5f, %9.5f \n",integSx,integSy,integSxy);



/*
be careful to process only half nsamp as the FFT only produces 1/2 nsamp Fourier coefficients
and this program pads the second half of the matrix with the imaginary component
*/


printf("call avespec \n");
        
 avespec(powerx, newpowerx, nfreq, df, flow, fhigh, &inum);
 avespec(powery, newpowery, nfreq, df, flow, fhigh, &inum);
 avespec(cospec_re, newcovar_re, nfreq, df, flow, fhigh, &inum);
 avespec(cospec_im, newcovar_im, nfreq, df, flow, fhigh, &inum);  



fprintf (fptr2, "Freq/U, nSxx, nSyy, nSxy_Re, nSxy_Im, Freq \n");

   for (j=1;j<=inum;j++)
        {

/* normalize spectra by frequency and plot as a function of wave number f/u */

         tempa=nfreq[j]*newpowerx[j]/varx;
         tempb=nfreq[j]*newpowery[j]/vary;
         tempc=nfreq[j]*newcovar_re[j]/covar;
         tempd=nfreq[j]*newcovar_im[j]/covar;
            
       fprintf (fptr2, "%9.5f, %9.5f, %9.5f, %9.5f, %9.5f, %9.5f\n",nfreq[j]/ubar,tempa, tempb, tempc, tempd,nfreq[j]);
        } 



        free_vector(data1,1,isamp2);
        free_vector(data2,1,isamp2);
        free_vector(powerx,1,isamp2);
        free_vector(powery,1,isamp2);
        free_vector(cospec_re,1,isamp2);
        free_vector(fftx,1,isamp2);
        free_vector(ffty,1,isamp2);
        free_vector(newpowerx,1,halfsamp);
        free_vector(newpowery,1,halfsamp);
        free_vector(newcovar_re,1,halfsamp);
        free_vector(newcovar_im,1,halfsamp);
        free_vector(nfreq,1,halfsamp);





printf("end of loop \n");

} /* end of while for reading file.raw */
        
fclose(fptr10);
fclose(fptr2);

} /* next iloop  */


/* end of main */
}

         

void lremv (float *xx, float nnn, float iswch, float dc, float slope)
{

int i;

  float fln = nnn;
  dc = 0;
  slope = 0;
  
  for (i = 1; i <= nnn; i ++)
  {
      /* compute sums for mean */
      
    (dc) += xx[i];
    (slope) += (xx[i] * (float)i);
  }


  (dc) /= fln;
  slope = 12.0 * (slope) / (fln* (fln*fln - 1.0)) - 6.0 * (dc) / (fln - 1.0);

  if (iswch > 1)
  {
    fln = dc - 0.5 * (fln + 1.0) * (slope);
    for (i = 1; i <= nnn; i ++) 
    xx[i] -= ((float)i * (slope) - fln);
   }

  if (iswch == 1)
  {
      /* compute deviations from the mean */
      for (i = 1; i <= nnn; i ++)
      xx[i] -= (dc);
  }
  
  return;
}


void INDATA(float *xd, float *yd, int i1, int i2, float *wndspd)
{
short int IXX[nch];       // make two-byte integer
long int freebyte, II, isamppl1;
int J, icount, i, j, handle;

float x_bar[nch];
float xdata[nisamp][nch];

float ws, tw, ce, se, ct, st, theta, eta, ce2,se2,ct2,st2, w,u,v;

printf("%d\n",nch);

isamppl1=isamp+1;



      II=1;

        handle=fileno(fptr1);
        freebyte=filelength(handle);

        for (j=1; j <= isamp; j++)
        for (i=1; i < nch; i++)
        xdata[j][i]=0.0;


while (! feof(fptr1) && II <=isamp) 
{


    /* read binary file, 2 byte words */
    
 fread(IXX,2,nch,fptr1);

/*
     Convert wind velocity from integer form and cm/s to float and m/s 
*/
        for(J=0; J< 3;J++)
        {
        xdata[II][J]=(long double)IXX[J]/100.;
        }

/*
      Convert mv integer to volt float
*/
       xdata[II][3]=(long double)IXX[3]/100.;
       xdata[II][4]=(long double)IXX[4]/1000.;
       xdata[II][5]=(long double)IXX[5]/1000.;
/*        xdata[II][6]=(long double)IXX[6]/1000.;  */
        
/*
    W, U AND V ARE CH. 0, 1 AND 2, RESPECTIVELY.
    T, Q, CO2 ARE CH 3, 4, 5
*/


II++;

}



fclose(fptr1); 


        icount=II;

               
       /* compute means and use this information to rotate the velocity coordinates */
        
        for (j=0;j<nch;j++)
{
        x_bar[j]=0;
               
        for (i=0; i < icount; i++)
        x_bar[j] += xdata[i][j];            /* compute mean       */
       

        x_bar[j] /= (float)icount;
}                      


/* compute angles for rotations based on the horizontal and total wind speed */

                ws = sqrt(x_bar[1]*x_bar[1] + x_bar[2]*x_bar[2]);
                tw = sqrt(x_bar[0]*x_bar[0] + ws*ws);

                *wndspd=ws;


                ce = x_bar[1]/ws;               /* cos eta      */
                se = x_bar[2]/ws;               /* sin eta      */
                ct = ws/tw;                     /* cos theta    */
                st = x_bar[0]/tw;               /* sin theta    */
                theta = 180.*acos(ct)/pi;                /* vertical       */
                eta = 180.*acos(ce)/pi;                  /* horizontal     */
                ce2 = ce*ce;
                se2 = se*se;
                ct2 = ct*ct;
                st2 = st*st;

        for (i=0; i < icount; i++)
        {
                w= xdata[i][0]*ct - xdata[i][1]*st*ce - xdata[i][2]*st*se;  /* rotated w */
                u= xdata[i][1]*ct*ce + xdata[i][2]*ct*se + xdata[i][0]*st;  /* rotated u */ 
                v= xdata[i][2]*ce - xdata[i][1]*se;                      /* rotated v */

             
             xdata[i][0]= w;
             xdata[i][1]=u;
             xdata[i][2]=v;

/* fill vectors with data to be rotated  */

        xd[i]=xdata[i][i1];
        yd[i]=xdata[i][i2];

         }



/*
end of first input read loop
*/

return;
}


void fftxy(float data1[], float data2[], float fft1[], float fft2[],
         unsigned long n, float *cosp, float *quadsp)
      {
        void four1(float data[], unsigned long nn, int isign);
        unsigned long nn3,nn2,nn,jj,j;
        float rep,rem,aip,aim;
        long int isamp_2, ii;

        float *an, *bn, *anconj, *bnconj; 
        
        nn3=1+(nn2=2+n+n);

        nn=n+n;
        
        an=vector(1,nn3);                /* code I added */
        bn=vector(1,nn3);
        anconj=vector(1,nn3);
        bnconj=vector(1,nn3);         
     

/*
this bit of code packs the two data sets into a 2N variable data set with
data sets 1 and 2 alternating
*/

        for (j=1,jj=2;j<=n;j++,jj+=2) {
                fft1[jj-1]=data1[j];
                fft1[jj]=data2[j];
        }


/* call the FFT program */

        four1(fft1,n,1);

        fft2[1]=fft1[2];
        fft1[2]=fft2[2]=0.0;

        an[1]=fft1[j];          
        anconj[1]=fft1[nn-1];
        bn[1]=fft1[2];
        bnconj[j]=fft1[nn];

ii=0;

        for (j=3;j<=n+1;j+=2)
     {

        ii +=1;

 /*
 the computations here correspond with the notes of Rick Eckman
 describe an FFT
 */


 /*
 identify the a(k), a(N-k), b(k) and b(N-k) coefficients before
 they get reassigned,

 a(1:1),b(1:2),a(2:3), b(2:4),..a(N/2,N-1), b(N/2:N)
 */
 
        an[ii]=fft1[j];          
        anconj[ii]=fft1[nn2-j];

        bn[ii]=fft1[j+1];
        bnconj[ii]=fft1[nn3-j];

        /* these assignments are based on my notes from Eckman */
        
                rep=0.5*(fft1[j]+fft1[nn2-j]);     /* re component of the transform of x, Fx */
                rem=0.5*(fft1[j]-fft1[nn2-j]);     /* im component of the transform of y, Fy */

                aip=0.5*(fft1[j+1]+fft1[nn3-j]);   /* re component of Fy */
                aim=0.5*(fft1[j+1]-fft1[nn3-j]);   /* im component of Fx */
                

        /* spectral coefficients for Fx */
        
                fft1[j]=rep;         /* re component of Fx */
                fft1[j+1]=aim;       /* im component of Fx */

        /* these represent the complex conjugates */

                fft1[nn2-j]=rep;      /* re component of the complex conjugate of Fx */
                fft1[nn3-j] = -aim;   /* im component of the complex conjugate of Fx */


        /* spectral coefficients for Fy */
        
                fft2[j]=aip;         /* re component of Fy */
                fft2[j+1] = -rem;    /* im component of the transform of y, Fy; note there is a neg sign */

        /* these represent the complex conjugates */

                fft2[nn2-j]=aip;     /* re component of the complex conjugate of Fy */
                fft2[nn3-j]=rem;     /* im component of the complex conjugate of Fy */

        }


           isamp_2=isamp/2;



/*
cospectrum and quadrature spectrum. The coefficients correspond with values in Fortran
code too
*/
  for (j=1;j<=isamp_2+1;j++)
  {
    cosp[j]=an[j]*bnconj[j]+anconj[j]*bn[j]; 
    quadsp[j]=an[j]*an[j]-anconj[j]*anconj[j]+bn[j]*bn[j]-bnconj[j]*bnconj[j];

  }


                free_vector(an,1,nn3);
                free_vector(anconj,1,nn3);
                free_vector(bn,1,nn3);
                free_vector(bnconj,1,n);
      }

/* (C) Copr. 1986-92 Numerical Recipes Software D..W"<$1##S. */

#include <math.h>
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

void four1(float data[], unsigned long nn, int isign)
{
        unsigned long n,mmax,m,j,istep,i;
        double wtemp,wr,wpr,wpi,wi,theta;
        float tempr,tempi;

        n=nn << 1;
        j=1;
        for (i=1;i<n;i+=2) {
                if (j > i) {
                        SWAP(data[j],data[i]);
                        SWAP(data[j+1],data[i+1]);
                }
                m=n >> 1;
                while (m >= 2 && j > m) {
                        j -= m;
                        m >>= 1;
                }
                j += m;
        }
        mmax=2;
        while (n > mmax) {
                istep=mmax << 1;
                theta=isign*(6.28318530717959/mmax);
                wtemp=sin(0.5*theta);
                wpr = -2.0*wtemp*wtemp;
                wpi=sin(theta);
                wr=1.0;
                wi=0.0;
                for (m=1;m<mmax;m+=2) {
                        for (i=m;i<=n;i+=istep) {
                                j=i+mmax;
                                tempr=wr*data[j]-wi*data[j+1];
                                tempi=wr*data[j+1]+wi*data[j];
                                data[j]=data[i]-tempr;
                                data[j+1]=data[i+1]-tempi;
                                data[i] += tempr;
                                data[i+1] += tempi;
                        }
                        wr=(wtemp=wr)*wpr-wi*wpi+wr;
                        wi=wi*wpr+wtemp*wpi+wi;
                }
                mmax=istep;
        }
}
#undef SWAP
/* (C) Copr. 1986-92 Numerical Recipes Software D..W"<$1##S. */



void avespec(float *powerin, float *powerout, float *frq, float df, float flow, float fhigh, int *inumb)
{

int istrt, istop,j, k, n,  inum, jp1, jp2, jm1, jm2;
long int isampm3, isampm2, isampm1;
float fbeg, sum1, sum2, freq, xpow;
float *xsmooth;  

xsmooth=vector(1,isamp); 

istrt=(int)(flow/df)+1;
istop=(int)(fhigh/df);
fbeg=(float)((int)(flow/df)*df);


isampm1=isamp-1;
isampm2=isamp-2;
isampm3=isamp-3;

xsmooth[1]=0.4*powerin[1]+0.4*powerin[2] + 0.2*powerin[3];
xsmooth[2]=0.4*powerin[2]+0.2*powerin[4] + 0.2*(powerin[1]+powerin[3]);
xsmooth[isampm1]=0.4*powerin[isampm1]+0.2*(powerin[isamp]+powerin[isampm2]) + 0.2*powerin[isampm3];
xsmooth[isamp]=0.4*powerin[isamp]+0.4*powerin[isampm1] + 0.2*powerin[isampm2];

for (j=3;j<=isampm3;j++)
{   jp1=j+1;
    jp2=j+2;
    jm1=j-1;
    jm2=j-2;
    xsmooth[j]=0.4*powerin[j]+.2*(powerin[jp1]+powerin[jm1])+.1*(powerin[jp2]+powerin[jm2]);
}

inum=0;
sum1=0;
sum2=0;
j=0;


for(k=istrt; k<=istop;k++)
{
    
j +=1;
n=k/10;
if (n < 1) n=1;

freq=fbeg+ df*(float)(k-istrt);

sum1 +=freq;
sum2 +=xsmooth[k];  /* powerin[] without smoothing */


if (j==n)
{
inum += 1;
freq=sum1/n;
xpow=sum2/n;

j=0;
sum1=0;
sum2=0;

powerout[inum]=xpow;
frq[inum]=freq;
}
 }

*inumb=inum;

free_vector(xsmooth,1,isamp);

return;
}