67 lines
1.6 KiB
C
67 lines
1.6 KiB
C
|
|
#include <math.h>
|
|
#include "nrutil.h"
|
|
#define MOD(a,b) while(a >= b) a -= b;
|
|
#define MACC 4
|
|
|
|
void fasper(x,y,n,ofac,hifac,wk1,wk2,nwk,nout,jmax,prob)
|
|
float *prob,hifac,ofac,wk1[],wk2[],x[],y[];
|
|
unsigned long *jmax,*nout,n,nwk;
|
|
{
|
|
void avevar(),realft(),spread();
|
|
unsigned long j,k,ndim,nfreq,nfreqt;
|
|
float ave,ck,ckk,cterm,cwt,den,df,effm,expy,fac,fndim,hc2wt;
|
|
float hs2wt,hypo,pmax,sterm,swt,var,xdif,xmax,xmin;
|
|
|
|
*nout=0.5*ofac*hifac*n;
|
|
nfreqt=ofac*hifac*n*MACC;
|
|
nfreq=64;
|
|
while (nfreq < nfreqt) nfreq <<= 1;
|
|
ndim=nfreq << 1;
|
|
if (ndim > nwk) nrerror("workspaces too small in fasper");
|
|
avevar(y,n,&ave,&var);
|
|
if (var == 0.0) nrerror("zero variance in fasper");
|
|
xmin=x[1];
|
|
xmax=xmin;
|
|
for (j=2;j<=n;j++) {
|
|
if (x[j] < xmin) xmin=x[j];
|
|
if (x[j] > xmax) xmax=x[j];
|
|
}
|
|
xdif=xmax-xmin;
|
|
for (j=1;j<=ndim;j++) wk1[j]=wk2[j]=0.0;
|
|
fac=ndim/(xdif*ofac);
|
|
fndim=ndim;
|
|
for (j=1;j<=n;j++) {
|
|
ck=(x[j]-xmin)*fac;
|
|
MOD(ck,fndim)
|
|
ckk=2.0*(ck++);
|
|
MOD(ckk,fndim)
|
|
++ckk;
|
|
spread(y[j]-ave,wk1,ndim,ck,MACC);
|
|
spread(1.0,wk2,ndim,ckk,MACC);
|
|
}
|
|
realft(wk1,ndim,1);
|
|
realft(wk2,ndim,1);
|
|
df=1.0/(xdif*ofac);
|
|
pmax = -1.0;
|
|
for (k=3,j=1;j<=(*nout);j++,k+=2) {
|
|
hypo=sqrt(wk2[k]*wk2[k]+wk2[k+1]*wk2[k+1]);
|
|
hc2wt=0.5*wk2[k]/hypo;
|
|
hs2wt=0.5*wk2[k+1]/hypo;
|
|
cwt=sqrt(0.5+hc2wt);
|
|
swt=SIGN(sqrt(0.5-hc2wt),hs2wt);
|
|
den=0.5*n+hc2wt*wk2[k]+hs2wt*wk2[k+1];
|
|
cterm=SQR(cwt*wk1[k]+swt*wk1[k+1])/den;
|
|
sterm=SQR(cwt*wk1[k+1]-swt*wk1[k])/(n-den);
|
|
wk1[j]=j*df;
|
|
wk2[j]=(cterm+sterm)/(2.0*var);
|
|
if (wk2[j] > pmax) pmax=wk2[(*jmax=j)];
|
|
}
|
|
expy=exp(-pmax);
|
|
effm=2.0*(*nout)/ofac;
|
|
*prob=effm*expy;
|
|
if (*prob > 0.01) *prob=1.0-pow(1.0-expy,effm);
|
|
}
|
|
#undef MOD
|
|
#undef MACC
|