diff options
Diffstat (limited to 'apps/plugins/pdbox/PDa/src/d_fftroutine.c')
-rw-r--r-- | apps/plugins/pdbox/PDa/src/d_fftroutine.c | 1000 |
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 | |||
1137 | typedef 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 | |||
1157 | void 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 | |||
1166 | static FFT_NET *firstnet; | ||
1167 | |||
1168 | /* prototypes */ | ||
1169 | |||
1170 | void net_alloc(FFT_NET *fft_net); | ||
1171 | void net_dealloc(FFT_NET *fft_net); | ||
1172 | int power_of_two(int n); | ||
1173 | void create_hanning(SAMPLE *window, int n, SAMPLE scale); | ||
1174 | void create_rectangular(SAMPLE *window, int n, SAMPLE scale); | ||
1175 | void short_to_float(short *short_buf, float *float_buf, int n); | ||
1176 | void load_registers(FFT_NET *fft_net, float *buf, int buf_form, | ||
1177 | int buf_scale, int trnsfrm_dir); | ||
1178 | void compute_fft(FFT_NET *fft_net); | ||
1179 | void store_registers(FFT_NET *fft_net, float *buf, int buf_form, | ||
1180 | int buf_scale, int debug); | ||
1181 | void build_fft_network(FFT_NET *fft_net, int n, int window_type); | ||
1182 | |||
1183 | /*****************************************************************************/ | ||
1184 | /* GENERALIZED FAST FOURIER TRANSFORM MODULE */ | ||
1185 | /*****************************************************************************/ | ||
1186 | |||
1187 | void 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 | |||
1269 | void 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 | |||
1293 | void 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 | |||
1475 | void 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 | |||
1609 | void 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 | |||
1776 | void 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 | |||
1875 | void 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 | |||
1910 | void 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 | |||
1934 | BOOL power_of_two(n) | ||
1935 | |||
1936 | int 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 | |||
1950 | void 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 | |||
1967 | void 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 | |||
1979 | void 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 | |||
1993 | void 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 | } | ||