1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include <math.h>
8 #include <stdlib.h>
9 #include <stdio.h>
10
11 #ifndef PI
12 # define PI 3.1415926536
13 #endif
14
15 /*
16 mode : -1 forward, 1 inverse.
17 idim : FFT dimension, max 1024.
18 xr : real array.
19 xi : imaginary part of the array.
20 */
21
22 /* --------------------------------- */
cfft(int mode,int idim,float * xr,float * xi)23 void cfft( int mode , int idim , float *xr , float *xi )
24 /* --------------------------------- */
25 {
26 #define IDMAX 1024
27 #define LIDMAX 10
28 static int idold = -999 ;
29 register int i0,i1,i2,i4,i5=0,m0, id = idim;
30 float f1,f3,f4,f5,al,co,si,md = mode;
31 static int n,m[LIDMAX];
32 static float f2,c[IDMAX/2],s[IDMAX/2];
33
34 /**** preparations if id has changed since last call ****/
35
36 if (idold != id) {
37 idold = id ;
38 /* check for power of 2 */
39 for( i4=4 ; i4 <= IDMAX ; i4 *= 2 ){
40 if( id == i4 )break ;
41 }
42 if( id != i4 ){
43 fprintf(stderr,"\n In cfft : illegal idim=%d\n",idim);
44 exit(1) ;
45 }
46 f2 = id;
47 n = log(f2)/log(2.) + .5;
48 m[n-1] = 1;
49 al = 2.*PI/f2;
50 co = cos(al);
51 si = sin(al);
52 c[0] = 1.;
53 s[0] = 0.;
54 for(i4 = 1; i4 < 512; i4++) {
55 c[i4] = c[i4-1]*co - s[i4-1]*si;
56 s[i4] = s[i4-1]*co + c[i4-1]*si;
57 }
58 for(i4 = n-1; i4 >= 1; i4--) m[i4-1] = m[i4]*2;
59 }
60
61 /**** Main loop starts here ****/
62
63 for(i0 = 0; i0 < n; i0++) {
64 i1 = 0;
65 m0 = m[i0];
66 for(i2 = 0; i2 < m[n-i0-1]; i2++) {
67 f4 = c[i1];
68 f5 = s[i1]*md;
69 for(i4 = 2*m0*i2; i4 < m0*(2*i2+1); i4++) {
70 f3 = xr[i4+m0]*f4 - xi[i4+m0]*f5;
71 f1 = xi[i4+m0]*f4 + xr[i4+m0]*f5;
72 xr[i4+m0] = xr[i4] - f3;
73 xr[i4] += f3;
74 xi[i4+m0] = xi[i4] - f1;
75 xi[i4] += f1;
76 }
77 for(i4 = 1; i4 < n; i4++) {
78 i5 = i4;
79 if (i1 < m[i4]) goto i1_plus1;
80 i1 -= m[i4];
81 }
82 i1_plus1: i1 += m[i5];
83 }
84 }
85 i1 = 0;
86 for(i4 = 0; i4 < id; i4++) {
87 if (i1 > i4) {
88 f3 = xr[i4];
89 f1 = xi[i4];
90 xr[i4] = xr[i1];
91 xi[i4] = xi[i1];
92 xr[i1] = f3;
93 xi[i1] = f1;
94 }
95 for(i2 = 0; i2 < n; i2++) {
96 i5 = i2;
97 if (i1 < m[i2]) goto i1_plus2;
98 i1 -= m[i2];
99 }
100 i1_plus2: i1 += m[i5];
101 }
102
103 if (md > 0.) {
104 f1 = 1./f2;
105 for(i4 = 0; i4 < id; i4++) {
106 xr[i4] *= f1;
107 xi[i4] *= f1;
108 }
109 }
110 return;
111 }
112
113 /*********************************************************************/
114
115 /*----------------------------------*/
cfft2d_cox(int mode,int nx,int ny,float * xr,float * xi)116 void cfft2d_cox( int mode , int nx,int ny , float *xr, float *xi )
117 /*----------------------------------*/
118 {
119 float *rbuf , *ibuf ;
120 register int ii , jj , jbase ;
121
122 rbuf = (float *)malloc( ny * sizeof(float) ) ;
123 ibuf = (float *)malloc( ny * sizeof(float) ) ;
124 if( rbuf == NULL || ibuf == NULL ){
125 fprintf(stderr,"malloc error in cfft2d_cox\n") ;
126 exit(1) ;
127 }
128
129 for( jj=0 ; jj < ny ; jj++ ){
130 jbase = nx * jj ;
131 cfft( mode , nx , &xr[jbase] , &xi[jbase] ) ;
132 }
133
134 for( ii=0 ; ii < nx ; ii++ ){
135 for( jj=0 ; jj < ny ; jj++ ){
136 rbuf[jj] = xr[ii + jj*nx] ;
137 ibuf[jj] = xi[ii + jj*nx] ;
138 }
139 cfft( mode , ny , rbuf , ibuf ) ;
140 for( jj=0 ; jj < ny ; jj++ ){
141 xr[ii+jj*nx] = rbuf[jj] ;
142 xi[ii+jj*nx] = ibuf[jj] ;
143 }
144 }
145
146 free(rbuf) ; free(ibuf) ;
147 return ;
148 }
149