widely known free FFT library.
Dependents: Peach_AudioChannelDividerAndCompensator
fftpack.c
- Committer:
- dokunewon
- Date:
- 2015-10-11
- Revision:
- 0:d7215946c769
File content as of revision 0:d7215946c769:
/* * This file is part of libfftpack. * * libfftpack is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * libfftpack is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with libfftpack; if not, write to the Free Software * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA */ /* * libfftpack is being developed at the Max-Planck-Institut fuer Astrophysik * and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt * (DLR). */ /* fftpack.c : A set of FFT routines in C. Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber (Version 4, 1985). C port by Martin Reinecke (2010) */ #include <math.h> #include <stdlib.h> #include <string.h> #include "mytype.h" #include "fftpack.h" #define WA(x,i) wa[(i)+(x)*ido] #define CH(a,b,c) ch[(a)+ido*((b)+l1*(c))] #define CC(a,b,c) cc[(a)+ido*((b)+cdim*(c))] #define PM(a,b,c,d) { a=c+d; b=c-d; } #define PMC(a,b,c,d) { a.r=c.r+d.r; a.i=c.i+d.i; b.r=c.r-d.r; b.i=c.i-d.i; } #define ADDC(a,b,c) { a.r=b.r+c.r; a.i=b.i+c.i; } #define SCALEC(a,b) { a.r*=b; a.i*=b; } #define CONJFLIPC(a) { FLOAT tmp_=a.r; a.r=-a.i; a.i=tmp_; } /* (a+ib) = conj(c+id) * (e+if) */ #define MULPM(a,b,c,d,e,f) { a=c*e+d*f; b=c*f-d*e; } typedef struct { FLOAT r,i; } cmplx; #define CONCAT(a,b) a ## b #define X(arg) CONCAT(passb,arg) #define BACKWARD #include "fftpack_inc.inc" #undef BACKWARD #undef X #define X(arg) CONCAT(passf,arg) #include "fftpack_inc.inc" #undef X #undef CC #undef CH #define CC(a,b,c) cc[(a)+ido*((b)+l1*(c))] #define CH(a,b,c) ch[(a)+ido*((b)+cdim*(c))] static void radf2 (size_t ido, size_t l1, const FLOAT *cc, FLOAT *ch, const FLOAT *wa) { const size_t cdim=2; size_t i, k, ic; FLOAT ti2, tr2; for (k=0; k<l1; k++) PM (CH(0,0,k),CH(ido-1,1,k),CC(0,k,0),CC(0,k,1)) if ((ido&1)==0) for (k=0; k<l1; k++) { CH( 0,1,k) = -CC(ido-1,k,1); CH(ido-1,0,k) = CC(ido-1,k,0); } if (ido<=2) return; for (k=0; k<l1; k++) for (i=2; i<ido; i+=2) { ic=ido-i; MULPM (tr2,ti2,WA(0,i-2),WA(0,i-1),CC(i-1,k,1),CC(i,k,1)) PM (CH(i-1,0,k),CH(ic-1,1,k),CC(i-1,k,0),tr2) PM (CH(i ,0,k),CH(ic ,1,k),ti2,CC(i ,k,0)) } } static void radf3(size_t ido, size_t l1, const FLOAT *cc, FLOAT *ch, const FLOAT *wa) { const size_t cdim=3; static const FLOAT taur=-0.5, taui=0.86602540378443864676; size_t i, k, ic; FLOAT ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3; for (k=0; k<l1; k++) { cr2=CC(0,k,1)+CC(0,k,2); CH(0,0,k) = CC(0,k,0)+cr2; CH(0,2,k) = taui*(CC(0,k,2)-CC(0,k,1)); CH(ido-1,1,k) = CC(0,k,0)+taur*cr2; } if (ido==1) return; for (k=0; k<l1; k++) for (i=2; i<ido; i+=2) { ic=ido-i; MULPM (dr2,di2,WA(0,i-2),WA(0,i-1),CC(i-1,k,1),CC(i,k,1)) MULPM (dr3,di3,WA(1,i-2),WA(1,i-1),CC(i-1,k,2),CC(i,k,2)) cr2=dr2+dr3; ci2=di2+di3; CH(i-1,0,k) = CC(i-1,k,0)+cr2; CH(i ,0,k) = CC(i ,k,0)+ci2; tr2 = CC(i-1,k,0)+taur*cr2; ti2 = CC(i ,k,0)+taur*ci2; tr3 = taui*(di2-di3); ti3 = taui*(dr3-dr2); PM(CH(i-1,2,k),CH(ic-1,1,k),tr2,tr3) PM(CH(i ,2,k),CH(ic ,1,k),ti3,ti2) } } static void radf4(size_t ido, size_t l1, const FLOAT *cc, FLOAT *ch, const FLOAT *wa) { const size_t cdim=4; static const FLOAT hsqt2=0.70710678118654752440; size_t i, k, ic; FLOAT ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4; for (k=0; k<l1; k++) { PM (tr1,CH(0,2,k),CC(0,k,3),CC(0,k,1)) PM (tr2,CH(ido-1,1,k),CC(0,k,0),CC(0,k,2)) PM (CH(0,0,k),CH(ido-1,3,k),tr2,tr1) } if ((ido&1)==0) for (k=0; k<l1; k++) { ti1=-hsqt2*(CC(ido-1,k,1)+CC(ido-1,k,3)); tr1= hsqt2*(CC(ido-1,k,1)-CC(ido-1,k,3)); PM (CH(ido-1,0,k),CH(ido-1,2,k),CC(ido-1,k,0),tr1) PM (CH( 0,3,k),CH( 0,1,k),ti1,CC(ido-1,k,2)) } if (ido<=2) return; for (k=0; k<l1; k++) for (i=2; i<ido; i+=2) { ic=ido-i; MULPM(cr2,ci2,WA(0,i-2),WA(0,i-1),CC(i-1,k,1),CC(i,k,1)) MULPM(cr3,ci3,WA(1,i-2),WA(1,i-1),CC(i-1,k,2),CC(i,k,2)) MULPM(cr4,ci4,WA(2,i-2),WA(2,i-1),CC(i-1,k,3),CC(i,k,3)) PM(tr1,tr4,cr4,cr2) PM(ti1,ti4,ci2,ci4) PM(tr2,tr3,CC(i-1,k,0),cr3) PM(ti2,ti3,CC(i ,k,0),ci3) PM(CH(i-1,0,k),CH(ic-1,3,k),tr2,tr1) PM(CH(i ,0,k),CH(ic ,3,k),ti1,ti2) PM(CH(i-1,2,k),CH(ic-1,1,k),tr3,ti4) PM(CH(i ,2,k),CH(ic ,1,k),tr4,ti3) } } static void radf5(size_t ido, size_t l1, const FLOAT *cc, FLOAT *ch, const FLOAT *wa) { const size_t cdim=5; static const FLOAT tr11= 0.3090169943749474241, ti11=0.95105651629515357212, tr12=-0.8090169943749474241, ti12=0.58778525229247312917; size_t i, k, ic; FLOAT ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3, dr4, dr5, cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5; for (k=0; k<l1; k++) { PM (cr2,ci5,CC(0,k,4),CC(0,k,1)) PM (cr3,ci4,CC(0,k,3),CC(0,k,2)) CH(0,0,k)=CC(0,k,0)+cr2+cr3; CH(ido-1,1,k)=CC(0,k,0)+tr11*cr2+tr12*cr3; CH(0,2,k)=ti11*ci5+ti12*ci4; CH(ido-1,3,k)=CC(0,k,0)+tr12*cr2+tr11*cr3; CH(0,4,k)=ti12*ci5-ti11*ci4; } if (ido==1) return; for (k=0; k<l1;++k) for (i=2; i<ido; i+=2) { ic=ido-i; MULPM (dr2,di2,WA(0,i-2),WA(0,i-1),CC(i-1,k,1),CC(i,k,1)) MULPM (dr3,di3,WA(1,i-2),WA(1,i-1),CC(i-1,k,2),CC(i,k,2)) MULPM (dr4,di4,WA(2,i-2),WA(2,i-1),CC(i-1,k,3),CC(i,k,3)) MULPM (dr5,di5,WA(3,i-2),WA(3,i-1),CC(i-1,k,4),CC(i,k,4)) PM(cr2,ci5,dr5,dr2) PM(ci2,cr5,di2,di5) PM(cr3,ci4,dr4,dr3) PM(ci3,cr4,di3,di4) CH(i-1,0,k)=CC(i-1,k,0)+cr2+cr3; CH(i ,0,k)=CC(i ,k,0)+ci2+ci3; tr2=CC(i-1,k,0)+tr11*cr2+tr12*cr3; ti2=CC(i ,k,0)+tr11*ci2+tr12*ci3; tr3=CC(i-1,k,0)+tr12*cr2+tr11*cr3; ti3=CC(i ,k,0)+tr12*ci2+tr11*ci3; MULPM(tr5,tr4,cr5,cr4,ti11,ti12) MULPM(ti5,ti4,ci5,ci4,ti11,ti12) PM(CH(i-1,2,k),CH(ic-1,1,k),tr2,tr5) PM(CH(i ,2,k),CH(ic ,1,k),ti5,ti2) PM(CH(i-1,4,k),CH(ic-1,3,k),tr3,tr4) PM(CH(i ,4,k),CH(ic ,3,k),ti4,ti3) } } #undef CH #undef CC #define CH(a,b,c) ch[(a)+ido*((b)+l1*(c))] #define CC(a,b,c) cc[(a)+ido*((b)+cdim*(c))] #define C1(a,b,c) cc[(a)+ido*((b)+l1*(c))] #define C2(a,b) cc[(a)+idl1*(b)] #define CH2(a,b) ch[(a)+idl1*(b)] static void radfg(size_t ido, size_t ip, size_t l1, size_t idl1, FLOAT *cc, FLOAT *ch, const FLOAT *wa) { const size_t cdim=ip; static const FLOAT twopi=6.28318530717958647692; size_t idij, ipph, i, j, k, l, j2, ic, jc, lc, ik; FLOAT ai1, ai2, ar1, ar2, arg; FLOAT *csarr; size_t aidx; ipph=(ip+1)/ 2; if(ido!=1) { memcpy(ch,cc,idl1*sizeof(FLOAT)); for(j=1; j<ip; j++) for(k=0; k<l1; k++) { CH(0,k,j)=C1(0,k,j); idij=(j-1)*ido+1; for(i=2; i<ido; i+=2,idij+=2) MULPM(CH(i-1,k,j),CH(i,k,j),wa[idij-1],wa[idij],C1(i-1,k,j),C1(i,k,j)) } for(j=1,jc=ip-1; j<ipph; j++,jc--) for(k=0; k<l1; k++) for(i=2; i<ido; i+=2) { PM(C1(i-1,k,j),C1(i ,k,jc),CH(i-1,k,jc),CH(i-1,k,j )) PM(C1(i ,k,j),C1(i-1,k,jc),CH(i ,k,j ),CH(i ,k,jc)) } } else memcpy(cc,ch,idl1*sizeof(FLOAT)); for(j=1,jc=ip-1; j<ipph; j++,jc--) for(k=0; k<l1; k++) PM(C1(0,k,j),C1(0,k,jc),CH(0,k,jc),CH(0,k,j)) csarr=RALLOC(FLOAT,2*ip); arg=twopi / ip; csarr[0]=1.; csarr[1]=0.; csarr[2]=csarr[2*ip-2]=cos(arg); csarr[3]=sin(arg); csarr[2*ip-1]=-csarr[3]; for (i=2; i<=ip/2; ++i) { csarr[2*i]=csarr[2*ip-2*i]=cos(i*arg); csarr[2*i+1]=sin(i*arg); csarr[2*ip-2*i+1]=-csarr[2*i+1]; } for(l=1,lc=ip-1; l<ipph; l++,lc--) { ar1=csarr[2*l]; ai1=csarr[2*l+1]; for(ik=0; ik<idl1; ik++) { CH2(ik,l)=C2(ik,0)+ar1*C2(ik,1); CH2(ik,lc)=ai1*C2(ik,ip-1); } aidx=2*l; for(j=2,jc=ip-2; j<ipph; j++,jc--) { aidx+=2*l; if (aidx>=2*ip) aidx-=2*ip; ar2=csarr[aidx]; ai2=csarr[aidx+1]; for(ik=0; ik<idl1; ik++) { CH2(ik,l )+=ar2*C2(ik,j ); CH2(ik,lc)+=ai2*C2(ik,jc); } } } DEALLOC(csarr); for(j=1; j<ipph; j++) for(ik=0; ik<idl1; ik++) CH2(ik,0)+=C2(ik,j); for(k=0; k<l1; k++) memcpy(&CC(0,0,k),&CH(0,k,0),ido*sizeof(FLOAT)); for(j=1; j<ipph; j++) { jc=ip-j; j2=2*j; for(k=0; k<l1; k++) { CC(ido-1,j2-1,k) = CH(0,k,j ); CC(0 ,j2 ,k) = CH(0,k,jc); } } if(ido==1) return; for(j=1; j<ipph; j++) { jc=ip-j; j2=2*j; for(k=0; k<l1; k++) for(i=2; i<ido; i+=2) { ic=ido-i; PM (CC(i-1,j2,k),CC(ic-1,j2-1,k),CH(i-1,k,j ),CH(i-1,k,jc)) PM (CC(i ,j2,k),CC(ic ,j2-1,k),CH(i ,k,jc),CH(i ,k,j )) } } } #undef CC #undef CH #define CH(a,b,c) ch[(a)+ido*((b)+l1*(c))] #define CC(a,b,c) cc[(a)+ido*((b)+cdim*(c))] static void radb2(size_t ido, size_t l1, const FLOAT *cc, FLOAT *ch, const FLOAT *wa) { const size_t cdim=2; size_t i, k, ic; FLOAT ti2, tr2; for (k=0; k<l1; k++) PM (CH(0,k,0),CH(0,k,1),CC(0,0,k),CC(ido-1,1,k)) if ((ido&1)==0) for (k=0; k<l1; k++) { CH(ido-1,k,0) = 2*CC(ido-1,0,k); CH(ido-1,k,1) = -2*CC(0 ,1,k); } if (ido<=2) return; for (k=0; k<l1;++k) for (i=2; i<ido; i+=2) { ic=ido-i; PM (CH(i-1,k,0),tr2,CC(i-1,0,k),CC(ic-1,1,k)) PM (ti2,CH(i ,k,0),CC(i ,0,k),CC(ic ,1,k)) MULPM (CH(i,k,1),CH(i-1,k,1),WA(0,i-2),WA(0,i-1),ti2,tr2) } } static void radb3(size_t ido, size_t l1, const FLOAT *cc, FLOAT *ch, const FLOAT *wa) { const size_t cdim=3; static const FLOAT taur=-0.5, taui=0.86602540378443864676; size_t i, k, ic; FLOAT ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2; for (k=0; k<l1; k++) { tr2=2*CC(ido-1,1,k); cr2=CC(0,0,k)+taur*tr2; CH(0,k,0)=CC(0,0,k)+tr2; ci3=2*taui*CC(0,2,k); PM (CH(0,k,2),CH(0,k,1),cr2,ci3); } if (ido==1) return; for (k=0; k<l1; k++) for (i=2; i<ido; i+=2) { ic=ido-i; tr2=CC(i-1,2,k)+CC(ic-1,1,k); ti2=CC(i ,2,k)-CC(ic ,1,k); cr2=CC(i-1,0,k)+taur*tr2; ci2=CC(i ,0,k)+taur*ti2; CH(i-1,k,0)=CC(i-1,0,k)+tr2; CH(i ,k,0)=CC(i ,0,k)+ti2; cr3=taui*(CC(i-1,2,k)-CC(ic-1,1,k)); ci3=taui*(CC(i ,2,k)+CC(ic ,1,k)); PM(dr3,dr2,cr2,ci3) PM(di2,di3,ci2,cr3) MULPM(CH(i,k,1),CH(i-1,k,1),WA(0,i-2),WA(0,i-1),di2,dr2) MULPM(CH(i,k,2),CH(i-1,k,2),WA(1,i-2),WA(1,i-1),di3,dr3) } } static void radb4(size_t ido, size_t l1, const FLOAT *cc, FLOAT *ch, const FLOAT *wa) { const size_t cdim=4; static const FLOAT sqrt2=1.41421356237309504880; size_t i, k, ic; FLOAT ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4; for (k=0; k<l1; k++) { PM (tr2,tr1,CC(0,0,k),CC(ido-1,3,k)) tr3=2*CC(ido-1,1,k); tr4=2*CC(0,2,k); PM (CH(0,k,0),CH(0,k,2),tr2,tr3) PM (CH(0,k,3),CH(0,k,1),tr1,tr4) } if ((ido&1)==0) for (k=0; k<l1; k++) { PM (ti1,ti2,CC(0 ,3,k),CC(0 ,1,k)) PM (tr2,tr1,CC(ido-1,0,k),CC(ido-1,2,k)) CH(ido-1,k,0)=tr2+tr2; CH(ido-1,k,1)=sqrt2*(tr1-ti1); CH(ido-1,k,2)=ti2+ti2; CH(ido-1,k,3)=-sqrt2*(tr1+ti1); } if (ido<=2) return; for (k=0; k<l1;++k) for (i=2; i<ido; i+=2) { ic=ido-i; PM (tr2,tr1,CC(i-1,0,k),CC(ic-1,3,k)) PM (ti1,ti2,CC(i ,0,k),CC(ic ,3,k)) PM (tr4,ti3,CC(i ,2,k),CC(ic ,1,k)) PM (tr3,ti4,CC(i-1,2,k),CC(ic-1,1,k)) PM (CH(i-1,k,0),cr3,tr2,tr3) PM (CH(i ,k,0),ci3,ti2,ti3) PM (cr4,cr2,tr1,tr4) PM (ci2,ci4,ti1,ti4) MULPM (CH(i,k,1),CH(i-1,k,1),WA(0,i-2),WA(0,i-1),ci2,cr2) MULPM (CH(i,k,2),CH(i-1,k,2),WA(1,i-2),WA(1,i-1),ci3,cr3) MULPM (CH(i,k,3),CH(i-1,k,3),WA(2,i-2),WA(2,i-1),ci4,cr4) } } static void radb5(size_t ido, size_t l1, const FLOAT *cc, FLOAT *ch, const FLOAT *wa) { const size_t cdim=5; static const FLOAT tr11= 0.3090169943749474241, ti11=0.95105651629515357212, tr12=-0.8090169943749474241, ti12=0.58778525229247312917; size_t i, k, ic; FLOAT ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5; for (k=0; k<l1; k++) { ti5=2*CC(0,2,k); ti4=2*CC(0,4,k); tr2=2*CC(ido-1,1,k); tr3=2*CC(ido-1,3,k); CH(0,k,0)=CC(0,0,k)+tr2+tr3; cr2=CC(0,0,k)+tr11*tr2+tr12*tr3; cr3=CC(0,0,k)+tr12*tr2+tr11*tr3; MULPM(ci5,ci4,ti5,ti4,ti11,ti12) PM(CH(0,k,4),CH(0,k,1),cr2,ci5) PM(CH(0,k,3),CH(0,k,2),cr3,ci4) } if (ido==1) return; for (k=0; k<l1;++k) for (i=2; i<ido; i+=2) { ic=ido-i; PM(tr2,tr5,CC(i-1,2,k),CC(ic-1,1,k)) PM(ti5,ti2,CC(i ,2,k),CC(ic ,1,k)) PM(tr3,tr4,CC(i-1,4,k),CC(ic-1,3,k)) PM(ti4,ti3,CC(i ,4,k),CC(ic ,3,k)) CH(i-1,k,0)=CC(i-1,0,k)+tr2+tr3; CH(i ,k,0)=CC(i ,0,k)+ti2+ti3; cr2=CC(i-1,0,k)+tr11*tr2+tr12*tr3; ci2=CC(i ,0,k)+tr11*ti2+tr12*ti3; cr3=CC(i-1,0,k)+tr12*tr2+tr11*tr3; ci3=CC(i ,0,k)+tr12*ti2+tr11*ti3; MULPM(cr5,cr4,tr5,tr4,ti11,ti12) MULPM(ci5,ci4,ti5,ti4,ti11,ti12) PM(dr4,dr3,cr3,ci4) PM(di3,di4,ci3,cr4) PM(dr5,dr2,cr2,ci5) PM(di2,di5,ci2,cr5) MULPM(CH(i,k,1),CH(i-1,k,1),WA(0,i-2),WA(0,i-1),di2,dr2) MULPM(CH(i,k,2),CH(i-1,k,2),WA(1,i-2),WA(1,i-1),di3,dr3) MULPM(CH(i,k,3),CH(i-1,k,3),WA(2,i-2),WA(2,i-1),di4,dr4) MULPM(CH(i,k,4),CH(i-1,k,4),WA(3,i-2),WA(3,i-1),di5,dr5) } } static void radbg(size_t ido, size_t ip, size_t l1, size_t idl1, FLOAT *cc, FLOAT *ch, const FLOAT *wa) { const size_t cdim=ip; static const FLOAT twopi=6.28318530717958647692; size_t idij, ipph, i, j, k, l, j2, ic, jc, lc, ik; FLOAT ai1, ai2, ar1, ar2, arg; FLOAT *csarr; size_t aidx; ipph=(ip+1)/ 2; for(k=0; k<l1; k++) memcpy(&CH(0,k,0),&CC(0,0,k),ido*sizeof(FLOAT)); for(j=1; j<ipph; j++) { jc=ip-j; j2=2*j; for(k=0; k<l1; k++) { CH(0,k,j )=2*CC(ido-1,j2-1,k); CH(0,k,jc)=2*CC(0 ,j2 ,k); } } if(ido!=1) for(j=1,jc=ip-1; j<ipph; j++,jc--) for(k=0; k<l1; k++) for(i=2; i<ido; i+=2) { ic=ido-i; PM (CH(i-1,k,j ),CH(i-1,k,jc),CC(i-1,2*j,k),CC(ic-1,2*j-1,k)) PM (CH(i ,k,jc),CH(i ,k,j ),CC(i ,2*j,k),CC(ic ,2*j-1,k)) } csarr=RALLOC(FLOAT,2*ip); arg=twopi/ip; csarr[0]=1.; csarr[1]=0.; csarr[2]=csarr[2*ip-2]=cos(arg); csarr[3]=sin(arg); csarr[2*ip-1]=-csarr[3]; for (i=2; i<=ip/2; ++i) { csarr[2*i]=csarr[2*ip-2*i]=cos(i*arg); csarr[2*i+1]=sin(i*arg); csarr[2*ip-2*i+1]=-csarr[2*i+1]; } for(l=1; l<ipph; l++) { lc=ip-l; ar1=csarr[2*l]; ai1=csarr[2*l+1]; for(ik=0; ik<idl1; ik++) { C2(ik,l)=CH2(ik,0)+ar1*CH2(ik,1); C2(ik,lc)=ai1*CH2(ik,ip-1); } aidx=2*l; for(j=2; j<ipph; j++) { jc=ip-j; aidx+=2*l; if (aidx>=2*ip) aidx-=2*ip; ar2=csarr[aidx]; ai2=csarr[aidx+1]; for(ik=0; ik<idl1; ik++) { C2(ik,l )+=ar2*CH2(ik,j ); C2(ik,lc)+=ai2*CH2(ik,jc); } } } DEALLOC(csarr); for(j=1; j<ipph; j++) for(ik=0; ik<idl1; ik++) CH2(ik,0)+=CH2(ik,j); for(j=1,jc=ip-1; j<ipph; j++,jc--) for(k=0; k<l1; k++) PM (CH(0,k,jc),CH(0,k,j),C1(0,k,j),C1(0,k,jc)) if(ido==1) return; for(j=1,jc=ip-1; j<ipph; j++,jc--) for(k=0; k<l1; k++) for(i=2; i<ido; i+=2) { PM (CH(i-1,k,jc),CH(i-1,k,j ),C1(i-1,k,j),C1(i ,k,jc)) PM (CH(i ,k,j ),CH(i ,k,jc),C1(i ,k,j),C1(i-1,k,jc)) } memcpy(cc,ch,idl1*sizeof(FLOAT)); for(j=1; j<ip; j++) for(k=0; k<l1; k++) { C1(0,k,j)=CH(0,k,j); idij=(j-1)*ido+1; for(i=2; i<ido; i+=2,idij+=2) MULPM (C1(i,k,j),C1(i-1,k,j),wa[idij-1],wa[idij],CH(i,k,j),CH(i-1,k,j)) } } #undef CC #undef CH #undef PM #undef MULPM /*---------------------------------------------------------------------- cfftf1, cfftb1, cfftf, cfftb, cffti1, cffti. Complex FFTs. ----------------------------------------------------------------------*/ static void cfft1(size_t n, cmplx c[], cmplx ch[], const cmplx wa[], const size_t ifac[], int isign) { size_t k1, l1=1, nf=ifac[1], iw=0; cmplx *p1=c, *p2=ch; for(k1=0; k1<nf; k1++) { size_t ip=ifac[k1+2]; size_t l2=ip*l1; size_t ido = n/l2; if(ip==4) (isign>0) ? passb4(ido, l1, p1, p2, wa+iw) : passf4(ido, l1, p1, p2, wa+iw); else if(ip==2) (isign>0) ? passb2(ido, l1, p1, p2, wa+iw) : passf2(ido, l1, p1, p2, wa+iw); else if(ip==3) (isign>0) ? passb3(ido, l1, p1, p2, wa+iw) : passf3(ido, l1, p1, p2, wa+iw); else if(ip==5) (isign>0) ? passb5(ido, l1, p1, p2, wa+iw) : passf5(ido, l1, p1, p2, wa+iw); else if(ip==6) (isign>0) ? passb6(ido, l1, p1, p2, wa+iw) : passf6(ido, l1, p1, p2, wa+iw); else (isign>0) ? passbg(ido, ip, l1, p1, p2, wa+iw) : passfg(ido, ip, l1, p1, p2, wa+iw); SWAP(p1,p2,cmplx *); l1=l2; iw+=(ip-1)*ido; } if (p1!=c) memcpy (c,p1,n*sizeof(cmplx)); } void cfftf(size_t n, FLOAT c[], FLOAT wsave[]) { if (n!=1) cfft1(n, (cmplx*)c, (cmplx*)wsave, (cmplx*)(wsave+2*n), (size_t*)(wsave+4*n),-1); } void cfftb(size_t n, FLOAT c[], FLOAT wsave[]) { if (n!=1) cfft1(n, (cmplx*)c, (cmplx*)wsave, (cmplx*)(wsave+2*n), (size_t*)(wsave+4*n),+1); } static void factorize (size_t n, const size_t *pf, size_t npf, size_t *ifac) { size_t nl=n, nf=0, ntry=0, j=0, i; startloop: j++; ntry = (j<=npf) ? pf[j-1] : ntry+2; do { size_t nq=nl / ntry; size_t nr=nl-ntry*nq; if (nr!=0) goto startloop; nf++; ifac[nf+1]=ntry; nl=nq; if ((ntry==2) && (nf!=1)) { for (i=nf+1; i>2; --i) ifac[i]=ifac[i-1]; ifac[2]=2; } } while(nl!=1); ifac[0]=n; ifac[1]=nf; } static void cffti1(size_t n, FLOAT wa[], size_t ifac[]) { static const size_t ntryh[5]={4,6,3,2,5}; static const FLOAT twopi=6.28318530717958647692; size_t j, k, fi; FLOAT argh=twopi/n; size_t i=0, l1=1; factorize (n,ntryh,5,ifac); for(k=1; k<=ifac[1]; k++) { size_t ip=ifac[k+1]; size_t ido=n/(l1*ip); for(j=1; j<ip; j++) { size_t is = i; FLOAT argld=j*l1*argh; wa[i ]=1; wa[i+1]=0; for(fi=1; fi<=ido; fi++) { FLOAT arg=fi*argld; i+=2; wa[i ]=cos(arg); wa[i+1]=sin(arg); } if(ip>6) { wa[is ]=wa[i ]; wa[is+1]=wa[i+1]; } } l1*=ip; } } void cffti(size_t n, FLOAT wsave[]) { if (n!=1) cffti1(n, wsave+2*n,(size_t*)(wsave+4*n)); } /*---------------------------------------------------------------------- rfftf1, rfftb1, rfftf, rfftb, rffti1, rffti. Real FFTs. ----------------------------------------------------------------------*/ static void rfftf1(size_t n, FLOAT c[], FLOAT ch[], const FLOAT wa[], const size_t ifac[]) { size_t k1, l1=n, nf=ifac[1], iw=n-1; FLOAT *p1=ch, *p2=c; for(k1=1; k1<=nf;++k1) { size_t ip=ifac[nf-k1+2]; size_t ido=n / l1; l1 /= ip; iw-=(ip-1)*ido; SWAP (p1,p2,FLOAT *); if(ip==4) radf4(ido, l1, p1, p2, wa+iw); else if(ip==2) radf2(ido, l1, p1, p2, wa+iw); else if(ip==3) radf3(ido, l1, p1, p2, wa+iw); else if(ip==5) radf5(ido, l1, p1, p2, wa+iw); else { if (ido==1) SWAP (p1,p2,FLOAT *); radfg(ido, ip, l1, ido*l1, p1, p2, wa+iw); SWAP (p1,p2,FLOAT *); } } if (p1==c) memcpy (c,ch,n*sizeof(FLOAT)); } static void rfftb1(size_t n, FLOAT c[], FLOAT ch[], const FLOAT wa[], const size_t ifac[]) { size_t k1, l1=1, nf=ifac[1], iw=0; FLOAT *p1=c, *p2=ch; for(k1=1; k1<=nf; k1++) { size_t ip = ifac[k1+1], ido= n/(ip*l1); if(ip==4) radb4(ido, l1, p1, p2, wa+iw); else if(ip==2) radb2(ido, l1, p1, p2, wa+iw); else if(ip==3) radb3(ido, l1, p1, p2, wa+iw); else if(ip==5) radb5(ido, l1, p1, p2, wa+iw); else { radbg(ido, ip, l1, ido*l1, p1, p2, wa+iw); if (ido!=1) SWAP (p1,p2,FLOAT *); } SWAP (p1,p2,FLOAT *); l1*=ip; iw+=(ip-1)*ido; } if (p1!=c) memcpy (c,ch,n*sizeof(FLOAT)); } void rfftf(size_t n, FLOAT r[], FLOAT wsave[]) { if(n!=1) rfftf1(n, r, wsave, wsave+n,(size_t*)(wsave+2*n)); } void rfftb(size_t n, FLOAT r[], FLOAT wsave[]) { if(n!=1) rfftb1(n, r, wsave, wsave+n,(size_t*)(wsave+2*n)); } static void rffti1(size_t n, FLOAT wa[], size_t ifac[]) { static const size_t ntryh[4]={4,2,3,5}; static const FLOAT twopi=6.28318530717958647692; size_t i, j, k, fi; FLOAT argh=twopi/n; size_t is=0, l1=1; factorize (n,ntryh,4,ifac); for (k=1; k<ifac[1]; k++) { size_t ip=ifac[k+1], ido=n/(l1*ip); for (j=1; j<ip; ++j) { FLOAT argld=j*l1*argh; for(i=is,fi=1; i<=ido+is-3; i+=2,++fi) { FLOAT arg=fi*argld; wa[i ]=cos(arg); wa[i+1]=sin(arg); } is+=ido; } l1*=ip; } } void rffti(size_t n, FLOAT wsave[]) { if (n!=1) rffti1(n, wsave+n,(size_t*)(wsave+2*n)); }