66 lines
1.4 KiB
C
66 lines
1.4 KiB
C
|
|
#include <math.h>
|
|
#include "nrutil.h"
|
|
#define M 64
|
|
#define NDFT 1024
|
|
#define MPOL 6
|
|
#define TWOPI (2.0*3.14159265)
|
|
|
|
void dftint(func,a,b,w,cosint,sinint)
|
|
float (*func)(),*cosint,*sinint,a,b,w;
|
|
{
|
|
void dftcor(),polint(),realft();
|
|
static int init=0;
|
|
int j,nn;
|
|
static float aold = -1.e30,bold = -1.e30,delta,(*funcold)();
|
|
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
|