71 lines
1.7 KiB
C
71 lines
1.7 KiB
C
|
|
#include <math.h>
|
|
#define NRANSI
|
|
#include "nrutil.h"
|
|
#define M 64
|
|
#define NDFT 1024
|
|
#define MPOL 6
|
|
#define TWOPI (2.0*3.14159265)
|
|
|
|
void dftint(float (*func)(float), float a, float b, float w, float *cosint,
|
|
float *sinint)
|
|
{
|
|
void dftcor(float w, float delta, float a, float b, float endpts[],
|
|
float *corre, float *corim, float *corfac);
|
|
void polint(float xa[], float ya[], int n, float x, float *y, float *dy);
|
|
void realft(float data[], unsigned long n, int isign);
|
|
static int init=0;
|
|
int j,nn;
|
|
static float aold = -1.e30,bold = -1.e30,delta,(*funcold)(float);
|
|
static float data[NDFT+1],endpts[9];
|
|
float c,cdft,cerr,corfac,corim,corre,en,s;
|
|
float sdft,serr,*cpol,*spol,*xpol;
|
|
|
|
cpol=vector(1,MPOL);
|
|
spol=vector(1,MPOL);
|
|
xpol=vector(1,MPOL);
|
|
if (init != 1 || a != aold || b != bold || func != funcold) {
|
|
init=1;
|
|
aold=a;
|
|
bold=b;
|
|
funcold=func;
|
|
delta=(b-a)/M;
|
|
for (j=1;j<=M+1;j++)
|
|
data[j]=(*func)(a+(j-1)*delta);
|
|
for (j=M+2;j<=NDFT;j++)
|
|
data[j]=0.0;
|
|
for (j=1;j<=4;j++) {
|
|
endpts[j]=data[j];
|
|
endpts[j+4]=data[M-3+j];
|
|
}
|
|
realft(data,NDFT,1);
|
|
data[2]=0.0;
|
|
}
|
|
en=w*delta*NDFT/TWOPI+1.0;
|
|
nn=IMIN(IMAX((int)(en-0.5*MPOL+1.0),1),NDFT/2-MPOL+1);
|
|
for (j=1;j<=MPOL;j++,nn++) {
|
|
cpol[j]=data[2*nn-1];
|
|
spol[j]=data[2*nn];
|
|
xpol[j]=nn;
|
|
}
|
|
polint(xpol,cpol,MPOL,en,&cdft,&cerr);
|
|
polint(xpol,spol,MPOL,en,&sdft,&serr);
|
|
dftcor(w,delta,a,b,endpts,&corre,&corim,&corfac);
|
|
cdft *= corfac;
|
|
sdft *= corfac;
|
|
cdft += corre;
|
|
sdft += corim;
|
|
c=delta*cos(w*a);
|
|
s=delta*sin(w*a);
|
|
*cosint=c*cdft-s*sdft;
|
|
*sinint=s*cdft+c*sdft;
|
|
free_vector(cpol,1,MPOL);
|
|
free_vector(spol,1,MPOL);
|
|
free_vector(xpol,1,MPOL);
|
|
}
|
|
#undef M
|
|
#undef NDFT
|
|
#undef MPOL
|
|
#undef TWOPI
|
|
#undef NRANSI
|