1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #ifdef WIN32
5 #include <sys\time.h>
6 #else
7 #include <sys/time.h>
8 #endif
9
10 #define PI 3.14159265359
11 #define MAXPOW 24
12
13 struct complex
14 {
15 double r;
16 double i;
17 };
18
19 int pow_2[MAXPOW];
20 int pow_4[MAXPOW];
21
twiddle(struct complex * W,int N,double stuff)22 void twiddle(struct complex *W, int N, double stuff)
23 {
24 W->r=cos(stuff*2.0*PI/(double)N);
25 W->i=-sin(stuff*2.0*PI/(double)N);
26 }
27
bit_reverse_reorder(struct complex * W,int N)28 void bit_reverse_reorder(struct complex *W, int N)
29 {
30 int bits, i, j, k;
31 double tempr, tempi;
32
33 for (i=0; i<MAXPOW; i++)
34 if (pow_2[i]==N) bits=i;
35
36 for (i=0; i<N; i++)
37 {
38 j=0;
39 for (k=0; k<bits; k++)
40 if (i&pow_2[k]) j+=pow_2[bits-k-1];
41
42 if (j>i) /** Only make "up" swaps */
43 {
44 tempr=W[i].r;
45 tempi=W[i].i;
46 W[i].r=W[j].r;
47 W[i].i=W[j].i;
48 W[j].r=tempr;
49 W[j].i=tempi;
50 }
51 }
52 }
bit_r4_reorder(struct complex * W,int N)53 void bit_r4_reorder(struct complex *W, int N)
54 {
55 int bits, i, j, k;
56 double tempr, tempi;
57
58 for (i=0; i<MAXPOW; i++)
59 if (pow_2[i]==N) bits=i;
60
61 for (i=0; i<N; i++)
62 {
63 j=0;
64 for (k=0; k<bits; k+=2)
65 {
66 if (i&pow_2[k]) j+=pow_2[bits-k-2];
67 if (i&pow_2[k+1]) j+=pow_2[bits-k-1];
68 }
69
70 if (j>i) /** Only make "up" swaps */
71 {
72 tempr=W[i].r;
73 tempi=W[i].i;
74 W[i].r=W[j].r;
75 W[i].i=W[j].i;
76 W[j].r=tempr;
77 W[j].i=tempi;
78 }
79 }
80 }
81
82 /** RADIX-4 FFT ALGORITHM */
radix4(struct complex * x,int N)83 void radix4(struct complex *x, int N)
84 {
85 int n2, k1, N1, N2;
86 struct complex W, bfly[4];
87
88 N1=4;
89 N2=N/4;
90
91 /** Do 4 Point DFT */
92 for (n2=0; n2<N2; n2++)
93 {
94 /** Don't hurt the butterfly */
95 bfly[0].r = (x[n2].r + x[N2 + n2].r + x[2*N2+n2].r + x[3*N2+n2].r);
96 bfly[0].i = (x[n2].i + x[N2 + n2].i + x[2*N2+n2].i + x[3*N2+n2].i);
97
98 bfly[1].r = (x[n2].r + x[N2 + n2].i - x[2*N2+n2].r - x[3*N2+n2].i);
99 bfly[1].i = (x[n2].i - x[N2 + n2].r - x[2*N2+n2].i + x[3*N2+n2].r);
100
101 bfly[2].r = (x[n2].r - x[N2 + n2].r + x[2*N2+n2].r - x[3*N2+n2].r);
102 bfly[2].i = (x[n2].i - x[N2 + n2].i + x[2*N2+n2].i - x[3*N2+n2].i);
103
104 bfly[3].r = (x[n2].r - x[N2 + n2].i - x[2*N2+n2].r + x[3*N2+n2].i);
105 bfly[3].i = (x[n2].i + x[N2 + n2].r - x[2*N2+n2].i - x[3*N2+n2].r);
106
107
108 /** In-place results */
109 for (k1=0; k1<N1; k1++)
110 {
111 twiddle(&W, N, (double)k1*(double)n2);
112 x[n2 + N2*k1].r = bfly[k1].r*W.r - bfly[k1].i*W.i;
113 x[n2 + N2*k1].i = bfly[k1].i*W.r + bfly[k1].r*W.i;
114 }
115 }
116
117 /** Don't recurse if we're down to one butterfly */
118 if (N2!=1)
119 for (k1=0; k1<N1; k1++)
120 {
121 radix4(&x[N2*k1], N2);
122 }
123 }
124
125 /** RADIX-2 FFT ALGORITHM */
radix2(struct complex * data,int N)126 void radix2(struct complex *data, int N)
127 {
128 int n2, k1, N1, N2;
129 struct complex W, bfly[2];
130
131 N1=2;
132 N2=N/2;
133
134 /** Do 2 Point DFT */
135 for (n2=0; n2<N2; n2++)
136 {
137 /** Don't hurt the butterfly */
138 twiddle(&W, N, (double)n2);
139 bfly[0].r = (data[n2].r + data[N2 + n2].r);
140 bfly[0].i = (data[n2].i + data[N2 + n2].i);
141 bfly[1].r = (data[n2].r - data[N2 + n2].r) * W.r -
142 ((data[n2].i - data[N2 + n2].i) * W.i);
143 bfly[1].i = (data[n2].i - data[N2 + n2].i) * W.r +
144 ((data[n2].r - data[N2 + n2].r) * W.i);
145
146 /** In-place results */
147 for (k1=0; k1<N1; k1++)
148 {
149 data[n2 + N2*k1].r = bfly[k1].r;
150 data[n2 + N2*k1].i = bfly[k1].i;
151 }
152 }
153
154 /** Don't recurse if we're down to one butterfly */
155 if (N2!=1)
156 for (k1=0; k1<N1; k1++)
157 radix2(&data[N2*k1], N2);
158 }
159
main(int argc,char * argv[])160 void main(int argc, char *argv[])
161 {
162 FILE *infile;
163 int N, radix, numsamp;
164 int i;
165 struct complex *data;
166 double freq, phase, fs, A;
167 int dotime;
168 struct timeval start, end;
169 long totaltime;
170
171 #ifdef GEN
172 if (argc!=9)
173 {
174 printf("usage:\n");
175 printf(" fft [A] [f] [phase] [fs] [num samp] [sequence length] [radix] [time]\n");
176 printf(" output: DFT\n");
177 exit(1);
178 }
179
180
181 sscanf(argv[1], "%lf", &A);
182 sscanf(argv[2], "%lf", &freq);
183 sscanf(argv[3], "%lf", &phase);
184 sscanf(argv[4], "%lf", &fs);
185 sscanf(argv[5], "%d", &numsamp);
186 sscanf(argv[6], "%d", &N);
187 sscanf(argv[7], "%d", &radix);
188 sscanf(argv[8], "%d", &dotime);
189 #endif
190 #ifndef GEN
191 if (argc<4)
192 {
193 printf("usage:\n");
194 printf(" fft [input file] [sequence length] [radix]\n");
195 printf(" output: DFT\n");
196 exit(1);
197 }
198 else if ((infile=fopen(argv[1], "r"))==NULL)
199 {
200 printf("Error reading input sequence file: %s\n", argv[1]);
201 exit(1);
202 }
203
204 sscanf(argv[2], "%d", &N);
205 sscanf(argv[3], "%d", &radix);
206 dotime=0;
207 #endif
208
209
210 /** Set up power of two arrays */
211 pow_2[0]=1;
212 for (i=1; i<MAXPOW; i++)
213 pow_2[i]=pow_2[i-1]*2;
214 pow_4[0]=1;
215 for (i=1; i<MAXPOW; i++)
216 pow_4[i]=pow_4[i-1]*4;
217
218 if ((data=malloc(sizeof(struct complex)*(size_t)N))==NULL)
219 {
220 fprintf(stderr, "Out of memory!\n");
221 exit(1);
222 }
223
224 /** Generate cosine **/
225 #ifdef GEN
226 for (i=0; i<N; i++)
227 {
228 data[i].r=0.0;
229 data[i].i=0.0;
230 }
231 for (i=0; i<numsamp; i++)
232 data[i].r=A*cos(2.0*PI*freq*i/fs - phase*PI/180);
233 #endif
234 #ifndef GEN
235 for (i=0; i<N; i++)
236 {
237 fscanf(infile, "%lf", &data[i].r);
238 data[i].i=0.0;
239 }
240 #endif
241
242 gettimeofday(&start, (struct timezone *) 0);
243 if (radix==2) radix2(data, N);
244 if (radix==4) radix4(data, N);
245 gettimeofday(&end, (struct timezone *) 0);
246 totaltime=(end.tv_sec*1000000 + end.tv_usec)-(start.tv_sec*1000000
247 + start.tv_usec);
248 if (radix==2) bit_reverse_reorder(data, N);
249 if (radix==4) bit_r4_reorder(data, N);
250
251 if (!dotime)
252 for (i=0; i<N; i++)
253 printf("%f\n", sqrt(data[i].r*data[i].r +
254 data[i].i*data[i].i));
255 else
256 printf("%ld us\n\n", totaltime);
257
258 #ifndef GEN
259 fclose(infile);
260 #endif
261 }
262
263
264
265