23 void wavefft(
double a[],
double b[],
int ntot,
int n,
int nspan,
int isn)
108 double c72=0.30901699437494742;
109 double s72=0.95105651629515357;
110 double s120=0.86602540378443865;
111 double pi=3.141592653589793;
113 double c1,
c2=0.,
c3=0., s1, s2=0., s3=0.;
115 if (isn < 0) { s72=-s72; s120=-s120; rad=-rad; inc=-inc; }
122 double radf=rad*double(jc)*0.5;
125 for (
int i=1; i<nf; i++) nfac[i]=0;
130 while ( ((
k/16)*16) ==
k ) { m++; nfac[
m]=4;
k/=16; }
135 while ( (
k%jj) == 0) { m++; nfac[
m]=
j;
k/=jj ;}
138 }
while ( jj <=
k ) ;
140 if (
k <= 4 ) { kt=
m; nfac[m+1]=
k;
if (
k != 1) m++; }
143 if ( ((
k/4)*4) ==
k ) { m++; nfac[
m]=2;
k/=4; }
148 if( (
k%j) == 0) { m++; nfac[
m]=
j;
k/=
j; }
154 { j=kt;
do { m++; nfac[
m]=nfac[
j]; j--; }
while ( j != 0 ) ; }
157 for (
int i=1; i<nf; i++)
158 if ( nfac[i] > maxf) { maxf=nfac[
i]; }
161 int kk, k1, k2, k3=0, k4;
163 double aa, bb, ak, bk, aj, bj;
164 double ajm, ajp, akm, akp, bjm, bjp, bkm, bkp;
168 double *at, *ck, *bt, *sk;
169 at=
new double[maxf+1];
170 bt=
new double[maxf+1];
171 ck=
new double[maxf+1];
172 sk=
new double[maxf+1];
178 sd=radf/double(kspan);
179 cd=sin(sd); cd*=
cd; cd*=2.;
185 if ( nfac[i] != 2 )
goto L400;
192 a[k2]=a[kk]-ak; b[k2]=b[kk]-bk;
193 a[kk]=a[kk]+ak; b[kk]=b[kk]+bk;
196 }
while (kk <= nn-1);
200 }
while ( kk <= jc-1);
201 if ( kk > (kspan -1) )
goto L800 ;
210 ak=a[kk]-a[k2]; bk=b[kk]-b[k2];
211 a[kk]=a[kk]+a[k2]; b[kk]=b[kk]+b[k2];
212 a[k2]=c1*ak-s1*bk; b[k2]=s1*ak+c1*bk;
233 kk=(k1-kspan+1)/2+jc-1;
235 }
while (kk <= jc+jc-1);
243 { k1=kk+kspan; k2=k1+kspan;
245 aj=a[k1]+a[k2]; bj=b[k1]+b[k2];
246 a[kk]=ak+aj; b[kk]=bk+bj;
247 ak=-0.5*aj+ak; bk=-0.5*bj+bk;
248 aj=(a[k1]-a[k2])*s120; bj=(b[k1]-b[k2])*s120;
249 a[k1]=ak-bj; b[k1]=bk+aj;
250 a[k2]=ak+bj; b[k2]=bk-aj;
253 }
while ( kk < nn -1 );
257 }
while ( kk <= kspan -1 );
263 if ( nfac[i] != 4 )
goto L600;
270 k1=kk+kspan; k2=k1+kspan; k3=k2+kspan;
271 akp=a[kk]+a[k2]; akm=a[kk]-a[k2];
272 ajp=a[k1]+a[k3]; ajm=a[k1]-a[k3];
275 bkp=b[kk]+b[k2]; bkm=b[kk]-b[k2];
276 bjp=b[k1]+b[k3]; bjm=b[k1]-b[k3];
279 if ( isn < 0 )
goto L450;
280 akp=akm-bjm; akm=akm+bjm;
281 bkp=bkm+ajm; bkm=bkm-ajm;
282 if ( s1 == 0 )
goto L460;
284 a[k1]=akp*c1-bkp*s1; b[k1]=akp*s1+bkp*
c1;
285 a[k2]=ajp*c2-bjp*s2; b[k2]=ajp*s2+bjp*
c2;
286 a[k3]=akm*
c3-bkm*s3; b[k3]=akm*s3+bkm*
c3;
289 if ( kk <= nt-1 )
goto L420;
293 c1=2.0-(c2*c2+s1*s1);
302 if ( kk <= kspan-1)
goto L420;
305 if ( kk <= jc-1 )
goto L410;
306 if ( kspan == jc )
goto L800;
311 akp=akm+bjm; akm=akm-bjm;
312 bkp=bkm-ajm; bkm=bkm+ajm;
314 if ( s1 != 0)
goto L430;
317 a[k1]=akp; b[k1]=bkp;
318 a[k2]=ajp; b[k2]=bjp;
319 a[k3]=akm; b[k3]=bkm;
322 if ( kk <= nt-1 )
goto L420;
332 { k1=kk+kspan; k2=k1+kspan; k3=k2+kspan; k4=k3+kspan;
333 akp=a[k1]+a[k4]; akm=a[k1]-a[k4];
334 bkp=b[k1]+b[k4]; bkm=b[k1]-b[k4];
335 ajp=a[k2]+a[k3]; ajm=a[k2]-a[k3];
336 bjp=b[k2]+b[k3]; bjm=b[k2]-b[k3];
338 a[kk]=aa+akp+ajp; b[kk]=bb+bkp+bjp;
339 ak=akp*c72+ajp*c2+aa; bk=bkp*c72+bjp*c2+bb;
340 aj=akm*s72+ajm*s2; bj=bkm*s72+bjm*s2;
341 a[k1]=ak-bj; a[k4]=ak+bj;
342 b[k1]=bk+aj; b[k4]=bk-aj;
343 ak=akp*c2+ajp*c72+aa; bk=bkp*c2+bjp*c72+bb;
344 aj=akm*s2-ajm*s72; bj=bkm*s2-bjm*s72;
345 a[k2]=ak-bj; a[k3]=ak+bj;
346 b[k2]=bk+aj; b[k3]=bk-aj;
348 }
while ( kk < nn-1 );
350 }
while ( kk <= kspan-1 ) ;
359 if (
k == 3 )
goto L320;
360 if (
k == 5 )
goto L510;
361 if (
k == jf )
goto L640;
364 c1=cos(s1); s1=sin(s1);
365 if ( jf > maxf )
goto L998 ;
366 ck[jf]=1.0; sk[jf]=0.0;
369 { ck[
j]=ck[
k]*c1+sk[
k]*s1;
370 sk[
j]=ck[
k]*s1-sk[
k]*
c1;
406 { k1+=kspan; k2-=kspan;
420 if ( jj > jf ) jj-=jf;
425 a[k1]=ak-bj; b[k1]=bk+aj;
426 a[k2]=ak+bj; b[k2]=bk-aj;
432 }
while ( kk <= nn-1 );
434 }
while ( kk <= kspan-1 );
437 if ( i == m )
goto L800;
451 a[kk]=c2*ak-s2*b[kk];
452 b[kk]=s2*ak+c2*b[kk];
454 }
while ( kk <= nt-1 );
460 }
while ( kk <= kspnn-1 );
464 c1=2.0-(c2*c2+s1*s1);
469 }
while ( kk <= kspan -1 );
473 }
while ( kk <= jc+jc-1 );
481 if ( kt == 0)
goto L890;
487 { np[j+1]=np[
j]/nfac[
j];
488 np[
k]=np[
k+1]*nfac[
j];
498 if ( n != ntot )
goto L850;
510 }
while ( k2 < ks-1 );
516 }
while ( k2 > np[j]-1 );
518 L840:
if ( kk < k2 )
goto L820;
521 if ( k2 < ks-1 )
goto L840;
522 if ( kk < ks-1 )
goto L830;
531 { ak=a[kk]; a[kk]=a[k2]; a[k2]=ak;
532 bk=b[kk]; b[kk]=b[k2]; b[k2]=bk;
534 }
while ( kk <
k-1 );
538 }
while ( kk < nt-1 );
542 }
while ( k2 < ks-1 );
544 do { k2-=np[
j]; j++; k2+=np[j+1]; }
while ( k2 > np[j]-1 );
548 if ( kk < k2 )
goto L850;
551 if ( k2 < ks-1 )
goto L880;
552 if ( kk < ks-1 )
goto L870;
555 if ( 2*kt+1 >= m ) {
delete []
np;
delete [] at;
delete [] bt;
556 delete [] ck ;
delete [] sk;
return;}
562 do { nfac[
j]=nfac[
j]*nfac[j+1]; j--; }
while ( j != kt );
568 delete []
np; np=
new int[maxp+1]; }
580 if ( jj-1 >= k2 )
goto L902;
587 if ( j <= nn )
goto L904;
592 do {
k=kk+1; kk=np[
k]-1; np[
k]=-kk-1; }
while ( kk != j-1 );
596 do { j++; kk=np[
j]-1; }
while ( kk < -1 );
598 if ( kk != j-1 )
goto L910;
600 if ( j != nn )
goto L914;
606 {
do { j--; }
while ( np[j] < 0 );
611 if ( jj > maxf ) kspan=maxf;
617 do { at[k2]=a[k1]; bt[k2]=b[k1]; k2++; k1-=inc; }
while ( k1 != kk );
623 do { a[k1]=a[k2]; b[k1]=b[k2]; k1-=inc; k2-=inc; }
while( k1 != kk );
629 do { a[k1]=at[k2]; b[k1]=bt[k2]; k2++; k1-=inc; }
while ( k1 != kk );
637 if ( nt >= 0 )
goto L924;
638 delete []
np;
delete [] at;
delete [] bt;
delete [] ck ;
delete [] sk;
642 delete []
np;
delete [] at;
delete [] bt;
delete [] ck ;
delete [] sk;
644 printf(
"Error: array bounds exceeded within subroutine wavefft.\n");
printf("total live time: non-zero lags = %10.1f \n", liveTot)
wavearray< double > a(hp.size())
void wavefft(double a[], double b[], int ntot, int n, int nspan, int isn)