summaryrefslogtreecommitdiff
path: root/apps/plugins/pdbox/PDa/src/d_fftroutine.c
diff options
context:
space:
mode:
Diffstat (limited to 'apps/plugins/pdbox/PDa/src/d_fftroutine.c')
-rw-r--r--apps/plugins/pdbox/PDa/src/d_fftroutine.c1000
1 files changed, 0 insertions, 1000 deletions
diff --git a/apps/plugins/pdbox/PDa/src/d_fftroutine.c b/apps/plugins/pdbox/PDa/src/d_fftroutine.c
index dde24fa358..41948132d7 100644
--- a/apps/plugins/pdbox/PDa/src/d_fftroutine.c
+++ b/apps/plugins/pdbox/PDa/src/d_fftroutine.c
@@ -999,1004 +999,4 @@ void pd_fft(float *buf, int npoints, int inverse)
999 buf, RECT, LINEAR, buf, RECT, LINEAR, 0); 999 buf, RECT, LINEAR, buf, RECT, LINEAR, 0);
1000 for (i = npoints << 1, fp = buf; i--; fp++) *fp *= renorm; 1000 for (i = npoints << 1, fp = buf; i--; fp++) *fp *= renorm;
1001} 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 1002
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-09-12 17:42:42 +0000