?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fix pdbox makefile to actually take part in dependency generation
[kugel-rb.git] / apps / plugins / pdbox / PDa / src / d_fftroutine.c
blob5e8e25a77edc9fba6f5c5fd3399e94e8128fdd1e
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 #ifdef ROCKBOX
88 #include "plugin.h"
89 #include "../../pdbox.h"
90 #else /* ROCKBOX */
91 #include <stdio.h>
92 #include <math.h>
93 #include <stdlib.h>
94 #endif /* ROCKBOX */
95
96 /* the following is needed only to declare pd_fft() as exportable in MSW */
97 #include "m_pd.h"
98
99 /* some basic definitions */
100 #ifndef BOOL
101 #define BOOL int
102 #define TRUE 1
103 #define FALSE 0
104 #endif
105
106 #define SAMPLE float /* data type used in calculation */
107
108 #define SHORT_SIZE sizeof(short)
109 #define INT_SIZE sizeof(int)
110 #define FLOAT_SIZE sizeof(float)
111 #define SAMPLE_SIZE sizeof(SAMPLE)
112 #define PNTR_SIZE sizeof(char *)
113
114 #define PI 3.1415927
115 #define TWO_PI 6.2831854
116
117 /* type definitions for I/O buffers */
118 #define REAL 0 /* real only */
119 #define IMAG 2 /* imaginary only */
120 #define RECT 8 /* real and imaginary */
121 #define MAG 16 /* magnitude only */
122 #define PHASE 32 /* phase only */
123 #define POLAR 64 /* magnitude and phase*/
124
125 /* scale definitions for I/O buffers */
126 #define LINEAR 0
127 #define DB 1 /* 20log10 */
128
129 /* transform direction definition */
130 #define FORWARD 1 /* Forward FFT */
131 #define INVERSE 2 /* Inverse FFT */
132
133 /* window type definitions */
134 #define HANNING 1
135 #define RECTANGULAR 0
136
137
138
139 /* network structure definition */
140
141 typedef struct Tfft_net {
142 int n;
143 int stages;
144 int bps;
145 int direction;
146 int window_type;
147 int *load_index;
148 SAMPLE *window, *inv_window;
149 SAMPLE *regr;
150 SAMPLE *regi;
151 SAMPLE **indexpr;
152 SAMPLE **indexpi;
153 SAMPLE **indexqr;
154 SAMPLE **indexqi;
155 SAMPLE *coeffr, *inv_coeffr;
156 SAMPLE *coeffi, *inv_coeffi;
157 struct Tfft_net *next;
158 } FFT_NET;
159
160
161 void cfft(int trnsfrm_dir, int npnt, int window,
162 float *source_buf, int source_form, int source_scale,
163 float *result_buf, int result_form, int result_scale, int debug);
164
165
166 /*****************************************************************************/
167 /* GLOBAL DECLARATIONS */
168 /*****************************************************************************/
169
170 static FFT_NET *firstnet;
171
172 /* prototypes */
173
174 void net_alloc(FFT_NET *fft_net);
175 void net_dealloc(FFT_NET *fft_net);
176 int power_of_two(int n);
177 void create_hanning(SAMPLE *window, int n, SAMPLE scale);
178 void create_rectangular(SAMPLE *window, int n, SAMPLE scale);
179 void short_to_float(short *short_buf, float *float_buf, int n);
180 void load_registers(FFT_NET *fft_net, float *buf, int buf_form,
181 int buf_scale, int trnsfrm_dir);
182 void compute_fft(FFT_NET *fft_net);
183 void store_registers(FFT_NET *fft_net, float *buf, int buf_form,
184 int buf_scale, int debug);
185 void build_fft_network(FFT_NET *fft_net, int n, int window_type);
186
187 /*****************************************************************************/
188 /* GENERALIZED FAST FOURIER TRANSFORM MODULE */
189 /*****************************************************************************/
190
191 void cfft(int trnsfrm_dir, int npnt, int window,
192 float *source_buf, int source_form, int source_scale,
193 float *result_buf, int result_form, int result_scale, int debug)
194
195 /* modifies: result_buf
196 effects: Computes npnt FFT specified by form, scale, and dir parameters.
197 Source samples (single precision float) are taken from soure_buf and
198 the transfrmd representation is stored in result_buf (single precision
199 float). The parameters are defined as follows:
200
201 trnsfrm_dir = FORWARD | INVERSE
202 npnt = 2^k for some any positive integer k
203 window = HANNING | RECTANGULAR
204 (RECT = real and imag parts, POLAR = magnitude and phase)
205 source_form = REAL | IMAG | RECT | POLAR
206 result_form = REAL | IMAG | RECT | MAG | PHASE | POLAR
207 xxxxxx_scale= LINEAR | DB ( 20log10 |mag| )
208
209 The input/output buffers are stored in a form appropriate to the type.
210 For example: REAL => {real, real, real ...},
211 MAG => {mag, mag, mag, ... },
212 RECT => {real, imag, real, imag, ... },
213 POLAR => {mag, phase, mag, phase, ... }.
214
215 To look at the magnitude (in db) of a 1024 point FFT of a real time
216 signal we have:
217
218 fft(FORWARD, 1024, RECTANGULAR, input, REAL, LINEAR, output, MAG, DB)
219
220 All possible input and output combinations are possible given the
221 choice of type and scale parameters.
222 */
223
224 {
225 FFT_NET *thisnet = (FFT_NET *)0;
226 FFT_NET *lastnet = (FFT_NET *)0;
227
228 /* A linked list of fft networks of different sizes is maintained to
229 avoid building with every call. The network is built on the first
230 call but reused for subsequent calls requesting the same size
231 transformation.
232 */
233
234 thisnet=firstnet;
235 while (thisnet) {
236 if (!(thisnet->n == npnt) || !(thisnet->window_type == window)) {
237 /* current net doesn't match size or window type */
238 lastnet=thisnet;
239 thisnet=thisnet->next;
240 continue; /* keep looking */
241 }
242
243 else { /* network matches desired size */
244 load_registers(thisnet, source_buf, source_form, source_scale,
245 trnsfrm_dir);
246 compute_fft(thisnet); /* do transformation */
247 store_registers(thisnet, result_buf, result_form, result_scale,debug);
248 return;
249 }
250 }
251
252 /* none of existing networks match required size*/
253
254 if (lastnet) { /* add new network to end of list */
255 thisnet = (FFT_NET *)malloc(sizeof(FFT_NET)); /* allocate */
256 thisnet->next = 0;
257 lastnet->next = thisnet; /* add to end of list */
258 }
259 else { /* first network to be created */
260 thisnet=firstnet=(FFT_NET *)malloc(sizeof(FFT_NET)); /* alloc. */
261 thisnet->next = 0;
262 }
263
264 /* build new network and compute transformation */
265 build_fft_network(thisnet, npnt, window);
266 load_registers(thisnet, source_buf, source_form, source_scale,
267 trnsfrm_dir);
268 compute_fft(thisnet);
269 store_registers(thisnet, result_buf, result_form, result_scale,debug);
270 return;
271 }
272
273 void fft_clear(void)
274
275 /* effects: Deallocates all preserved FFT networks. Should be used when
276 finished with all computations.
277 */
278
279 {
280 FFT_NET *thisnet, *nextnet;
281
282 if (firstnet) {
283 thisnet=firstnet;
284 do {
285 nextnet = thisnet->next;
286 net_dealloc(thisnet);
287 free((char *)thisnet);
288 } while ((thisnet = nextnet));
289 }
290 }
291
292
293 /*****************************************************************************/
294 /* NETWORK CONSTRUCTION */
295 /*****************************************************************************/
296
297 void build_fft_network(FFT_NET *fft_net, int n, int window_type)
298
299
300 /* modifies:fft_net
301 effects: Constructs the fft network as described in fft.h. Butterfly
302 coefficients, read/write indicies, bit reversed load indicies,
303 and array allocations are computed.
304 */
305
306 {
307 int cntr, i, j, s;
308 int stages, bps;
309 int **p, **q, *pp, *qp;
310 SAMPLE two_pi_div_n = TWO_PI / n;
311
312
313 /* network definition */
314 fft_net->n = n;
315 fft_net->bps = bps = n/2;
316 for (i = 0, j = n; j > 1; j >>= 1, i++);
317 fft_net->stages = stages = i;
318 fft_net->direction = FORWARD;
319 fft_net->window_type = window_type;
320 fft_net->next = (FFT_NET *)0;
321
322 /* allocate registers, index, coefficient arrays */
323 net_alloc(fft_net);
324
325
326 /* create appropriate windows */
327 if (window_type==HANNING) {
328 create_hanning(fft_net->window, n, 1.);
329 create_hanning(fft_net->inv_window, n, 1./n);
330 }
331 else {
332 create_rectangular(fft_net->window, n, 1.);
333 create_rectangular(fft_net->inv_window, n, 1./n);
334 }
335
336
337 /* calculate butterfly coefficients */ {
338
339 int num_diff_coeffs, power_inc, power;
340 SAMPLE *coeffpr = fft_net->coeffr;
341 SAMPLE *coeffpi = fft_net->coeffi;
342 SAMPLE *inv_coeffpr = fft_net->inv_coeffr;
343 SAMPLE *inv_coeffpi = fft_net->inv_coeffi;
344
345 /* stage one coeffs are 1 + 0j */
346 for (i = 0; i < bps; i++) {
347 *coeffpr = *inv_coeffpr = 1.;
348 *coeffpi = *inv_coeffpi = 0.;
349 coeffpr++; inv_coeffpr++;
350 coeffpi++; inv_coeffpi++;
351 }
352
353 /* stage 2 to last stage coeffs need calculation */
354 /* (1<<r <=> 2^r */
355 for (s = 2; s <= stages; s++) {
356
357 num_diff_coeffs = n / (1 << (stages - s + 1));
358 power_inc = 1 << (stages -s);
359 cntr = 0;
360
361 for (i = bps/num_diff_coeffs; i > 0; i--) {
362
363 power = 0;
364
365 for (j = num_diff_coeffs; j > 0; j--) {
366 *coeffpr = cos(two_pi_div_n*power);
367 *inv_coeffpr = cos(two_pi_div_n*power);
368 /* AAA change these signs */ *coeffpi = -sin(two_pi_div_n*power);
369 /* change back */ *inv_coeffpi = sin(two_pi_div_n*power);
370 power += power_inc;
371 coeffpr++; inv_coeffpr++;
372 coeffpi++; inv_coeffpi++;
373 }
374 }
375 }
376 }
377
378 /* calculate network indicies: stage exchange indicies are
379 calculated and then used as offset values from the base
380 register locations. The final addresses are then stored in
381 fft_net.
382 */ {
383
384 int index, inc;
385 SAMPLE **indexpr = fft_net->indexpr;
386 SAMPLE **indexpi = fft_net->indexpi;
387 SAMPLE **indexqr = fft_net->indexqr;
388 SAMPLE **indexqi = fft_net->indexqi;
389 SAMPLE *regr = fft_net->regr;
390 SAMPLE *regi = fft_net->regi;
391
392
393 /* allocate temporary 2d stage exchange index, 1d temp
394 load index */
395 p = (int **)malloc(stages * PNTR_SIZE);
396 q = (int **)malloc(stages * PNTR_SIZE);
397
398 for (s = 0; s < stages; s++) {
399 p[s] = (int *)malloc(bps * INT_SIZE);
400 q[s] = (int *)malloc(bps * INT_SIZE);
401 }
402
403 /* calculate stage exchange indicies: */
404 for (s = 0; s < stages; s++) {
405 pp = p[s];
406 qp = q[s];
407 inc = 1 << s;
408 cntr = 1 << (stages-s-1);
409 i = j = index = 0;
410
411 do {
412 do {
413 qp[i] = index + inc;
414 pp[i++] = index++;
415 } while (++j < inc);
416 index = qp[i-1] + 1;
417 j = 0;
418 } while (--cntr);
419 }
420
421 /* compute actual address values using indicies as offsets */
422 for (s = 0; s < stages; s++) {
423 for (i = 0; i < bps; i++) {
424 *indexpr++ = regr + p[s][i];
425 *indexpi++ = regi + p[s][i];
426 *indexqr++ = regr + q[s][i];
427 *indexqi++ = regi + q[s][i];
428 }
429 }
430 }
431
432
433 /* calculate load indicies (bit reverse ordering) */
434 /* bit reverse ordering achieved by passing normal
435 order indicies backwards through the network */
436
437 /* init to normal order indicies */ {
438 int *load_index,*load_indexp;
439 int *temp_indexp, *temp_index;
440 temp_index=temp_indexp=(int *)malloc(n * INT_SIZE);
441
442 i = 0; j = n;
443 load_index = load_indexp = fft_net->load_index;
444
445 while (j--)
446 *load_indexp++ = i++;
447
448 /* pass indicies backwards through net */
449 for (s = stages - 1; s > 0; s--) {
450 pp = p[s];
451 qp = q[s];
452
453 for (i = 0; i < bps; i++) {
454 temp_index[pp[i]]=load_index[2*i];
455 temp_index[qp[i]]=load_index[2*i+1];
456 }
457 j = n;
458 load_indexp = load_index;
459 temp_indexp = temp_index;
460 while (j--)
461 *load_indexp++ = *temp_indexp++;
462 }
463
464 /* free all temporary arrays */
465 free((char *)temp_index);
466 for (s = 0; s < stages; s++) {
467 free((char *)p[s]);free((char *)q[s]);
468 }
469 free((char *)p);free((char *)q);
470 }
471 }
472
473
474 \f
475 /*****************************************************************************/
476 /* REGISTER LOAD AND STORE */
477 /*****************************************************************************/
478
479 void load_registers(FFT_NET *fft_net, float *buf, int buf_form,
480 int buf_scale, int trnsfrm_dir)
481
482 /* effects: Multiplies the input buffer with the appropriate window and
483 stores the resulting values in the initial registers of the
484 network. Input buffer must contain values appropriate to form.
485 For RECT, the buffer contains real num. followed by imag num,
486 and for POLAR, it contains magnitude followed by phase. Pure
487 inputs are listed normally. Both LINEAR and DB scales are
488 interpreted.
489 */
490
491 {
492 int *load_index = fft_net->load_index;
493 #ifdef ROCKBOX
494 SAMPLE *window = NULL;
495 int index, i = 0;
496 #else
497 SAMPLE *window;
498 int index, i = 0, n = fft_net->n;
499 #endif
500
501 if (trnsfrm_dir==FORWARD) window = fft_net->window;
502 else if (trnsfrm_dir==INVERSE) window = fft_net->inv_window;
503 else {
504 #ifndef ROCKBOX
505 fprintf(stderr, "load_registers:illegal transform direction\n");
506 exit(0);
507 #endif
508 }
509 fft_net->direction = trnsfrm_dir;
510
511 switch(buf_scale) {
512 case LINEAR: {
513
514 switch (buf_form) {
515 case REAL: { /* pure REAL */
516 while (i < fft_net->n) {
517 index = load_index[i];
518 fft_net->regr[i]=(SAMPLE)buf[index] * window[index];
519 fft_net->regi[i]=0.;
520 i++;
521 }
522 } break;
523
524 case IMAG: { /* pure IMAGinary */
525 while (i < fft_net->n) {
526 index = load_index[i];
527 fft_net->regr[i]=0;
528 fft_net->regi[i]=(SAMPLE)buf[index] * window[index];
529 i++;
530 }
531 } break;
532
533 case RECT: { /* both REAL and IMAGinary */
534 while (i < fft_net->n) {
535 index = load_index[i];
536 fft_net->regr[i]=(SAMPLE)buf[index*2] * window[index];
537 fft_net->regi[i]=(SAMPLE)buf[index*2+1] * window[index];
538 i++;
539 }
540 } break;
541
542 case POLAR: { /* magnitude followed by phase */
543 while (i < fft_net->n) {
544 index = load_index[i];
545 fft_net->regr[i]=(SAMPLE)(buf[index*2] * cos(buf[index*2+1]))
546 * window[index];
547 fft_net->regi[i]=(SAMPLE)(buf[index*2] * sin(buf[index*2+1]))
548 * window[index];
549 i++;
550 }
551 } break;
552
553 default: {
554 #ifndef ROCKBOX
555 fprintf(stderr, "load_registers:illegal input form\n");
556 exit(0);
557 #endif
558 } break;
559 }
560 } break;
561
562 case DB: {
563
564 switch (buf_form) {
565 case REAL: { /* log pure REAL */
566 while (i < fft_net->n) {
567 index = load_index[i];
568 fft_net->regr[i]=(SAMPLE)pow(10., (1./20.)*buf[index])
569 * window[index]; /* window scaling after linearization */
570 fft_net->regi[i]=0.;
571 i++;
572 }
573 } break;
574
575 case IMAG: { /* log pure IMAGinary */
576 while (i < fft_net->n) {
577 index = load_index[i];
578 fft_net->regr[i]=0.;
579 fft_net->regi[i]=(SAMPLE)pow(10., (1./20.)*buf[index])
580 * window[index];
581 i++;
582 }
583 } break;
584
585 case RECT: { /* log REAL and log IMAGinary */
586 while (i < fft_net->n) {
587 index = load_index[i];
588 fft_net->regr[i]=(SAMPLE)pow(10., (1./20.)*buf[index*2])
589 * window[index];
590 fft_net->regi[i]=(SAMPLE)pow(10., (1./20.)*buf[index*2+1])
591 * window[index];
592 i++;
593 }
594 } break;
595
596 case POLAR: { /* log mag followed by phase */
597 while (i < fft_net->n) {
598 index = load_index[i];
599 fft_net->regr[i]=(SAMPLE)(pow(10., (1./20.)*buf[index*2])
600 * cos(buf[index*2+1])) * window[index];
601 fft_net->regi[i]=(SAMPLE)(pow(10., (1./20.)*buf[index*2])
602 * sin(buf[index*2+1])) * window[index];
603 i++;
604 }
605 } break;
606
607 default: {
608 #ifndef ROCKBOX
609 fprintf(stderr, "load_registers:illegal input form\n");
610 exit(0);
611 #endif
612 } break;
613 }
614 } break;
615
616 default: {
617 #ifndef ROCKBOX
618 fprintf(stderr, "load_registers:illegal input scale\n");
619 exit(0);
620 #endif
621 } break;
622 }
623 }
624
625
626 void store_registers(FFT_NET *fft_net, float *buf, int buf_form,
627 int buf_scale, int debug)
628
629 /* modifies: buf
630 effects: Writes the final contents of the network registers into buf in
631 either linear or db scale, polar or rectangular form. If any of
632 the pure forms(REAL, IMAG, MAG, or PHASE) are used then only the
633 corresponding part of the registers is stored in buf.
634 */
635
636 {
637 #ifdef ROCKBOX
638 (void) debug;
639 #endif
640 int i;
641 #ifdef ROCKBOX
642 SAMPLE real, imag;
643 #else
644 SAMPLE real, imag, mag, phase;
645 #endif
646 int n;
647
648 i = 0;
649 n = fft_net->n;
650
651 switch (buf_scale) {
652 case LINEAR: {
653
654 switch (buf_form) {
655 case REAL: { /* pure REAL */
656 do {
657 *buf++ = (float)fft_net->regr[i];
658 } while (++i < n);
659 } break;
660
661 case IMAG: { /* pure IMAGinary */
662 do {
663 *buf++ = (float)fft_net->regi[i];
664 } while (++i < n);
665 } break;
666
667 case RECT: { /* both REAL and IMAGinary */
668 do {
669 *buf++ = (float)fft_net->regr[i];
670 *buf++ = (float)fft_net->regi[i];
671 } while (++i < n);
672 } break;
673
674 case MAG: { /* magnitude only */
675 do {
676 real = fft_net->regr[i];
677 imag = fft_net->regi[i];
678 *buf++ = (float)sqrt(real*real+imag*imag);
679 } while (++i < n);
680 } break;
681
682 case PHASE: { /* phase only */
683 do {
684 real = fft_net->regr[i];
685 imag = fft_net->regi[i];
686 if (real > .00001)
687 *buf++ = (float)atan2(imag, real);
688 else { /* deal with bad case */
689 #ifdef ROCKBOX
690 if (imag > 0){ *buf++ = PI / 2.;
691 }
692 else if (imag < 0){ *buf++ = -PI / 2.;
693 }
694 else { *buf++ = 0;
695 }
696 #else
697 if (imag > 0){ *buf++ = PI / 2.;
698 if(debug) fprintf(stderr,"real=0 and imag > 0\n");}
699 else if (imag < 0){ *buf++ = -PI / 2.;
700 if(debug) fprintf(stderr,"real=0 and imag < 0\n");}
701 else { *buf++ = 0;
702 if(debug) fprintf(stderr,"real=0 and imag=0\n");}
703 #endif
704 }
705 } while (++i < n);
706 } break;
707
708 case POLAR: { /* magnitude and phase */
709 do {
710 real = fft_net->regr[i];
711 imag = fft_net->regi[i];
712 *buf++ = (float)sqrt(real*real+imag*imag);
713 if (real) /* a hack to avoid div by zero */
714 *buf++ = (float)atan2(imag, real);
715 else { /* deal with bad case */
716 if (imag > 0) *buf++ = PI / 2.;
717 else if (imag < 0) *buf++ = -PI / 2.;
718 else *buf++ = 0;
719 }
720 } while (++i < n);
721 } break;
722
723 default: {
724 #ifndef ROCKBOX
725 fprintf(stderr, "store_registers:illegal output form\n");
726 exit(0);
727 #endif
728 } break;
729 }
730 } break;
731
732 case DB: {
733
734 switch (buf_form) {
735 case REAL: { /* real only */
736 do {
737 *buf++ = (float)20.*log10(fft_net->regr[i]);
738 } while (++i < n);
739 } break;
740
741 case IMAG: { /* imag only */
742 do {
743 *buf++ = (float)20.*log10(fft_net->regi[i]);
744 } while (++i < n);
745 } break;
746
747 case RECT: { /* real and imag */
748 do {
749 *buf++ = (float)20.*log10(fft_net->regr[i]);
750 *buf++ = (float)20.*log10(fft_net->regi[i]);
751 } while (++i < n);
752 } break;
753
754 case MAG: { /* magnitude only */
755 do {
756 real = fft_net->regr[i];
757 imag = fft_net->regi[i];
758 *buf++ = (float)20.*log10(sqrt(real*real+imag*imag));
759 } while (++i < n);
760 } break;
761
762 case PHASE: { /* phase only */
763 do {
764 real = fft_net->regr[i];
765 imag = fft_net->regi[i];
766 if (real)
767 *buf++ = (float)atan2(imag, real);
768 else { /* deal with bad case */
769 if (imag > 0) *buf++ = PI / 2.;
770 else if (imag < 0) *buf++ = -PI / 2.;
771 else *buf++ = 0;
772 }
773 } while (++i < n);
774 } break;
775
776 case POLAR: { /* magnitude and phase */
777 do {
778 real = fft_net->regr[i];
779 imag = fft_net->regi[i];
780 *buf++ = (float)20.*log10(sqrt(real*real+imag*imag));
781 if (real)
782 *buf++ = (float)atan2(imag, real);
783 else { /* deal with bad case */
784 if (imag > 0) *buf++ = PI / 2.;
785 else if (imag < 0) *buf++ = -PI / 2.;
786 else *buf++ = 0;
787 }
788 } while (++i < n);
789 } break;
790
791 default: {
792 #ifndef ROCKBOX
793 fprintf(stderr, "store_registers:illegal output form\n");
794 exit(0);
795 #endif
796 } break;
797 }
798 } break;
799
800 default: {
801 #ifndef ROCKBOX
802 fprintf(stderr, "store_registers:illegal output scale\n");
803 exit(0);
804 #endif
805 } break;
806 }
807 }
808
809
810 \f
811 /*****************************************************************************/
812 /* COMPUTE TRANSFORMATION */
813 /*****************************************************************************/
814
815 void compute_fft(FFT_NET *fft_net)
816
817
818 /* modifies: fft_net
819 effects: Passes the values (already loaded) in the registers through
820 the network, multiplying with appropriate coefficients at each
821 stage. The fft result will be in the registers at the end of
822 the computation. The direction of the transformation is indicated
823 by the network flag 'direction'. The form of the computation is:
824
825 X(pn) = X(p) + C*X(q)
826 X(qn) = X(p) - C*X(q)
827
828 where X(pn,qn) represents the output of the registers at each stage.
829 The calculations are actually done in place. Register pointers are
830 used to speed up the calculations.
831
832 Register and coefficient addresses involved in the calculations
833 are stored sequentially and are accessed as such. fft_net->indexp,
834 indexq contain pointers to the relevant addresses, and fft_net->coeffs,
835 inv_coeffs points to the appropriate coefficients at each stage of the
836 computation.
837 */
838
839 {
840 SAMPLE **xpr, **xpi, **xqr, **xqi, *cr, *ci;
841 int i;
842 SAMPLE tpr, tpi, tqr, tqi;
843 int bps = fft_net->bps;
844 int cnt = bps * (fft_net->stages - 1);
845
846 /* predetermined register addresses and coefficients */
847 xpr = fft_net->indexpr;
848 xpi = fft_net->indexpi;
849 xqr = fft_net->indexqr;
850 xqi = fft_net->indexqi;
851
852 if (fft_net->direction==FORWARD) { /* FORWARD FFT coefficients */
853 cr = fft_net->coeffr;
854 ci = fft_net->coeffi;
855 }
856 else { /* INVERSE FFT coefficients */
857 cr = fft_net->inv_coeffr;
858 ci = fft_net->inv_coeffi;
859 }
860
861 /* stage one coefficients are 1 + 0j so C*X(q)=X(q) */
862 /* bps mults can be avoided */
863
864 for (i = 0; i < bps; i++) {
865
866 /* add X(p) and X(q) */
867 tpr = **xpr + **xqr;
868 tpi = **xpi + **xqi;
869 tqr = **xpr - **xqr;
870 tqi = **xpi - **xqi;
871
872 /* exchange register with temp */
873 **xpr = tpr;
874 **xpi = tpi;
875 **xqr = tqr;
876 **xqi = tqi;
877
878 /* next set of register for calculations: */
879 xpr++; xpi++; xqr++; xqi++; cr++; ci++;
880
881 }
882
883 for (i = 0; i < cnt; i++) {
884
885 /* mult X(q) by coeff C */
886 tqr = **xqr * *cr - **xqi * *ci;
887 tqi = **xqr * *ci + **xqi * *cr;
888
889 /* exchange register with temp */
890 **xqr = tqr;
891 **xqi = tqi;
892
893 /* add X(p) and X(q) */
894 tpr = **xpr + **xqr;
895 tpi = **xpi + **xqi;
896 tqr = **xpr - **xqr;
897 tqi = **xpi - **xqi;
898
899 /* exchange register with temp */
900 **xpr = tpr;
901 **xpi = tpi;
902 **xqr = tqr;
903 **xqi = tqi;
904 /* next set of register for calculations: */
905 xpr++; xpi++; xqr++; xqi++; cr++; ci++;
906 }
907 }
908
909
910 /****************************************************************************/
911 /* SUPPORT MODULES */
912 /****************************************************************************/
913
914 void net_alloc(FFT_NET *fft_net)
915
916
917 /* effects: Allocates appropriate two dimensional arrays and assigns
918 correct internal pointers.
919 */
920
921 {
922
923 int stages, bps, n;
924
925 n = fft_net->n;
926 stages = fft_net->stages;
927 bps = fft_net->bps;
928
929
930 /* two dimensional arrays with elements stored sequentially */
931
932 fft_net->load_index = (int *)malloc(n * INT_SIZE);
933 fft_net->regr = (SAMPLE *)malloc(n * SAMPLE_SIZE);
934 fft_net->regi = (SAMPLE *)malloc(n * SAMPLE_SIZE);
935 fft_net->coeffr = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
936 fft_net->coeffi = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
937 fft_net->inv_coeffr = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
938 fft_net->inv_coeffi = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
939 fft_net->indexpr = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
940 fft_net->indexpi = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
941 fft_net->indexqr = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
942 fft_net->indexqi = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
943
944 /* one dimensional load window */
945 fft_net->window = (SAMPLE *)malloc(n * SAMPLE_SIZE);
946 fft_net->inv_window = (SAMPLE *)malloc(n * SAMPLE_SIZE);
947 }
948
949 void net_dealloc(FFT_NET *fft_net)
950
951
952 /* effects: Deallocates given FFT network.
953 */
954
955 {
956
957 free((char *)fft_net->load_index);
958 free((char *)fft_net->regr);
959 free((char *)fft_net->regi);
960 free((char *)fft_net->coeffr);
961 free((char *)fft_net->coeffi);
962 free((char *)fft_net->inv_coeffr);
963 free((char *)fft_net->inv_coeffi);
964 free((char *)fft_net->indexpr);
965 free((char *)fft_net->indexpi);
966 free((char *)fft_net->indexqr);
967 free((char *)fft_net->indexqi);
968 free((char *)fft_net->window);
969 free((char *)fft_net->inv_window);
970 }
971
972
973 BOOL power_of_two(n)
974
975 int n;
976
977 /* effects: Returns TRUE if n is a power of two, otherwise FALSE.
978 */
979
980 {
981 int i;
982
983 for (i = n; i > 1; i >>= 1)
984 if (i & 1) return FALSE; /* more than one bit high */
985 return TRUE;
986 }
987
988
989 void create_hanning(SAMPLE *window, int n, SAMPLE scale)
990
991 /* effects: Fills the buffer window with a hanning window of the appropriate
992 size scaled by scale.
993 */
994
995 {
996 SAMPLE a, pi_div_n = PI/n;
997 int k;
998
999 for (k=1; k <= n; k++) {
1000 a = sin(k * pi_div_n);
1001 *window++ = scale * a * a;
1002 }
1003 }
1004
1005
1006 void create_rectangular(SAMPLE *window, int n, SAMPLE scale)
1007
1008 /* effects: Fills the buffer window with a rectangular window of the
1009 appropriate size of height scale.
1010 */
1011
1012 {
1013 while (n--)
1014 *window++ = scale;
1015 }
1016
1017
1018 void short_to_float(short *short_buf, float *float_buf, int n)
1019
1020 /* effects; Converts short_buf to floats and stores them in float_buf.
1021 */
1022
1023 {
1024 while (n--) {
1025 *float_buf++ = (float)*short_buf++;
1026 }
1027 }
1028
1029
1030 /* here's the meat: */
1031
1032 void pd_fft(float *buf, int npoints, int inverse)
1033 {
1034 double renorm;
1035 #ifdef ROCKBOX
1036 float *fp;
1037 #else
1038 float *fp, *fp2;
1039 #endif
1040 int i;
1041 renorm = (inverse ? npoints : 1.);
1042 cfft((inverse ? INVERSE : FORWARD), npoints, RECTANGULAR,
1043 buf, RECT, LINEAR, buf, RECT, LINEAR, 0);
1044 for (i = npoints << 1, fp = buf; i--; fp++) *fp *= renorm;
1045 }
1046
My public Rockbox git repo