BRL-CAD
fftfast.c
Go to the documentation of this file.
1 /* F F T F A S T . C
2  * BRL-CAD
3  *
4  * Copyright (c) 2004-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 /** @file libfft/fftfast.c
21  *
22  * Complex Number and FFT Library
23  *
24  * "Fast" Version - Function calls to complex math routines removed.
25  * Uses pre-computed sine/cosine tables.
26  *
27  * The FFT is:
28  *
29  * N-1
30  * Xf(k) = Sum x(n)(cos(2PI(nk/N)) - isin(2PI(nk/N)))
31  * n=0
32  *
33  */
34 
35 #include "common.h"
36 
37 #include <stdlib.h>
38 #include <stdio.h> /* for stderr */
39 
40 #include "fft.h"
41 
42 
43 #define MAXSIZE 65536 /* Needed for sin/cos tables */
44 int _init_size = 0; /* Internal: shows last initialized size */
45 
46 void scramble(int numpoints, COMPLEX *dat);
47 void butterflies(int numpoints, int inverse, COMPLEX *dat);
48 int init_sintab(int size);
49 
50 
51 /*
52  * Forward Complex Fourier Transform
53  */
54 void
55 cfft(COMPLEX *dat, int num)
56 {
57  /* Check for trig table initialization */
58  if (num != _init_size) {
59  if (init_sintab(num) == 0) {
60  /* Can't do requested size */
61  return;
62  }
63  }
64 
65  scramble(num, dat);
66  butterflies(num, -1, dat);
67 }
68 
69 /*
70  * Inverse Complex Fourier Transform
71  */
72 void
73 icfft(COMPLEX *dat, int num)
74 {
75  /* Check for trig table initialization */
76  if (num != _init_size) {
77  if (init_sintab(num) == 0) {
78  /* Can't do requested size */
79  return;
80  }
81  }
82 
83  scramble(num, dat);
84  butterflies(num, 1, dat);
85 }
86 
87 /******************* INTERNAL FFT ROUTINES ********************/
88 
89 /* The trig tables */
90 double *sintab;
91 double *costab;
92 
93 /*
94  * Internal routine to initialize the sine/cosine table for
95  * transforms of a given size. Checks size for power of two
96  * and within table limits.
97  *
98  * Note that once initialized for one size it ready for one
99  * smaller than that also, but it is convenient to do power of
100  * two checking here so we change the _init_size every time
101  * (We *could* pick up where ever we left off by keeping a
102  * _max_init_size but forget that for now).
103  *
104  * Note that we need sin and cos values for +/- (M_PI * (m / col))
105  * where col = 1, 2, 4, ..., N/2
106  * m = 0, 1, 2, ..., col-1
107  *
108  * Thus we can subscript by: table[(m / col) * N/2]
109  * or with twice as many values by: table[m + col]
110  * We chose the later. (but N.B. this doesn't allow sub
111  * _init_size requests to use existing numbers!)
112  */
113 int
114 init_sintab(int size)
115 {
116  double theta;
117  int col, m;
118 
119  /*
120  * Check whether the requested size is within our compiled
121  * limit and make sure it's a power of two.
122  */
123  if (size > MAXSIZE) {
124  fprintf(stderr, "fft: Only compiled for max size of %d\n", MAXSIZE);
125  fprintf(stderr, "fft: Can't do the requested %d\n", size);
126  return 0;
127  }
128  for (m = size; (m & 1) == 0; m >>= 1)
129  ;
130  if (m != 1) {
131  fprintf(stderr, "fft: Can only do powers of two, not %d\n", size);
132  fprintf(stderr, "fft: What do you think this is, a Winograd transform?\n");
133  return 0;
134  }
135 
136  /* Get some buffer space */
137  if (sintab != NULL) free(sintab);
138  if (costab != NULL) free(costab);
139  /* should not use bu_calloc() as libfft is not dependent upon libbu */
140  sintab = (double *)calloc(sizeof(*sintab), size);
141  costab = (double *)calloc(sizeof(*costab), size);
142 
143  /*
144  * Size is okay. Set up tables.
145  */
146  for (col = 1; col < size; col <<= 1) {
147  for (m = 0; m < col; m++) {
148  theta = M_PI * (double)m / (double)col;
149  sintab[ m + col ] = sin(theta);
150  costab[ m + col ] = cos(theta);
151  }
152  }
153 
154  /*
155  * Mark size and return success.
156  */
157  _init_size = size;
158 /* fprintf(stderr, "fft: table init, size = %d\n", size);*/
159  return 1;
160 }
161 
162 /*
163  * This section implements the Cooley-Tukey Complex
164  * Fourier transform.
165  * Reference: Nov 84 Dr. Dobbs [#97], and
166  * "Musical Applications of Microprocessors", Hal Chamberlin.
167  */
168 
169 /*
170  * SCRAMBLE - put data in bit reversed order
171  *
172  * Note: Could speed this up with pointers if necessary,
173  * but the butterflies take much longer.
174  */
175 void
176 scramble(int numpoints, COMPLEX *dat)
177 {
178  register int i, j, m;
179  COMPLEX temp;
180 
181  j = 0;
182  for (i = 0; i < numpoints; i++, j += m) {
183  if (i < j) {
184  /* Switch nodes i and j */
185  temp.re = dat[j].re;
186  temp.im = dat[j].im;
187  dat[j].re = dat[i].re;
188  dat[j].im = dat[i].im;
189  dat[i].re = temp.re;
190  dat[i].im = temp.im;
191  }
192  m = numpoints/2;
193  while (m-1 < j) {
194  j -= m;
195  m = (m + 1) / 2;
196  }
197  }
198 }
199 
200 void
201 butterflies(int numpoints, int inverse, COMPLEX *dat)
202 {
203  register COMPLEX *node1, *node2;
204  register int step, column, m;
205  COMPLEX w, temp;
206 
207  /*
208  * For each column of the butterfly
209  */
210  for (column = 1; column < numpoints; column = step) {
211  step = 2 * column; /* step is size of "cross-hatch" */
212  /*
213  * For each principle value of W (roots on units).
214  */
215  for (m = 0; m < column; m++) {
216  /*
217  * Do these by table lookup:
218  * theta = M_PI*(inverse*m)/column;
219  * w.re = cos(theta);
220  * w.im = sin(theta);
221  */
222  w.re = costab[ column + m ];
223  w.im = sintab[ column + m ] * inverse;
224  /* Do all pairs of nodes */
225  for (node1 = &dat[m]; node1 < &dat[numpoints]; node1 += step) {
226  node2 = node1 + column;
227  /*
228  * Want to compute:
229  * dat[node2] = dat[node1] - w * dat[node2];
230  * dat[node1] = dat[node1] + w * dat[node2];
231  *
232  * We do all this with pointers now.
233  */
234 
235  /*cmult(&w, &dat[node2], &temp);*/
236  temp.re = w.re*node2->re - w.im*node2->im;
237  temp.im = w.im*node2->re + w.re*node2->im;
238 
239  /*csub(&dat[node1], &temp, &dat[node2]);*/
240  node2->re = node1->re - temp.re;
241  node2->im = node1->im - temp.im;
242 
243  /*cadd(&dat[node1], &temp, &dat[node1]);*/
244  node1->re += temp.re;
245  node1->im += temp.im;
246  }
247  }
248  }
249 
250  /* Scale Data (on forward transform only) */
251  /*
252  * Technically speaking this gives us the periodogram. XXX
253  * The canonical definition does the scaling only
254  * after the inverse xform. Our method may hurt certain
255  * other forms of analysis, e.g. cepstrum.
256  * **** We Now Do It The Canonical Way! ****
257  */
258  if (inverse > 0) {
259  for (node1 = &dat[0]; node1 < &dat[numpoints]; node1++) {
260  /* cdiv(&dat[i], &const, &dat[i]); */
261  node1->re /= (double)numpoints;
262  node1->im /= (double)numpoints;
263  }
264  }
265 }
266 
267 /**** COMPLEX ARITHMETIC ROUTINES ****/
268 /**** NO LONGER USED BY TRANSFORMS ***/
269 /*
270  * CADD - Complex ADDition
271  */
272 void
273 cadd(COMPLEX *result, COMPLEX *val1, COMPLEX *val2)
274 {
275  result->re = val1->re + val2->re;
276  result->im = val1->im + val2->im;
277 }
278 
279 /*
280  * CSUB - Complex SUBtraction
281  */
282 void
283 csub(COMPLEX *result, COMPLEX *val1, COMPLEX *val2)
284 {
285  result->re = val1->re - val2->re;
286  result->im = val1->im - val2->im;
287 }
288 
289 /*
290  * CMULT - Complex MULTiply
291  */
292 void
293 cmult(COMPLEX *result, COMPLEX *val1, COMPLEX *val2)
294 {
295  result->re = val1->re*val2->re - val1->im*val2->im;
296  result->im = val1->im*val2->re + val1->re*val2->im;
297 }
298 
299 /*
300  * CDIV - Complex DIVide
301  */
302 void
303 cdiv(COMPLEX *result, COMPLEX *val1, COMPLEX *val2)
304 {
305  double denom;
306 
307  denom = val2->re*val2->re + val2->im*val2->im;
308  result->re = (val1->re*val2->re + val1->im*val2->im)/denom;
309  result->im = (val1->im*val2->re - val1->re*val2->im)/denom;
310 }
311 
312 /*
313  * Local Variables:
314  * mode: C
315  * tab-width: 8
316  * indent-tabs-mode: t
317  * c-file-style: "stroustrup"
318  * End:
319  * ex: shiftwidth=4 tabstop=8
320  */
double im
Definition: fft.h:56
void butterflies(int numpoints, int inverse, COMPLEX *dat)
Definition: fftfast.c:201
void cfft(COMPLEX *dat, int num)
Definition: fftfast.c:55
#define M_PI
Definition: fft.h:35
void cdiv(COMPLEX *result, COMPLEX *val1, COMPLEX *val2)
Definition: fftfast.c:303
Header file for the BRL-CAD common definitions.
void cmult(COMPLEX *result, COMPLEX *val1, COMPLEX *val2)
Definition: fftfast.c:293
int _init_size
Definition: fftfast.c:44
Definition: fft.h:54
#define MAXSIZE
Definition: fftfast.c:43
int init_sintab(int size)
Definition: fftfast.c:114
void icfft(COMPLEX *dat, int num)
Definition: fftfast.c:73
double re
Definition: fft.h:55
void csub(COMPLEX *result, COMPLEX *val1, COMPLEX *val2)
Definition: fftfast.c:283
double * sintab
Definition: fftfast.c:90
void scramble(int numpoints, COMPLEX *dat)
Definition: fftfast.c:176
double * costab
Definition: fftfast.c:91
void cadd(COMPLEX *result, COMPLEX *val1, COMPLEX *val2)
Definition: fftfast.c:273