summaryrefslogtreecommitdiff
path: root/apps/plugins/pdbox/PDa/src/d_fftroutine.c
diff options
context:
space:
mode:
authorPeter D'Hoye <peter.dhoye@gmail.com>2009-05-22 21:58:48 +0000
committerPeter D'Hoye <peter.dhoye@gmail.com>2009-05-22 21:58:48 +0000
commit513389b4c1bc8afe4b2dc9947c534bfeb105e3da (patch)
tree10e673b35651ac567fed2eda0c679c7ade64cbc6 /apps/plugins/pdbox/PDa/src/d_fftroutine.c
parent95fa7f6a2ef466444fbe3fe87efc6d5db6b77b36 (diff)
downloadrockbox-513389b4c1bc8afe4b2dc9947c534bfeb105e3da.tar.gz
rockbox-513389b4c1bc8afe4b2dc9947c534bfeb105e3da.zip
Add FS #10214. Initial commit of the original PDa code for the GSoC Pure Data plugin project of Wincent Balin. Stripped some non-sourcefiles and added a rockbox readme that needs a bit more info from Wincent. Is added to CATEGORIES and viewers, but not yet to SUBDIRS (ie doesn't build yet)
git-svn-id: svn://svn.rockbox.org/rockbox/trunk@21044 a1c6a512-1295-4272-9138-f99709370657
Diffstat (limited to 'apps/plugins/pdbox/PDa/src/d_fftroutine.c')
-rw-r--r--apps/plugins/pdbox/PDa/src/d_fftroutine.c2002
1 files changed, 2002 insertions, 0 deletions
diff --git a/apps/plugins/pdbox/PDa/src/d_fftroutine.c b/apps/plugins/pdbox/PDa/src/d_fftroutine.c
new file mode 100644
index 0000000000..dde24fa358
--- /dev/null
+++ b/apps/plugins/pdbox/PDa/src/d_fftroutine.c
@@ -0,0 +1,2002 @@
1/*****************************************************************************/
2/* */
3/* Fast Fourier Transform */
4/* Network Abstraction, Definitions */
5/* Kevin Peterson, MIT Media Lab, EMS */
6/* UROP - Fall '86 */
7/* REV: 6/12/87(KHP) - To incorporate link list of different sized networks */
8/* */
9/*****************************************************************************/
10
11/*****************************************************************************/
12/* added debug option 5/91 brown@nadia */
13/* change sign at AAA */
14/* */
15/* Fast Fourier Transform */
16/* FFT Network Interaction and Support Modules */
17/* Kevin Peterson, MIT Media Lab, EMS */
18/* UROP - Fall '86 */
19/* REV: 6/12/87(KHP) - Generalized to one procedure call with typed I/O */
20/* */
21/*****************************************************************************/
22
23/* Overview:
24
25 My realization of the FFT involves a representation of a network of
26 "butterfly" elements that takes a set of 'N' sound samples as input and
27 computes the discrete Fourier transform. This network consists of a
28 series of stages (log2 N), each stage consisting of N/2 parallel butterfly
29 elements. Consecutive stages are connected by specific, predetermined flow
30 paths, (see Oppenheim, Schafer for details) and each butterfly element has
31 an associated multiplicative coefficient.
32
33 FFT NETWORK:
34 -----------
35 ____ _ ____ _ ____ _ ____ _ ____
36 o--| |o-| |-o| |o-| |-o| |o-| |-o| |o-| |-o| |--o
37 |reg1| | | |W^r1| | | |reg1| | | |W^r1| | | |reg1|
38 | | | | | | | | | | | | | | | | | | .....
39 | | | | | | | | | | | | | | | | | |
40 o--|____|o-| |-o|____|o-| |-o|____|o-| |-o|____|o-| |-o|____|--o
41 | | | | | | | |
42 | | | | | | | |
43 ____ | | ____ | | ____ | | ____ | | ____
44 o--| |o-| |-o| |o-| |-o| |o-| |-o| |o-| |-o| |--o
45 |reg2| | | |W^r2| | | |reg2| | | |W^r2| | | |reg2|
46 | | | | | | | | | | | | | | | | | | .....
47 | | | | | | | | | | | | | | | | | |
48 o--|____|o-| |-o|____|o-| |-o|____|o-| |-o|____|o-| |-o|____|--o
49 | | | | | | | |
50 | | | | | | | |
51 : : : : : : : : :
52 : : : : : : : : :
53 : : : : : : : : :
54 : : : : : : : : :
55 : : : : : : : : :
56
57 ____ | | ____ | | ____ | | ____ | | ____
58 o--| |o-| |-o| |o-| |-o| |o-| |-o| |o-| |-o| |--o
59 |reg | | | |W^r | | | |reg | | | |W^r | | | |reg |
60 | N/2| | | | N/2| | | | N/2| | | | N/2| | | | N/2| .....
61 | | | | | | | | | | | | | | | | | |
62 o--|____|o-|_|-o|____|o-|_|-o|____|o-|_|-o|____|o-|_|-o|____|--o
63
64 ^ ^ ^ ^
65 Initial | Bttrfly | Rd/Wrt | Bttrfly | Rd/Wrt
66 Buffer | | Register | | Register
67 |____________|____________|____________|
68 |
69 |
70 Interconnect
71 Paths
72
73 The use of "in-place" computation permits one to use only one set of
74 registers realized by an array of complex number structures. To describe
75 the coefficients for each butterfly I am using a two dimensional array
76 (stage, butterfly) of complex numbers. The predetermined stage connections
77 will be described in a two dimensional array of indicies. These indicies
78 will be used to determine the order of reading at each stage of the
79 computation.
80*/
81
82
83/*****************************************************************************/
84/* INCLUDE FILES */
85/*****************************************************************************/
86
87#include <stdio.h>
88#include <math.h>
89#include <stdlib.h>
90
91 /* the following is needed only to declare pd_fft() as exportable in MSW */
92#include "m_pd.h"
93
94/* some basic definitions */
95#ifndef BOOL
96#define BOOL int
97#define TRUE 1
98#define FALSE 0
99#endif
100
101#define SAMPLE float /* data type used in calculation */
102
103#define SHORT_SIZE sizeof(short)
104#define INT_SIZE sizeof(int)
105#define FLOAT_SIZE sizeof(float)
106#define SAMPLE_SIZE sizeof(SAMPLE)
107#define PNTR_SIZE sizeof(char *)
108
109#define PI 3.1415927
110#define TWO_PI 6.2831854
111
112/* type definitions for I/O buffers */
113#define REAL 0 /* real only */
114#define IMAG 2 /* imaginary only */
115#define RECT 8 /* real and imaginary */
116#define MAG 16 /* magnitude only */
117#define PHASE 32 /* phase only */
118#define POLAR 64 /* magnitude and phase*/
119
120/* scale definitions for I/O buffers */
121#define LINEAR 0
122#define DB 1 /* 20log10 */
123
124/* transform direction definition */
125#define FORWARD 1 /* Forward FFT */
126#define INVERSE 2 /* Inverse FFT */
127
128/* window type definitions */
129#define HANNING 1
130#define RECTANGULAR 0
131
132
133
134/* network structure definition */
135
136typedef struct Tfft_net {
137 int n;
138 int stages;
139 int bps;
140 int direction;
141 int window_type;
142 int *load_index;
143 SAMPLE *window, *inv_window;
144 SAMPLE *regr;
145 SAMPLE *regi;
146 SAMPLE **indexpr;
147 SAMPLE **indexpi;
148 SAMPLE **indexqr;
149 SAMPLE **indexqi;
150 SAMPLE *coeffr, *inv_coeffr;
151 SAMPLE *coeffi, *inv_coeffi;
152 struct Tfft_net *next;
153} FFT_NET;
154
155
156void cfft(int trnsfrm_dir, int npnt, int window,
157 float *source_buf, int source_form, int source_scale,
158 float *result_buf, int result_form, int result_scale, int debug);
159
160
161/*****************************************************************************/
162/* GLOBAL DECLARATIONS */
163/*****************************************************************************/
164
165static FFT_NET *firstnet;
166
167/* prototypes */
168
169void net_alloc(FFT_NET *fft_net);
170void net_dealloc(FFT_NET *fft_net);
171int power_of_two(int n);
172void create_hanning(SAMPLE *window, int n, SAMPLE scale);
173void create_rectangular(SAMPLE *window, int n, SAMPLE scale);
174void short_to_float(short *short_buf, float *float_buf, int n);
175void load_registers(FFT_NET *fft_net, float *buf, int buf_form,
176 int buf_scale, int trnsfrm_dir);
177void compute_fft(FFT_NET *fft_net);
178void store_registers(FFT_NET *fft_net, float *buf, int buf_form,
179 int buf_scale, int debug);
180void build_fft_network(FFT_NET *fft_net, int n, int window_type);
181
182/*****************************************************************************/
183/* GENERALIZED FAST FOURIER TRANSFORM MODULE */
184/*****************************************************************************/
185
186void cfft(int trnsfrm_dir, int npnt, int window,
187 float *source_buf, int source_form, int source_scale,
188 float *result_buf, int result_form, int result_scale, int debug)
189
190/* modifies: result_buf
191 effects: Computes npnt FFT specified by form, scale, and dir parameters.
192 Source samples (single precision float) are taken from soure_buf and
193 the transfrmd representation is stored in result_buf (single precision
194 float). The parameters are defined as follows:
195
196 trnsfrm_dir = FORWARD | INVERSE
197 npnt = 2^k for some any positive integer k
198 window = HANNING | RECTANGULAR
199 (RECT = real and imag parts, POLAR = magnitude and phase)
200 source_form = REAL | IMAG | RECT | POLAR
201 result_form = REAL | IMAG | RECT | MAG | PHASE | POLAR
202 xxxxxx_scale= LINEAR | DB ( 20log10 |mag| )
203
204 The input/output buffers are stored in a form appropriate to the type.
205 For example: REAL => {real, real, real ...},
206 MAG => {mag, mag, mag, ... },
207 RECT => {real, imag, real, imag, ... },
208 POLAR => {mag, phase, mag, phase, ... }.
209
210 To look at the magnitude (in db) of a 1024 point FFT of a real time
211 signal we have:
212
213 fft(FORWARD, 1024, RECTANGULAR, input, REAL, LINEAR, output, MAG, DB)
214
215 All possible input and output combinations are possible given the
216 choice of type and scale parameters.
217*/
218
219{
220 FFT_NET *thisnet = (FFT_NET *)0;
221 FFT_NET *lastnet = (FFT_NET *)0;
222
223 /* A linked list of fft networks of different sizes is maintained to
224 avoid building with every call. The network is built on the first
225 call but reused for subsequent calls requesting the same size
226 transformation.
227 */
228
229 thisnet=firstnet;
230 while (thisnet) {
231 if (!(thisnet->n == npnt) || !(thisnet->window_type == window)) {
232 /* current net doesn't match size or window type */
233 lastnet=thisnet;
234 thisnet=thisnet->next;
235 continue; /* keep looking */
236 }
237
238 else { /* network matches desired size */
239 load_registers(thisnet, source_buf, source_form, source_scale,
240 trnsfrm_dir);
241 compute_fft(thisnet); /* do transformation */
242 store_registers(thisnet, result_buf, result_form, result_scale,debug);
243 return;
244 }
245 }
246
247 /* none of existing networks match required size*/
248
249 if (lastnet) { /* add new network to end of list */
250 thisnet = (FFT_NET *)malloc(sizeof(FFT_NET)); /* allocate */
251 thisnet->next = 0;
252 lastnet->next = thisnet; /* add to end of list */
253 }
254 else { /* first network to be created */
255 thisnet=firstnet=(FFT_NET *)malloc(sizeof(FFT_NET)); /* alloc. */
256 thisnet->next = 0;
257 }
258
259 /* build new network and compute transformation */
260 build_fft_network(thisnet, npnt, window);
261 load_registers(thisnet, source_buf, source_form, source_scale,
262 trnsfrm_dir);
263 compute_fft(thisnet);
264 store_registers(thisnet, result_buf, result_form, result_scale,debug);
265 return;
266}
267
268void fft_clear(void)
269
270/* effects: Deallocates all preserved FFT networks. Should be used when
271 finished with all computations.
272*/
273
274{
275 FFT_NET *thisnet, *nextnet;
276
277 if (firstnet) {
278 thisnet=firstnet;
279 do {
280 nextnet = thisnet->next;
281 net_dealloc(thisnet);
282 free((char *)thisnet);
283 } while (thisnet = nextnet);
284 }
285}
286
287
288/*****************************************************************************/
289/* NETWORK CONSTRUCTION */
290/*****************************************************************************/
291
292void build_fft_network(FFT_NET *fft_net, int n, int window_type)
293
294
295/* modifies:fft_net
296 effects: Constructs the fft network as described in fft.h. Butterfly
297 coefficients, read/write indicies, bit reversed load indicies,
298 and array allocations are computed.
299*/
300
301{
302 int cntr, i, j, s;
303 int stages, bps;
304 int **p, **q, *pp, *qp;
305 SAMPLE two_pi_div_n = TWO_PI / n;
306
307
308 /* network definition */
309 fft_net->n = n;
310 fft_net->bps = bps = n/2;
311 for (i = 0, j = n; j > 1; j >>= 1, i++);
312 fft_net->stages = stages = i;
313 fft_net->direction = FORWARD;
314 fft_net->window_type = window_type;
315 fft_net->next = (FFT_NET *)0;
316
317 /* allocate registers, index, coefficient arrays */
318 net_alloc(fft_net);
319
320
321 /* create appropriate windows */
322 if (window_type==HANNING) {
323 create_hanning(fft_net->window, n, 1.);
324 create_hanning(fft_net->inv_window, n, 1./n);
325 }
326 else {
327 create_rectangular(fft_net->window, n, 1.);
328 create_rectangular(fft_net->inv_window, n, 1./n);
329 }
330
331
332 /* calculate butterfly coefficients */ {
333
334 int num_diff_coeffs, power_inc, power;
335 SAMPLE *coeffpr = fft_net->coeffr;
336 SAMPLE *coeffpi = fft_net->coeffi;
337 SAMPLE *inv_coeffpr = fft_net->inv_coeffr;
338 SAMPLE *inv_coeffpi = fft_net->inv_coeffi;
339
340 /* stage one coeffs are 1 + 0j */
341 for (i = 0; i < bps; i++) {
342 *coeffpr = *inv_coeffpr = 1.;
343 *coeffpi = *inv_coeffpi = 0.;
344 coeffpr++; inv_coeffpr++;
345 coeffpi++; inv_coeffpi++;
346 }
347
348 /* stage 2 to last stage coeffs need calculation */
349 /* (1<<r <=> 2^r */
350 for (s = 2; s <= stages; s++) {
351
352 num_diff_coeffs = n / (1 << (stages - s + 1));
353 power_inc = 1 << (stages -s);
354 cntr = 0;
355
356 for (i = bps/num_diff_coeffs; i > 0; i--) {
357
358 power = 0;
359
360 for (j = num_diff_coeffs; j > 0; j--) {
361 *coeffpr = cos(two_pi_div_n*power);
362 *inv_coeffpr = cos(two_pi_div_n*power);
363/* AAA change these signs */ *coeffpi = -sin(two_pi_div_n*power);
364/* change back */ *inv_coeffpi = sin(two_pi_div_n*power);
365 power += power_inc;
366 coeffpr++; inv_coeffpr++;
367 coeffpi++; inv_coeffpi++;
368 }
369 }
370 }
371 }
372
373 /* calculate network indicies: stage exchange indicies are
374 calculated and then used as offset values from the base
375 register locations. The final addresses are then stored in
376 fft_net.
377 */ {
378
379 int index, inc;
380 SAMPLE **indexpr = fft_net->indexpr;
381 SAMPLE **indexpi = fft_net->indexpi;
382 SAMPLE **indexqr = fft_net->indexqr;
383 SAMPLE **indexqi = fft_net->indexqi;
384 SAMPLE *regr = fft_net->regr;
385 SAMPLE *regi = fft_net->regi;
386
387
388 /* allocate temporary 2d stage exchange index, 1d temp
389 load index */
390 p = (int **)malloc(stages * PNTR_SIZE);
391 q = (int **)malloc(stages * PNTR_SIZE);
392
393 for (s = 0; s < stages; s++) {
394 p[s] = (int *)malloc(bps * INT_SIZE);
395 q[s] = (int *)malloc(bps * INT_SIZE);
396 }
397
398 /* calculate stage exchange indicies: */
399 for (s = 0; s < stages; s++) {
400 pp = p[s];
401 qp = q[s];
402 inc = 1 << s;
403 cntr = 1 << (stages-s-1);
404 i = j = index = 0;
405
406 do {
407 do {
408 qp[i] = index + inc;
409 pp[i++] = index++;
410 } while (++j < inc);
411 index = qp[i-1] + 1;
412 j = 0;
413 } while (--cntr);
414 }
415
416 /* compute actual address values using indicies as offsets */
417 for (s = 0; s < stages; s++) {
418 for (i = 0; i < bps; i++) {
419 *indexpr++ = regr + p[s][i];
420 *indexpi++ = regi + p[s][i];
421 *indexqr++ = regr + q[s][i];
422 *indexqi++ = regi + q[s][i];
423 }
424 }
425 }
426
427
428 /* calculate load indicies (bit reverse ordering) */
429 /* bit reverse ordering achieved by passing normal
430 order indicies backwards through the network */
431
432 /* init to normal order indicies */ {
433 int *load_index,*load_indexp;
434 int *temp_indexp, *temp_index;
435 temp_index=temp_indexp=(int *)malloc(n * INT_SIZE);
436
437 i = 0; j = n;
438 load_index = load_indexp = fft_net->load_index;
439
440 while (j--)
441 *load_indexp++ = i++;
442
443 /* pass indicies backwards through net */
444 for (s = stages - 1; s > 0; s--) {
445 pp = p[s];
446 qp = q[s];
447
448 for (i = 0; i < bps; i++) {
449 temp_index[pp[i]]=load_index[2*i];
450 temp_index[qp[i]]=load_index[2*i+1];
451 }
452 j = n;
453 load_indexp = load_index;
454 temp_indexp = temp_index;
455 while (j--)
456 *load_indexp++ = *temp_indexp++;
457 }
458
459 /* free all temporary arrays */
460 free((char *)temp_index);
461 for (s = 0; s < stages; s++) {
462 free((char *)p[s]);free((char *)q[s]);
463 }
464 free((char *)p);free((char *)q);
465 }
466}
467
468
469
470/*****************************************************************************/
471/* REGISTER LOAD AND STORE */
472/*****************************************************************************/
473
474void load_registers(FFT_NET *fft_net, float *buf, int buf_form,
475 int buf_scale, int trnsfrm_dir)
476
477/* effects: Multiplies the input buffer with the appropriate window and
478 stores the resulting values in the initial registers of the
479 network. Input buffer must contain values appropriate to form.
480 For RECT, the buffer contains real num. followed by imag num,
481 and for POLAR, it contains magnitude followed by phase. Pure
482 inputs are listed normally. Both LINEAR and DB scales are
483 interpreted.
484*/
485
486{
487 int *load_index = fft_net->load_index;
488 SAMPLE *window;
489 int index, i = 0, n = fft_net->n;
490
491 if (trnsfrm_dir==FORWARD) window = fft_net->window;
492 else if (trnsfrm_dir==INVERSE) window = fft_net->inv_window;
493 else {
494 fprintf(stderr, "load_registers:illegal transform direction\n");
495 exit(0);
496 }
497 fft_net->direction = trnsfrm_dir;
498
499 switch(buf_scale) {
500 case LINEAR: {
501
502 switch (buf_form) {
503 case REAL: { /* pure REAL */
504 while (i < fft_net->n) {
505 index = load_index[i];
506 fft_net->regr[i]=(SAMPLE)buf[index] * window[index];
507 fft_net->regi[i]=0.;
508 i++;
509 }
510 } break;
511
512 case IMAG: { /* pure IMAGinary */
513 while (i < fft_net->n) {
514 index = load_index[i];
515 fft_net->regr[i]=0;
516 fft_net->regi[i]=(SAMPLE)buf[index] * window[index];
517 i++;
518 }
519 } break;
520
521 case RECT: { /* both REAL and IMAGinary */
522 while (i < fft_net->n) {
523 index = load_index[i];
524 fft_net->regr[i]=(SAMPLE)buf[index*2] * window[index];
525 fft_net->regi[i]=(SAMPLE)buf[index*2+1] * window[index];
526 i++;
527 }
528 } break;
529
530 case POLAR: { /* magnitude followed by phase */
531 while (i < fft_net->n) {
532 index = load_index[i];
533 fft_net->regr[i]=(SAMPLE)(buf[index*2] * cos(buf[index*2+1]))
534 * window[index];
535 fft_net->regi[i]=(SAMPLE)(buf[index*2] * sin(buf[index*2+1]))
536 * window[index];
537 i++;
538 }
539 } break;
540
541 default: {
542 fprintf(stderr, "load_registers:illegal input form\n");
543 exit(0);
544 } break;
545 }
546 } break;
547
548 case DB: {
549
550 switch (buf_form) {
551 case REAL: { /* log pure REAL */
552 while (i < fft_net->n) {
553 index = load_index[i];
554 fft_net->regr[i]=(SAMPLE)pow(10., (1./20.)*buf[index])
555 * window[index]; /* window scaling after linearization */
556 fft_net->regi[i]=0.;
557 i++;
558 }
559 } break;
560
561 case IMAG: { /* log pure IMAGinary */
562 while (i < fft_net->n) {
563 index = load_index[i];
564 fft_net->regr[i]=0.;
565 fft_net->regi[i]=(SAMPLE)pow(10., (1./20.)*buf[index])
566 * window[index];
567 i++;
568 }
569 } break;
570
571 case RECT: { /* log REAL and log IMAGinary */
572 while (i < fft_net->n) {
573 index = load_index[i];
574 fft_net->regr[i]=(SAMPLE)pow(10., (1./20.)*buf[index*2])
575 * window[index];
576 fft_net->regi[i]=(SAMPLE)pow(10., (1./20.)*buf[index*2+1])
577 * window[index];
578 i++;
579 }
580 } break;
581
582 case POLAR: { /* log mag followed by phase */
583 while (i < fft_net->n) {
584 index = load_index[i];
585 fft_net->regr[i]=(SAMPLE)(pow(10., (1./20.)*buf[index*2])
586 * cos(buf[index*2+1])) * window[index];
587 fft_net->regi[i]=(SAMPLE)(pow(10., (1./20.)*buf[index*2])
588 * sin(buf[index*2+1])) * window[index];
589 i++;
590 }
591 } break;
592
593 default: {
594 fprintf(stderr, "load_registers:illegal input form\n");
595 exit(0);
596 } break;
597 }
598 } break;
599
600 default: {
601 fprintf(stderr, "load_registers:illegal input scale\n");
602 exit(0);
603 } break;
604 }
605}
606
607
608void store_registers(FFT_NET *fft_net, float *buf, int buf_form,
609 int buf_scale, int debug)
610
611/* modifies: buf
612 effects: Writes the final contents of the network registers into buf in
613 either linear or db scale, polar or rectangular form. If any of
614 the pure forms(REAL, IMAG, MAG, or PHASE) are used then only the
615 corresponding part of the registers is stored in buf.
616*/
617
618{
619 int i;
620 SAMPLE real, imag, mag, phase;
621 int n;
622
623 i = 0;
624 n = fft_net->n;
625
626 switch (buf_scale) {
627 case LINEAR: {
628
629 switch (buf_form) {
630 case REAL: { /* pure REAL */
631 do {
632 *buf++ = (float)fft_net->regr[i];
633 } while (++i < n);
634 } break;
635
636 case IMAG: { /* pure IMAGinary */
637 do {
638 *buf++ = (float)fft_net->regi[i];
639 } while (++i < n);
640 } break;
641
642 case RECT: { /* both REAL and IMAGinary */
643 do {
644 *buf++ = (float)fft_net->regr[i];
645 *buf++ = (float)fft_net->regi[i];
646 } while (++i < n);
647 } break;
648
649 case MAG: { /* magnitude only */
650 do {
651 real = fft_net->regr[i];
652 imag = fft_net->regi[i];
653 *buf++ = (float)sqrt(real*real+imag*imag);
654 } while (++i < n);
655 } break;
656
657 case PHASE: { /* phase only */
658 do {
659 real = fft_net->regr[i];
660 imag = fft_net->regi[i];
661 if (real > .00001)
662 *buf++ = (float)atan2(imag, real);
663 else { /* deal with bad case */
664 if (imag > 0){ *buf++ = PI / 2.;
665 if(debug) fprintf(stderr,"real=0 and imag > 0\n");}
666 else if (imag < 0){ *buf++ = -PI / 2.;
667 if(debug) fprintf(stderr,"real=0 and imag < 0\n");}
668 else { *buf++ = 0;
669 if(debug) fprintf(stderr,"real=0 and imag=0\n");}
670 }
671 } while (++i < n);
672 } break;
673
674 case POLAR: { /* magnitude and phase */
675 do {
676 real = fft_net->regr[i];
677 imag = fft_net->regi[i];
678 *buf++ = (float)sqrt(real*real+imag*imag);
679 if (real) /* a hack to avoid div by zero */
680 *buf++ = (float)atan2(imag, real);
681 else { /* deal with bad case */
682 if (imag > 0) *buf++ = PI / 2.;
683 else if (imag < 0) *buf++ = -PI / 2.;
684 else *buf++ = 0;
685 }
686 } while (++i < n);
687 } break;
688
689 default: {
690 fprintf(stderr, "store_registers:illegal output form\n");
691 exit(0);
692 } break;
693 }
694 } break;
695
696 case DB: {
697
698 switch (buf_form) {
699 case REAL: { /* real only */
700 do {
701 *buf++ = (float)20.*log10(fft_net->regr[i]);
702 } while (++i < n);
703 } break;
704
705 case IMAG: { /* imag only */
706 do {
707 *buf++ = (float)20.*log10(fft_net->regi[i]);
708 } while (++i < n);
709 } break;
710
711 case RECT: { /* real and imag */
712 do {
713 *buf++ = (float)20.*log10(fft_net->regr[i]);
714 *buf++ = (float)20.*log10(fft_net->regi[i]);
715 } while (++i < n);
716 } break;
717
718 case MAG: { /* magnitude only */
719 do {
720 real = fft_net->regr[i];
721 imag = fft_net->regi[i];
722 *buf++ = (float)20.*log10(sqrt(real*real+imag*imag));
723 } while (++i < n);
724 } break;
725
726 case PHASE: { /* phase only */
727 do {
728 real = fft_net->regr[i];
729 imag = fft_net->regi[i];
730 if (real)
731 *buf++ = (float)atan2(imag, real);
732 else { /* deal with bad case */
733 if (imag > 0) *buf++ = PI / 2.;
734 else if (imag < 0) *buf++ = -PI / 2.;
735 else *buf++ = 0;
736 }
737 } while (++i < n);
738 } break;
739
740 case POLAR: { /* magnitude and phase */
741 do {
742 real = fft_net->regr[i];
743 imag = fft_net->regi[i];
744 *buf++ = (float)20.*log10(sqrt(real*real+imag*imag));
745 if (real)
746 *buf++ = (float)atan2(imag, real);
747 else { /* deal with bad case */
748 if (imag > 0) *buf++ = PI / 2.;
749 else if (imag < 0) *buf++ = -PI / 2.;
750 else *buf++ = 0;
751 }
752 } while (++i < n);
753 } break;
754
755 default: {
756 fprintf(stderr, "store_registers:illegal output form\n");
757 exit(0);
758 } break;
759 }
760 } break;
761
762 default: {
763 fprintf(stderr, "store_registers:illegal output scale\n");
764 exit(0);
765 } break;
766 }
767}
768
769
770
771/*****************************************************************************/
772/* COMPUTE TRANSFORMATION */
773/*****************************************************************************/
774
775void compute_fft(FFT_NET *fft_net)
776
777
778/* modifies: fft_net
779 effects: Passes the values (already loaded) in the registers through
780 the network, multiplying with appropriate coefficients at each
781 stage. The fft result will be in the registers at the end of
782 the computation. The direction of the transformation is indicated
783 by the network flag 'direction'. The form of the computation is:
784
785 X(pn) = X(p) + C*X(q)
786 X(qn) = X(p) - C*X(q)
787
788 where X(pn,qn) represents the output of the registers at each stage.
789 The calculations are actually done in place. Register pointers are
790 used to speed up the calculations.
791
792 Register and coefficient addresses involved in the calculations
793 are stored sequentially and are accessed as such. fft_net->indexp,
794 indexq contain pointers to the relevant addresses, and fft_net->coeffs,
795 inv_coeffs points to the appropriate coefficients at each stage of the
796 computation.
797*/
798
799{
800 SAMPLE **xpr, **xpi, **xqr, **xqi, *cr, *ci;
801 int i;
802 SAMPLE tpr, tpi, tqr, tqi;
803 int bps = fft_net->bps;
804 int cnt = bps * (fft_net->stages - 1);
805
806 /* predetermined register addresses and coefficients */
807 xpr = fft_net->indexpr;
808 xpi = fft_net->indexpi;
809 xqr = fft_net->indexqr;
810 xqi = fft_net->indexqi;
811
812 if (fft_net->direction==FORWARD) { /* FORWARD FFT coefficients */
813 cr = fft_net->coeffr;
814 ci = fft_net->coeffi;
815 }
816 else { /* INVERSE FFT coefficients */
817 cr = fft_net->inv_coeffr;
818 ci = fft_net->inv_coeffi;
819 }
820
821 /* stage one coefficients are 1 + 0j so C*X(q)=X(q) */
822 /* bps mults can be avoided */
823
824 for (i = 0; i < bps; i++) {
825
826 /* add X(p) and X(q) */
827 tpr = **xpr + **xqr;
828 tpi = **xpi + **xqi;
829 tqr = **xpr - **xqr;
830 tqi = **xpi - **xqi;
831
832 /* exchange register with temp */
833 **xpr = tpr;
834 **xpi = tpi;
835 **xqr = tqr;
836 **xqi = tqi;
837
838 /* next set of register for calculations: */
839 xpr++; xpi++; xqr++; xqi++; cr++; ci++;
840
841 }
842
843 for (i = 0; i < cnt; i++) {
844
845 /* mult X(q) by coeff C */
846 tqr = **xqr * *cr - **xqi * *ci;
847 tqi = **xqr * *ci + **xqi * *cr;
848
849 /* exchange register with temp */
850 **xqr = tqr;
851 **xqi = tqi;
852
853 /* add X(p) and X(q) */
854 tpr = **xpr + **xqr;
855 tpi = **xpi + **xqi;
856 tqr = **xpr - **xqr;
857 tqi = **xpi - **xqi;
858
859 /* exchange register with temp */
860 **xpr = tpr;
861 **xpi = tpi;
862 **xqr = tqr;
863 **xqi = tqi;
864 /* next set of register for calculations: */
865 xpr++; xpi++; xqr++; xqi++; cr++; ci++;
866 }
867}
868
869
870/****************************************************************************/
871/* SUPPORT MODULES */
872/****************************************************************************/
873
874void net_alloc(FFT_NET *fft_net)
875
876
877/* effects: Allocates appropriate two dimensional arrays and assigns
878 correct internal pointers.
879*/
880
881{
882
883 int stages, bps, n;
884
885 n = fft_net->n;
886 stages = fft_net->stages;
887 bps = fft_net->bps;
888
889
890 /* two dimensional arrays with elements stored sequentially */
891
892 fft_net->load_index = (int *)malloc(n * INT_SIZE);
893 fft_net->regr = (SAMPLE *)malloc(n * SAMPLE_SIZE);
894 fft_net->regi = (SAMPLE *)malloc(n * SAMPLE_SIZE);
895 fft_net->coeffr = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
896 fft_net->coeffi = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
897 fft_net->inv_coeffr = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
898 fft_net->inv_coeffi = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
899 fft_net->indexpr = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
900 fft_net->indexpi = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
901 fft_net->indexqr = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
902 fft_net->indexqi = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
903
904 /* one dimensional load window */
905 fft_net->window = (SAMPLE *)malloc(n * SAMPLE_SIZE);
906 fft_net->inv_window = (SAMPLE *)malloc(n * SAMPLE_SIZE);
907}
908
909void net_dealloc(FFT_NET *fft_net)
910
911
912/* effects: Deallocates given FFT network.
913*/
914
915{
916
917 free((char *)fft_net->load_index);
918 free((char *)fft_net->regr);
919 free((char *)fft_net->regi);
920 free((char *)fft_net->coeffr);
921 free((char *)fft_net->coeffi);
922 free((char *)fft_net->inv_coeffr);
923 free((char *)fft_net->inv_coeffi);
924 free((char *)fft_net->indexpr);
925 free((char *)fft_net->indexpi);
926 free((char *)fft_net->indexqr);
927 free((char *)fft_net->indexqi);
928 free((char *)fft_net->window);
929 free((char *)fft_net->inv_window);
930}
931
932
933BOOL power_of_two(n)
934
935int n;
936
937/* effects: Returns TRUE if n is a power of two, otherwise FALSE.
938*/
939
940{
941 int i;
942
943 for (i = n; i > 1; i >>= 1)
944 if (i & 1) return FALSE; /* more than one bit high */
945 return TRUE;
946}
947
948
949void create_hanning(SAMPLE *window, int n, SAMPLE scale)
950
951/* effects: Fills the buffer window with a hanning window of the appropriate
952 size scaled by scale.
953*/
954
955{
956 SAMPLE a, pi_div_n = PI/n;
957 int k;
958
959 for (k=1; k <= n; k++) {
960 a = sin(k * pi_div_n);
961 *window++ = scale * a * a;
962 }
963}
964
965
966void create_rectangular(SAMPLE *window, int n, SAMPLE scale)
967
968/* effects: Fills the buffer window with a rectangular window of the
969 appropriate size of height scale.
970*/
971
972{
973 while (n--)
974 *window++ = scale;
975}
976
977
978void short_to_float(short *short_buf, float *float_buf, int n)
979
980/* effects; Converts short_buf to floats and stores them in float_buf.
981*/
982
983{
984 while (n--) {
985 *float_buf++ = (float)*short_buf++;
986 }
987}
988
989
990/* here's the meat: */
991
992void pd_fft(float *buf, int npoints, int inverse)
993{
994 double renorm;
995 float *fp, *fp2;
996 int i;
997 renorm = (inverse ? npoints : 1.);
998 cfft((inverse ? INVERSE : FORWARD), npoints, RECTANGULAR,
999 buf, RECT, LINEAR, buf, RECT, LINEAR, 0);
1000 for (i = npoints << 1, fp = buf; i--; fp++) *fp *= renorm;
1001}
1002/*****************************************************************************/
1003/* */
1004/* Fast Fourier Transform */
1005/* Network Abstraction, Definitions */
1006/* Kevin Peterson, MIT Media Lab, EMS */
1007/* UROP - Fall '86 */
1008/* REV: 6/12/87(KHP) - To incorporate link list of different sized networks */
1009/* */
1010/*****************************************************************************/
1011
1012/*****************************************************************************/
1013/* added debug option 5/91 brown@nadia */
1014/* change sign at AAA */
1015/* */
1016/* Fast Fourier Transform */
1017/* FFT Network Interaction and Support Modules */
1018/* Kevin Peterson, MIT Media Lab, EMS */
1019/* UROP - Fall '86 */
1020/* REV: 6/12/87(KHP) - Generalized to one procedure call with typed I/O */
1021/* */
1022/*****************************************************************************/
1023
1024/* Overview:
1025
1026 My realization of the FFT involves a representation of a network of
1027 "butterfly" elements that takes a set of 'N' sound samples as input and
1028 computes the discrete Fourier transform. This network consists of a
1029 series of stages (log2 N), each stage consisting of N/2 parallel butterfly
1030 elements. Consecutive stages are connected by specific, predetermined flow
1031 paths, (see Oppenheim, Schafer for details) and each butterfly element has
1032 an associated multiplicative coefficient.
1033
1034 FFT NETWORK:
1035 -----------
1036 ____ _ ____ _ ____ _ ____ _ ____
1037 o--| |o-| |-o| |o-| |-o| |o-| |-o| |o-| |-o| |--o
1038 |reg1| | | |W^r1| | | |reg1| | | |W^r1| | | |reg1|
1039 | | | | | | | | | | | | | | | | | | .....
1040 | | | | | | | | | | | | | | | | | |
1041 o--|____|o-| |-o|____|o-| |-o|____|o-| |-o|____|o-| |-o|____|--o
1042 | | | | | | | |
1043 | | | | | | | |
1044 ____ | | ____ | | ____ | | ____ | | ____
1045 o--| |o-| |-o| |o-| |-o| |o-| |-o| |o-| |-o| |--o
1046 |reg2| | | |W^r2| | | |reg2| | | |W^r2| | | |reg2|
1047 | | | | | | | | | | | | | | | | | | .....
1048 | | | | | | | | | | | | | | | | | |
1049 o--|____|o-| |-o|____|o-| |-o|____|o-| |-o|____|o-| |-o|____|--o
1050 | | | | | | | |
1051 | | | | | | | |
1052 : : : : : : : : :
1053 : : : : : : : : :
1054 : : : : : : : : :
1055 : : : : : : : : :
1056 : : : : : : : : :
1057
1058 ____ | | ____ | | ____ | | ____ | | ____
1059 o--| |o-| |-o| |o-| |-o| |o-| |-o| |o-| |-o| |--o
1060 |reg | | | |W^r | | | |reg | | | |W^r | | | |reg |
1061 | N/2| | | | N/2| | | | N/2| | | | N/2| | | | N/2| .....
1062 | | | | | | | | | | | | | | | | | |
1063 o--|____|o-|_|-o|____|o-|_|-o|____|o-|_|-o|____|o-|_|-o|____|--o
1064
1065 ^ ^ ^ ^
1066 Initial | Bttrfly | Rd/Wrt | Bttrfly | Rd/Wrt
1067 Buffer | | Register | | Register
1068 |____________|____________|____________|
1069 |
1070 |
1071 Interconnect
1072 Paths
1073
1074 The use of "in-place" computation permits one to use only one set of
1075 registers realized by an array of complex number structures. To describe
1076 the coefficients for each butterfly I am using a two dimensional array
1077 (stage, butterfly) of complex numbers. The predetermined stage connections
1078 will be described in a two dimensional array of indicies. These indicies
1079 will be used to determine the order of reading at each stage of the
1080 computation.
1081*/
1082
1083
1084/*****************************************************************************/
1085/* INCLUDE FILES */
1086/*****************************************************************************/
1087
1088#include <stdio.h>
1089#include <math.h>
1090#include <stdlib.h>
1091
1092 /* the following is needed only to declare pd_fft() as exportable in MSW */
1093#include "m_pd.h"
1094
1095/* some basic definitions */
1096#ifndef BOOL
1097#define BOOL int
1098#define TRUE 1
1099#define FALSE 0
1100#endif
1101
1102#define SAMPLE float /* data type used in calculation */
1103
1104#define SHORT_SIZE sizeof(short)
1105#define INT_SIZE sizeof(int)
1106#define FLOAT_SIZE sizeof(float)
1107#define SAMPLE_SIZE sizeof(SAMPLE)
1108#define PNTR_SIZE sizeof(char *)
1109
1110#define PI 3.1415927
1111#define TWO_PI 6.2831854
1112
1113/* type definitions for I/O buffers */
1114#define REAL 0 /* real only */
1115#define IMAG 2 /* imaginary only */
1116#define RECT 8 /* real and imaginary */
1117#define MAG 16 /* magnitude only */
1118#define PHASE 32 /* phase only */
1119#define POLAR 64 /* magnitude and phase*/
1120
1121/* scale definitions for I/O buffers */
1122#define LINEAR 0
1123#define DB 1 /* 20log10 */
1124
1125/* transform direction definition */
1126#define FORWARD 1 /* Forward FFT */
1127#define INVERSE 2 /* Inverse FFT */
1128
1129/* window type definitions */
1130#define HANNING 1
1131#define RECTANGULAR 0
1132
1133
1134
1135/* network structure definition */
1136
1137typedef struct Tfft_net {
1138 int n;
1139 int stages;
1140 int bps;
1141 int direction;
1142 int window_type;
1143 int *load_index;
1144 SAMPLE *window, *inv_window;
1145 SAMPLE *regr;
1146 SAMPLE *regi;
1147 SAMPLE **indexpr;
1148 SAMPLE **indexpi;
1149 SAMPLE **indexqr;
1150 SAMPLE **indexqi;
1151 SAMPLE *coeffr, *inv_coeffr;
1152 SAMPLE *coeffi, *inv_coeffi;
1153 struct Tfft_net *next;
1154} FFT_NET;
1155
1156
1157void cfft(int trnsfrm_dir, int npnt, int window,
1158 float *source_buf, int source_form, int source_scale,
1159 float *result_buf, int result_form, int result_scale, int debug);
1160
1161
1162/*****************************************************************************/
1163/* GLOBAL DECLARATIONS */
1164/*****************************************************************************/
1165
1166static FFT_NET *firstnet;
1167
1168/* prototypes */
1169
1170void net_alloc(FFT_NET *fft_net);
1171void net_dealloc(FFT_NET *fft_net);
1172int power_of_two(int n);
1173void create_hanning(SAMPLE *window, int n, SAMPLE scale);
1174void create_rectangular(SAMPLE *window, int n, SAMPLE scale);
1175void short_to_float(short *short_buf, float *float_buf, int n);
1176void load_registers(FFT_NET *fft_net, float *buf, int buf_form,
1177 int buf_scale, int trnsfrm_dir);
1178void compute_fft(FFT_NET *fft_net);
1179void store_registers(FFT_NET *fft_net, float *buf, int buf_form,
1180 int buf_scale, int debug);
1181void build_fft_network(FFT_NET *fft_net, int n, int window_type);
1182
1183/*****************************************************************************/
1184/* GENERALIZED FAST FOURIER TRANSFORM MODULE */
1185/*****************************************************************************/
1186
1187void cfft(int trnsfrm_dir, int npnt, int window,
1188 float *source_buf, int source_form, int source_scale,
1189 float *result_buf, int result_form, int result_scale, int debug)
1190
1191/* modifies: result_buf
1192 effects: Computes npnt FFT specified by form, scale, and dir parameters.
1193 Source samples (single precision float) are taken from soure_buf and
1194 the transfrmd representation is stored in result_buf (single precision
1195 float). The parameters are defined as follows:
1196
1197 trnsfrm_dir = FORWARD | INVERSE
1198 npnt = 2^k for some any positive integer k
1199 window = HANNING | RECTANGULAR
1200 (RECT = real and imag parts, POLAR = magnitude and phase)
1201 source_form = REAL | IMAG | RECT | POLAR
1202 result_form = REAL | IMAG | RECT | MAG | PHASE | POLAR
1203 xxxxxx_scale= LINEAR | DB ( 20log10 |mag| )
1204
1205 The input/output buffers are stored in a form appropriate to the type.
1206 For example: REAL => {real, real, real ...},
1207 MAG => {mag, mag, mag, ... },
1208 RECT => {real, imag, real, imag, ... },
1209 POLAR => {mag, phase, mag, phase, ... }.
1210
1211 To look at the magnitude (in db) of a 1024 point FFT of a real time
1212 signal we have:
1213
1214 fft(FORWARD, 1024, RECTANGULAR, input, REAL, LINEAR, output, MAG, DB)
1215
1216 All possible input and output combinations are possible given the
1217 choice of type and scale parameters.
1218*/
1219
1220{
1221 FFT_NET *thisnet = (FFT_NET *)0;
1222 FFT_NET *lastnet = (FFT_NET *)0;
1223
1224 /* A linked list of fft networks of different sizes is maintained to
1225 avoid building with every call. The network is built on the first
1226 call but reused for subsequent calls requesting the same size
1227 transformation.
1228 */
1229
1230 thisnet=firstnet;
1231 while (thisnet) {
1232 if (!(thisnet->n == npnt) || !(thisnet->window_type == window)) {
1233 /* current net doesn't match size or window type */
1234 lastnet=thisnet;
1235 thisnet=thisnet->next;
1236 continue; /* keep looking */
1237 }
1238
1239 else { /* network matches desired size */
1240 load_registers(thisnet, source_buf, source_form, source_scale,
1241 trnsfrm_dir);
1242 compute_fft(thisnet); /* do transformation */
1243 store_registers(thisnet, result_buf, result_form, result_scale,debug);
1244 return;
1245 }
1246 }
1247
1248 /* none of existing networks match required size*/
1249
1250 if (lastnet) { /* add new network to end of list */
1251 thisnet = (FFT_NET *)malloc(sizeof(FFT_NET)); /* allocate */
1252 thisnet->next = 0;
1253 lastnet->next = thisnet; /* add to end of list */
1254 }
1255 else { /* first network to be created */
1256 thisnet=firstnet=(FFT_NET *)malloc(sizeof(FFT_NET)); /* alloc. */
1257 thisnet->next = 0;
1258 }
1259
1260 /* build new network and compute transformation */
1261 build_fft_network(thisnet, npnt, window);
1262 load_registers(thisnet, source_buf, source_form, source_scale,
1263 trnsfrm_dir);
1264 compute_fft(thisnet);
1265 store_registers(thisnet, result_buf, result_form, result_scale,debug);
1266 return;
1267}
1268
1269void fft_clear(void)
1270
1271/* effects: Deallocates all preserved FFT networks. Should be used when
1272 finished with all computations.
1273*/
1274
1275{
1276 FFT_NET *thisnet, *nextnet;
1277
1278 if (firstnet) {
1279 thisnet=firstnet;
1280 do {
1281 nextnet = thisnet->next;
1282 net_dealloc(thisnet);
1283 free((char *)thisnet);
1284 } while (thisnet = nextnet);
1285 }
1286}
1287
1288
1289/*****************************************************************************/
1290/* NETWORK CONSTRUCTION */
1291/*****************************************************************************/
1292
1293void build_fft_network(FFT_NET *fft_net, int n, int window_type)
1294
1295
1296/* modifies:fft_net
1297 effects: Constructs the fft network as described in fft.h. Butterfly
1298 coefficients, read/write indicies, bit reversed load indicies,
1299 and array allocations are computed.
1300*/
1301
1302{
1303 int cntr, i, j, s;
1304 int stages, bps;
1305 int **p, **q, *pp, *qp;
1306 SAMPLE two_pi_div_n = TWO_PI / n;
1307
1308
1309 /* network definition */
1310 fft_net->n = n;
1311 fft_net->bps = bps = n/2;
1312 for (i = 0, j = n; j > 1; j >>= 1, i++);
1313 fft_net->stages = stages = i;
1314 fft_net->direction = FORWARD;
1315 fft_net->window_type = window_type;
1316 fft_net->next = (FFT_NET *)0;
1317
1318 /* allocate registers, index, coefficient arrays */
1319 net_alloc(fft_net);
1320
1321
1322 /* create appropriate windows */
1323 if (window_type==HANNING) {
1324 create_hanning(fft_net->window, n, 1.);
1325 create_hanning(fft_net->inv_window, n, 1./n);
1326 }
1327 else {
1328 create_rectangular(fft_net->window, n, 1.);
1329 create_rectangular(fft_net->inv_window, n, 1./n);
1330 }
1331
1332
1333 /* calculate butterfly coefficients */ {
1334
1335 int num_diff_coeffs, power_inc, power;
1336 SAMPLE *coeffpr = fft_net->coeffr;
1337 SAMPLE *coeffpi = fft_net->coeffi;
1338 SAMPLE *inv_coeffpr = fft_net->inv_coeffr;
1339 SAMPLE *inv_coeffpi = fft_net->inv_coeffi;
1340
1341 /* stage one coeffs are 1 + 0j */
1342 for (i = 0; i < bps; i++) {
1343 *coeffpr = *inv_coeffpr = 1.;
1344 *coeffpi = *inv_coeffpi = 0.;
1345 coeffpr++; inv_coeffpr++;
1346 coeffpi++; inv_coeffpi++;
1347 }
1348
1349 /* stage 2 to last stage coeffs need calculation */
1350 /* (1<<r <=> 2^r */
1351 for (s = 2; s <= stages; s++) {
1352
1353 num_diff_coeffs = n / (1 << (stages - s + 1));
1354 power_inc = 1 << (stages -s);
1355 cntr = 0;
1356
1357 for (i = bps/num_diff_coeffs; i > 0; i--) {
1358
1359 power = 0;
1360
1361 for (j = num_diff_coeffs; j > 0; j--) {
1362 *coeffpr = cos(two_pi_div_n*power);
1363 *inv_coeffpr = cos(two_pi_div_n*power);
1364/* AAA change these signs */ *coeffpi = -sin(two_pi_div_n*power);
1365/* change back */ *inv_coeffpi = sin(two_pi_div_n*power);
1366 power += power_inc;
1367 coeffpr++; inv_coeffpr++;
1368 coeffpi++; inv_coeffpi++;
1369 }
1370 }
1371 }
1372 }
1373
1374 /* calculate network indicies: stage exchange indicies are
1375 calculated and then used as offset values from the base
1376 register locations. The final addresses are then stored in
1377 fft_net.
1378 */ {
1379
1380 int index, inc;
1381 SAMPLE **indexpr = fft_net->indexpr;
1382 SAMPLE **indexpi = fft_net->indexpi;
1383 SAMPLE **indexqr = fft_net->indexqr;
1384 SAMPLE **indexqi = fft_net->indexqi;
1385 SAMPLE *regr = fft_net->regr;
1386 SAMPLE *regi = fft_net->regi;
1387
1388
1389 /* allocate temporary 2d stage exchange index, 1d temp
1390 load index */
1391 p = (int **)malloc(stages * PNTR_SIZE);
1392 q = (int **)malloc(stages * PNTR_SIZE);
1393
1394 for (s = 0; s < stages; s++) {
1395 p[s] = (int *)malloc(bps * INT_SIZE);
1396 q[s] = (int *)malloc(bps * INT_SIZE);
1397 }
1398
1399 /* calculate stage exchange indicies: */
1400 for (s = 0; s < stages; s++) {
1401 pp = p[s];
1402 qp = q[s];
1403 inc = 1 << s;
1404 cntr = 1 << (stages-s-1);
1405 i = j = index = 0;
1406
1407 do {
1408 do {
1409 qp[i] = index + inc;
1410 pp[i++] = index++;
1411 } while (++j < inc);
1412 index = qp[i-1] + 1;
1413 j = 0;
1414 } while (--cntr);
1415 }
1416
1417 /* compute actual address values using indicies as offsets */
1418 for (s = 0; s < stages; s++) {
1419 for (i = 0; i < bps; i++) {
1420 *indexpr++ = regr + p[s][i];
1421 *indexpi++ = regi + p[s][i];
1422 *indexqr++ = regr + q[s][i];
1423 *indexqi++ = regi + q[s][i];
1424 }
1425 }
1426 }
1427
1428
1429 /* calculate load indicies (bit reverse ordering) */
1430 /* bit reverse ordering achieved by passing normal
1431 order indicies backwards through the network */
1432
1433 /* init to normal order indicies */ {
1434 int *load_index,*load_indexp;
1435 int *temp_indexp, *temp_index;
1436 temp_index=temp_indexp=(int *)malloc(n * INT_SIZE);
1437
1438 i = 0; j = n;
1439 load_index = load_indexp = fft_net->load_index;
1440
1441 while (j--)
1442 *load_indexp++ = i++;
1443
1444 /* pass indicies backwards through net */
1445 for (s = stages - 1; s > 0; s--) {
1446 pp = p[s];
1447 qp = q[s];
1448
1449 for (i = 0; i < bps; i++) {
1450 temp_index[pp[i]]=load_index[2*i];
1451 temp_index[qp[i]]=load_index[2*i+1];
1452 }
1453 j = n;
1454 load_indexp = load_index;
1455 temp_indexp = temp_index;
1456 while (j--)
1457 *load_indexp++ = *temp_indexp++;
1458 }
1459
1460 /* free all temporary arrays */
1461 free((char *)temp_index);
1462 for (s = 0; s < stages; s++) {
1463 free((char *)p[s]);free((char *)q[s]);
1464 }
1465 free((char *)p);free((char *)q);
1466 }
1467}
1468
1469
1470
1471/*****************************************************************************/
1472/* REGISTER LOAD AND STORE */
1473/*****************************************************************************/
1474
1475void load_registers(FFT_NET *fft_net, float *buf, int buf_form,
1476 int buf_scale, int trnsfrm_dir)
1477
1478/* effects: Multiplies the input buffer with the appropriate window and
1479 stores the resulting values in the initial registers of the
1480 network. Input buffer must contain values appropriate to form.
1481 For RECT, the buffer contains real num. followed by imag num,
1482 and for POLAR, it contains magnitude followed by phase. Pure
1483 inputs are listed normally. Both LINEAR and DB scales are
1484 interpreted.
1485*/
1486
1487{
1488 int *load_index = fft_net->load_index;
1489 SAMPLE *window;
1490 int index, i = 0, n = fft_net->n;
1491
1492 if (trnsfrm_dir==FORWARD) window = fft_net->window;
1493 else if (trnsfrm_dir==INVERSE) window = fft_net->inv_window;
1494 else {
1495 fprintf(stderr, "load_registers:illegal transform direction\n");
1496 exit(0);
1497 }
1498 fft_net->direction = trnsfrm_dir;
1499
1500 switch(buf_scale) {
1501 case LINEAR: {
1502
1503 switch (buf_form) {
1504 case REAL: { /* pure REAL */
1505 while (i < fft_net->n) {
1506 index = load_index[i];
1507 fft_net->regr[i]=(SAMPLE)buf[index] * window[index];
1508 fft_net->regi[i]=0.;
1509 i++;
1510 }
1511 } break;
1512
1513 case IMAG: { /* pure IMAGinary */
1514 while (i < fft_net->n) {
1515 index = load_index[i];
1516 fft_net->regr[i]=0;
1517 fft_net->regi[i]=(SAMPLE)buf[index] * window[index];
1518 i++;
1519 }
1520 } break;
1521
1522 case RECT: { /* both REAL and IMAGinary */
1523 while (i < fft_net->n) {
1524 index = load_index[i];
1525 fft_net->regr[i]=(SAMPLE)buf[index*2] * window[index];
1526 fft_net->regi[i]=(SAMPLE)buf[index*2+1] * window[index];
1527 i++;
1528 }
1529 } break;
1530
1531 case POLAR: { /* magnitude followed by phase */
1532 while (i < fft_net->n) {
1533 index = load_index[i];
1534 fft_net->regr[i]=(SAMPLE)(buf[index*2] * cos(buf[index*2+1]))
1535 * window[index];
1536 fft_net->regi[i]=(SAMPLE)(buf[index*2] * sin(buf[index*2+1]))
1537 * window[index];
1538 i++;
1539 }
1540 } break;
1541
1542 default: {
1543 fprintf(stderr, "load_registers:illegal input form\n");
1544 exit(0);
1545 } break;
1546 }
1547 } break;
1548
1549 case DB: {
1550
1551 switch (buf_form) {
1552 case REAL: { /* log pure REAL */
1553 while (i < fft_net->n) {
1554 index = load_index[i];
1555 fft_net->regr[i]=(SAMPLE)pow(10., (1./20.)*buf[index])
1556 * window[index]; /* window scaling after linearization */
1557 fft_net->regi[i]=0.;
1558 i++;
1559 }
1560 } break;
1561
1562 case IMAG: { /* log pure IMAGinary */
1563 while (i < fft_net->n) {
1564 index = load_index[i];
1565 fft_net->regr[i]=0.;
1566 fft_net->regi[i]=(SAMPLE)pow(10., (1./20.)*buf[index])
1567 * window[index];
1568 i++;
1569 }
1570 } break;
1571
1572 case RECT: { /* log REAL and log IMAGinary */
1573 while (i < fft_net->n) {
1574 index = load_index[i];
1575 fft_net->regr[i]=(SAMPLE)pow(10., (1./20.)*buf[index*2])
1576 * window[index];
1577 fft_net->regi[i]=(SAMPLE)pow(10., (1./20.)*buf[index*2+1])
1578 * window[index];
1579 i++;
1580 }
1581 } break;
1582
1583 case POLAR: { /* log mag followed by phase */
1584 while (i < fft_net->n) {
1585 index = load_index[i];
1586 fft_net->regr[i]=(SAMPLE)(pow(10., (1./20.)*buf[index*2])
1587 * cos(buf[index*2+1])) * window[index];
1588 fft_net->regi[i]=(SAMPLE)(pow(10., (1./20.)*buf[index*2])
1589 * sin(buf[index*2+1])) * window[index];
1590 i++;
1591 }
1592 } break;
1593
1594 default: {
1595 fprintf(stderr, "load_registers:illegal input form\n");
1596 exit(0);
1597 } break;
1598 }
1599 } break;
1600
1601 default: {
1602 fprintf(stderr, "load_registers:illegal input scale\n");
1603 exit(0);
1604 } break;
1605 }
1606}
1607
1608
1609void store_registers(FFT_NET *fft_net, float *buf, int buf_form,
1610 int buf_scale, int debug)
1611
1612/* modifies: buf
1613 effects: Writes the final contents of the network registers into buf in
1614 either linear or db scale, polar or rectangular form. If any of
1615 the pure forms(REAL, IMAG, MAG, or PHASE) are used then only the
1616 corresponding part of the registers is stored in buf.
1617*/
1618
1619{
1620 int i;
1621 SAMPLE real, imag, mag, phase;
1622 int n;
1623
1624 i = 0;
1625 n = fft_net->n;
1626
1627 switch (buf_scale) {
1628 case LINEAR: {
1629
1630 switch (buf_form) {
1631 case REAL: { /* pure REAL */
1632 do {
1633 *buf++ = (float)fft_net->regr[i];
1634 } while (++i < n);
1635 } break;
1636
1637 case IMAG: { /* pure IMAGinary */
1638 do {
1639 *buf++ = (float)fft_net->regi[i];
1640 } while (++i < n);
1641 } break;
1642
1643 case RECT: { /* both REAL and IMAGinary */
1644 do {
1645 *buf++ = (float)fft_net->regr[i];
1646 *buf++ = (float)fft_net->regi[i];
1647 } while (++i < n);
1648 } break;
1649
1650 case MAG: { /* magnitude only */
1651 do {
1652 real = fft_net->regr[i];
1653 imag = fft_net->regi[i];
1654 *buf++ = (float)sqrt(real*real+imag*imag);
1655 } while (++i < n);
1656 } break;
1657
1658 case PHASE: { /* phase only */
1659 do {
1660 real = fft_net->regr[i];
1661 imag = fft_net->regi[i];
1662 if (real > .00001)
1663 *buf++ = (float)atan2(imag, real);
1664 else { /* deal with bad case */
1665 if (imag > 0){ *buf++ = PI / 2.;
1666 if(debug) fprintf(stderr,"real=0 and imag > 0\n");}
1667 else if (imag < 0){ *buf++ = -PI / 2.;
1668 if(debug) fprintf(stderr,"real=0 and imag < 0\n");}
1669 else { *buf++ = 0;
1670 if(debug) fprintf(stderr,"real=0 and imag=0\n");}
1671 }
1672 } while (++i < n);
1673 } break;
1674
1675 case POLAR: { /* magnitude and phase */
1676 do {
1677 real = fft_net->regr[i];
1678 imag = fft_net->regi[i];
1679 *buf++ = (float)sqrt(real*real+imag*imag);
1680 if (real) /* a hack to avoid div by zero */
1681 *buf++ = (float)atan2(imag, real);
1682 else { /* deal with bad case */
1683 if (imag > 0) *buf++ = PI / 2.;
1684 else if (imag < 0) *buf++ = -PI / 2.;
1685 else *buf++ = 0;
1686 }
1687 } while (++i < n);
1688 } break;
1689
1690 default: {
1691 fprintf(stderr, "store_registers:illegal output form\n");
1692 exit(0);
1693 } break;
1694 }
1695 } break;
1696
1697 case DB: {
1698
1699 switch (buf_form) {
1700 case REAL: { /* real only */
1701 do {
1702 *buf++ = (float)20.*log10(fft_net->regr[i]);
1703 } while (++i < n);
1704 } break;
1705
1706 case IMAG: { /* imag only */
1707 do {
1708 *buf++ = (float)20.*log10(fft_net->regi[i]);
1709 } while (++i < n);
1710 } break;
1711
1712 case RECT: { /* real and imag */
1713 do {
1714 *buf++ = (float)20.*log10(fft_net->regr[i]);
1715 *buf++ = (float)20.*log10(fft_net->regi[i]);
1716 } while (++i < n);
1717 } break;
1718
1719 case MAG: { /* magnitude only */
1720 do {
1721 real = fft_net->regr[i];
1722 imag = fft_net->regi[i];
1723 *buf++ = (float)20.*log10(sqrt(real*real+imag*imag));
1724 } while (++i < n);
1725 } break;
1726
1727 case PHASE: { /* phase only */
1728 do {
1729 real = fft_net->regr[i];
1730 imag = fft_net->regi[i];
1731 if (real)
1732 *buf++ = (float)atan2(imag, real);
1733 else { /* deal with bad case */
1734 if (imag > 0) *buf++ = PI / 2.;
1735 else if (imag < 0) *buf++ = -PI / 2.;
1736 else *buf++ = 0;
1737 }
1738 } while (++i < n);
1739 } break;
1740
1741 case POLAR: { /* magnitude and phase */
1742 do {
1743 real = fft_net->regr[i];
1744 imag = fft_net->regi[i];
1745 *buf++ = (float)20.*log10(sqrt(real*real+imag*imag));
1746 if (real)
1747 *buf++ = (float)atan2(imag, real);
1748 else { /* deal with bad case */
1749 if (imag > 0) *buf++ = PI / 2.;
1750 else if (imag < 0) *buf++ = -PI / 2.;
1751 else *buf++ = 0;
1752 }
1753 } while (++i < n);
1754 } break;
1755
1756 default: {
1757 fprintf(stderr, "store_registers:illegal output form\n");
1758 exit(0);
1759 } break;
1760 }
1761 } break;
1762
1763 default: {
1764 fprintf(stderr, "store_registers:illegal output scale\n");
1765 exit(0);
1766 } break;
1767 }
1768}
1769
1770
1771
1772/*****************************************************************************/
1773/* COMPUTE TRANSFORMATION */
1774/*****************************************************************************/
1775
1776void compute_fft(FFT_NET *fft_net)
1777
1778
1779/* modifies: fft_net
1780 effects: Passes the values (already loaded) in the registers through
1781 the network, multiplying with appropriate coefficients at each
1782 stage. The fft result will be in the registers at the end of
1783 the computation. The direction of the transformation is indicated
1784 by the network flag 'direction'. The form of the computation is:
1785
1786 X(pn) = X(p) + C*X(q)
1787 X(qn) = X(p) - C*X(q)
1788
1789 where X(pn,qn) represents the output of the registers at each stage.
1790 The calculations are actually done in place. Register pointers are
1791 used to speed up the calculations.
1792
1793 Register and coefficient addresses involved in the calculations
1794 are stored sequentially and are accessed as such. fft_net->indexp,
1795 indexq contain pointers to the relevant addresses, and fft_net->coeffs,
1796 inv_coeffs points to the appropriate coefficients at each stage of the
1797 computation.
1798*/
1799
1800{
1801 SAMPLE **xpr, **xpi, **xqr, **xqi, *cr, *ci;
1802 int i;
1803 SAMPLE tpr, tpi, tqr, tqi;
1804 int bps = fft_net->bps;
1805 int cnt = bps * (fft_net->stages - 1);
1806
1807 /* predetermined register addresses and coefficients */
1808 xpr = fft_net->indexpr;
1809 xpi = fft_net->indexpi;
1810 xqr = fft_net->indexqr;
1811 xqi = fft_net->indexqi;
1812
1813 if (fft_net->direction==FORWARD) { /* FORWARD FFT coefficients */
1814 cr = fft_net->coeffr;
1815 ci = fft_net->coeffi;
1816 }
1817 else { /* INVERSE FFT coefficients */
1818 cr = fft_net->inv_coeffr;
1819 ci = fft_net->inv_coeffi;
1820 }
1821
1822 /* stage one coefficients are 1 + 0j so C*X(q)=X(q) */
1823 /* bps mults can be avoided */
1824
1825 for (i = 0; i < bps; i++) {
1826
1827 /* add X(p) and X(q) */
1828 tpr = **xpr + **xqr;
1829 tpi = **xpi + **xqi;
1830 tqr = **xpr - **xqr;
1831 tqi = **xpi - **xqi;
1832
1833 /* exchange register with temp */
1834 **xpr = tpr;
1835 **xpi = tpi;
1836 **xqr = tqr;
1837 **xqi = tqi;
1838
1839 /* next set of register for calculations: */
1840 xpr++; xpi++; xqr++; xqi++; cr++; ci++;
1841
1842 }
1843
1844 for (i = 0; i < cnt; i++) {
1845
1846 /* mult X(q) by coeff C */
1847 tqr = **xqr * *cr - **xqi * *ci;
1848 tqi = **xqr * *ci + **xqi * *cr;
1849
1850 /* exchange register with temp */
1851 **xqr = tqr;
1852 **xqi = tqi;
1853
1854 /* add X(p) and X(q) */
1855 tpr = **xpr + **xqr;
1856 tpi = **xpi + **xqi;
1857 tqr = **xpr - **xqr;
1858 tqi = **xpi - **xqi;
1859
1860 /* exchange register with temp */
1861 **xpr = tpr;
1862 **xpi = tpi;
1863 **xqr = tqr;
1864 **xqi = tqi;
1865 /* next set of register for calculations: */
1866 xpr++; xpi++; xqr++; xqi++; cr++; ci++;
1867 }
1868}
1869
1870
1871/****************************************************************************/
1872/* SUPPORT MODULES */
1873/****************************************************************************/
1874
1875void net_alloc(FFT_NET *fft_net)
1876
1877
1878/* effects: Allocates appropriate two dimensional arrays and assigns
1879 correct internal pointers.
1880*/
1881
1882{
1883
1884 int stages, bps, n;
1885
1886 n = fft_net->n;
1887 stages = fft_net->stages;
1888 bps = fft_net->bps;
1889
1890
1891 /* two dimensional arrays with elements stored sequentially */
1892
1893 fft_net->load_index = (int *)malloc(n * INT_SIZE);
1894 fft_net->regr = (SAMPLE *)malloc(n * SAMPLE_SIZE);
1895 fft_net->regi = (SAMPLE *)malloc(n * SAMPLE_SIZE);
1896 fft_net->coeffr = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
1897 fft_net->coeffi = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
1898 fft_net->inv_coeffr = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
1899 fft_net->inv_coeffi = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
1900 fft_net->indexpr = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
1901 fft_net->indexpi = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
1902 fft_net->indexqr = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
1903 fft_net->indexqi = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
1904
1905 /* one dimensional load window */
1906 fft_net->window = (SAMPLE *)malloc(n * SAMPLE_SIZE);
1907 fft_net->inv_window = (SAMPLE *)malloc(n * SAMPLE_SIZE);
1908}
1909
1910void net_dealloc(FFT_NET *fft_net)
1911
1912
1913/* effects: Deallocates given FFT network.
1914*/
1915
1916{
1917
1918 free((char *)fft_net->load_index);
1919 free((char *)fft_net->regr);
1920 free((char *)fft_net->regi);
1921 free((char *)fft_net->coeffr);
1922 free((char *)fft_net->coeffi);
1923 free((char *)fft_net->inv_coeffr);
1924 free((char *)fft_net->inv_coeffi);
1925 free((char *)fft_net->indexpr);
1926 free((char *)fft_net->indexpi);
1927 free((char *)fft_net->indexqr);
1928 free((char *)fft_net->indexqi);
1929 free((char *)fft_net->window);
1930 free((char *)fft_net->inv_window);
1931}
1932
1933
1934BOOL power_of_two(n)
1935
1936int n;
1937
1938/* effects: Returns TRUE if n is a power of two, otherwise FALSE.
1939*/
1940
1941{
1942 int i;
1943
1944 for (i = n; i > 1; i >>= 1)
1945 if (i & 1) return FALSE; /* more than one bit high */
1946 return TRUE;
1947}
1948
1949
1950void create_hanning(SAMPLE *window, int n, SAMPLE scale)
1951
1952/* effects: Fills the buffer window with a hanning window of the appropriate
1953 size scaled by scale.
1954*/
1955
1956{
1957 SAMPLE a, pi_div_n = PI/n;
1958 int k;
1959
1960 for (k=1; k <= n; k++) {
1961 a = sin(k * pi_div_n);
1962 *window++ = scale * a * a;
1963 }
1964}
1965
1966
1967void create_rectangular(SAMPLE *window, int n, SAMPLE scale)
1968
1969/* effects: Fills the buffer window with a rectangular window of the
1970 appropriate size of height scale.
1971*/
1972
1973{
1974 while (n--)
1975 *window++ = scale;
1976}
1977
1978
1979void short_to_float(short *short_buf, float *float_buf, int n)
1980
1981/* effects; Converts short_buf to floats and stores them in float_buf.
1982*/
1983
1984{
1985 while (n--) {
1986 *float_buf++ = (float)*short_buf++;
1987 }
1988}
1989
1990
1991/* here's the meat: */
1992
1993void pd_fft(float *buf, int npoints, int inverse)
1994{
1995 double renorm;
1996 float *fp, *fp2;
1997 int i;
1998 renorm = (inverse ? npoints : 1.);
1999 cfft((inverse ? INVERSE : FORWARD), npoints, RECTANGULAR,
2000 buf, RECT, LINEAR, buf, RECT, LINEAR, 0);
2001 for (i = npoints << 1, fp = buf; i--; fp++) *fp *= renorm;
2002}
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-05 14:07:47 +0000