summaryrefslogtreecommitdiff
path: root/src/sndobj/rfftw/cfft.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/sndobj/rfftw/cfft.c')
-rw-r--r--src/sndobj/rfftw/cfft.c939
1 files changed, 939 insertions, 0 deletions
diff --git a/src/sndobj/rfftw/cfft.c b/src/sndobj/rfftw/cfft.c
new file mode 100644
index 0000000..e87a0b9
--- /dev/null
+++ b/src/sndobj/rfftw/cfft.c
@@ -0,0 +1,939 @@
+/*
+ * Copyright (c) 1997-1999 Massachusetts Institute of Technology
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ */
+
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+
+/* some basic definitions */
+#ifndef BOOL
+#define BOOL int
+#define TRUE 1
+#define FALSE 0
+#endif
+
+#define SAMPLE float /* data type used in calculation */
+
+#define SHORT_SIZE sizeof(short)
+#define INT_SIZE sizeof(int)
+#define FLOAT_SIZE sizeof(float)
+#define SAMPLE_SIZE sizeof(SAMPLE)
+#define PNTR_SIZE sizeof(char *)
+
+#define PI 3.1415927
+#define TWO_PI 6.2831854
+
+/* type definitions for I/O buffers */
+#define REAL 0 /* real only */
+#define IMAG 2 /* imaginary only */
+#define RECT 8 /* real and imaginary */
+#define MAG 16 /* magnitude only */
+#define PHASE 32 /* phase only */
+#define POLAR 64 /* magnitude and phase*/
+
+/* scale definitions for I/O buffers */
+#define LINEAR 0
+#define DB 1 /* 20log10 */
+
+/* transform direction definition */
+#define FORWARD 1 /* Forward FFT */
+#define INVERSE 2 /* Inverse FFT */
+
+/* window type definitions */
+#define HANNING 1
+#define RECTANGULAR 0
+
+
+
+/* network structure definition */
+
+typedef struct Tfft_net {
+ int n;
+ int stages;
+ int bps;
+ int direction;
+ int window_type;
+ int *load_index;
+ SAMPLE *window, *inv_window;
+ SAMPLE *regr;
+ SAMPLE *regi;
+ SAMPLE **indexpr;
+ SAMPLE **indexpi;
+ SAMPLE **indexqr;
+ SAMPLE **indexqi;
+ SAMPLE *coeffr, *inv_coeffr;
+ SAMPLE *coeffi, *inv_coeffi;
+ struct Tfft_net *next;
+} FFT_NET;
+
+
+void cfft(int trnsfrm_dir, int npnt, int window,
+ float *source_buf, int source_form, int source_scale,
+ float *result_buf, int result_form, int result_scale, int debug);
+
+
+/*****************************************************************************/
+/* GLOBAL DECLARATIONS */
+/*****************************************************************************/
+
+static FFT_NET *firstnet;
+
+/* prototypes */
+
+void net_alloc(FFT_NET *fft_net);
+void net_dealloc(FFT_NET *fft_net);
+int power_of_two(int n);
+void create_hanning(SAMPLE *window, int n, SAMPLE scale);
+void create_rectangular(SAMPLE *window, int n, SAMPLE scale);
+void short_to_float(short *short_buf, float *float_buf, int n);
+void load_registers(FFT_NET *fft_net, float *buf, int buf_form,
+ int buf_scale, int trnsfrm_dir);
+void compute_fft(FFT_NET *fft_net);
+void store_registers(FFT_NET *fft_net, float *buf, int buf_form,
+ int buf_scale, int debug);
+void build_fft_network(FFT_NET *fft_net, int n, int window_type);
+
+/*****************************************************************************/
+/* GENERALIZED FAST FOURIER TRANSFORM MODULE */
+/*****************************************************************************/
+
+void cfft(int trnsfrm_dir, int npnt, int window,
+ float *source_buf, int source_form, int source_scale,
+ float *result_buf, int result_form, int result_scale, int debug)
+
+/* modifies: result_buf
+ effects: Computes npnt FFT specified by form, scale, and dir parameters.
+ Source samples (single precision float) are taken from soure_buf and
+ the transfrmd representation is stored in result_buf (single precision
+ float). The parameters are defined as follows:
+
+ trnsfrm_dir = FORWARD | INVERSE
+ npnt = 2^k for some any positive integer k
+ window = HANNING | RECTANGULAR
+ (RECT = real and imag parts, POLAR = magnitude and phase)
+ source_form = REAL | IMAG | RECT | POLAR
+ result_form = REAL | IMAG | RECT | MAG | PHASE | POLAR
+ xxxxxx_scale= LINEAR | DB ( 20log10 |mag| )
+
+ The input/output buffers are stored in a form appropriate to the type.
+ For example: REAL => {real, real, real ...},
+ MAG => {mag, mag, mag, ... },
+ RECT => {real, imag, real, imag, ... },
+ POLAR => {mag, phase, mag, phase, ... }.
+
+ To look at the magnitude (in db) of a 1024 point FFT of a real time
+ signal we have:
+
+ fft(FORWARD, 1024, RECTANGULAR, input, REAL, LINEAR, output, MAG, DB)
+
+ All possible input and output combinations are possible given the
+ choice of type and scale parameters.
+*/
+
+{
+ FFT_NET *thisnet = (FFT_NET *)0;
+ FFT_NET *lastnet = (FFT_NET *)0;
+
+ /* A linked list of fft networks of different sizes is maintained to
+ avoid building with every call. The network is built on the first
+ call but reused for subsequent calls requesting the same size
+ transformation.
+ */
+
+ thisnet=firstnet;
+ while (thisnet) {
+ if (!(thisnet->n == npnt) || !(thisnet->window_type == window)) {
+ /* current net doesn't match size or window type */
+ lastnet=thisnet;
+ thisnet=thisnet->next;
+ continue; /* keep looking */
+ }
+
+ else { /* network matches desired size */
+ load_registers(thisnet, source_buf, source_form, source_scale,
+ trnsfrm_dir);
+ compute_fft(thisnet); /* do transformation */
+ store_registers(thisnet, result_buf, result_form, result_scale,debug);
+ return;
+ }
+ }
+
+ /* none of existing networks match required size*/
+
+ if (lastnet) { /* add new network to end of list */
+ thisnet = (FFT_NET *)malloc(sizeof(FFT_NET)); /* allocate */
+ thisnet->next = 0;
+ lastnet->next = thisnet; /* add to end of list */
+ }
+ else { /* first network to be created */
+ thisnet=firstnet=(FFT_NET *)malloc(sizeof(FFT_NET)); /* alloc. */
+ thisnet->next = 0;
+ }
+
+ /* build new network and compute transformation */
+ build_fft_network(thisnet, npnt, window);
+ load_registers(thisnet, source_buf, source_form, source_scale,
+ trnsfrm_dir);
+ compute_fft(thisnet);
+ store_registers(thisnet, result_buf, result_form, result_scale,debug);
+ return;
+}
+
+void fft_clear(void)
+
+/* effects: Deallocates all preserved FFT networks. Should be used when
+ finished with all computations.
+*/
+
+{
+ FFT_NET *thisnet, *nextnet;
+
+ if (firstnet) {
+ thisnet=firstnet;
+ do {
+ nextnet = thisnet->next;
+ net_dealloc(thisnet);
+ free((char *)thisnet);
+ } while (thisnet = nextnet);
+ }
+}
+
+
+/*****************************************************************************/
+/* NETWORK CONSTRUCTION */
+/*****************************************************************************/
+
+void build_fft_network(FFT_NET *fft_net, int n, int window_type)
+
+
+/* modifies:fft_net
+ effects: Constructs the fft network as described in fft.h. Butterfly
+ coefficients, read/write indicies, bit reversed load indicies,
+ and array allocations are computed.
+*/
+
+{
+ int cntr, i, j, s;
+ int stages, bps;
+ int **p, **q, *pp, *qp;
+ SAMPLE two_pi_div_n = TWO_PI / n;
+
+
+ /* network definition */
+ fft_net->n = n;
+ fft_net->bps = bps = n/2;
+ for (i = 0, j = n; j > 1; j >>= 1, i++);
+ fft_net->stages = stages = i;
+ fft_net->direction = FORWARD;
+ fft_net->window_type = window_type;
+ fft_net->next = (FFT_NET *)0;
+
+ /* allocate registers, index, coefficient arrays */
+ net_alloc(fft_net);
+
+
+ /* create appropriate windows */
+ if (window_type==HANNING) {
+ create_hanning(fft_net->window, n, 1.);
+ create_hanning(fft_net->inv_window, n, 1./n);
+ }
+ else {
+ create_rectangular(fft_net->window, n, 1.);
+ create_rectangular(fft_net->inv_window, n, 1./n);
+ }
+
+
+ /* calculate butterfly coefficients */ {
+
+ int num_diff_coeffs, power_inc, power;
+ SAMPLE *coeffpr = fft_net->coeffr;
+ SAMPLE *coeffpi = fft_net->coeffi;
+ SAMPLE *inv_coeffpr = fft_net->inv_coeffr;
+ SAMPLE *inv_coeffpi = fft_net->inv_coeffi;
+
+ /* stage one coeffs are 1 + 0j */
+ for (i = 0; i < bps; i++) {
+ *coeffpr = *inv_coeffpr = 1.;
+ *coeffpi = *inv_coeffpi = 0.;
+ coeffpr++; inv_coeffpr++;
+ coeffpi++; inv_coeffpi++;
+ }
+
+ /* stage 2 to last stage coeffs need calculation */
+ /* (1<<r <=> 2^r */
+ for (s = 2; s <= stages; s++) {
+
+ num_diff_coeffs = n / (1 << (stages - s + 1));
+ power_inc = 1 << (stages -s);
+ cntr = 0;
+
+ for (i = bps/num_diff_coeffs; i > 0; i--) {
+
+ power = 0;
+
+ for (j = num_diff_coeffs; j > 0; j--) {
+ *coeffpr = cos(two_pi_div_n*power);
+ *inv_coeffpr = cos(two_pi_div_n*power);
+/* AAA change these signs */ *coeffpi = -sin(two_pi_div_n*power);
+/* change back */ *inv_coeffpi = sin(two_pi_div_n*power);
+ power += power_inc;
+ coeffpr++; inv_coeffpr++;
+ coeffpi++; inv_coeffpi++;
+ }
+ }
+ }
+ }
+
+ /* calculate network indicies: stage exchange indicies are
+ calculated and then used as offset values from the base
+ register locations. The final addresses are then stored in
+ fft_net.
+ */ {
+
+ int index, inc;
+ SAMPLE **indexpr = fft_net->indexpr;
+ SAMPLE **indexpi = fft_net->indexpi;
+ SAMPLE **indexqr = fft_net->indexqr;
+ SAMPLE **indexqi = fft_net->indexqi;
+ SAMPLE *regr = fft_net->regr;
+ SAMPLE *regi = fft_net->regi;
+
+
+ /* allocate temporary 2d stage exchange index, 1d temp
+ load index */
+ p = (int **)malloc(stages * PNTR_SIZE);
+ q = (int **)malloc(stages * PNTR_SIZE);
+
+ for (s = 0; s < stages; s++) {
+ p[s] = (int *)malloc(bps * INT_SIZE);
+ q[s] = (int *)malloc(bps * INT_SIZE);
+ }
+
+ /* calculate stage exchange indicies: */
+ for (s = 0; s < stages; s++) {
+ pp = p[s];
+ qp = q[s];
+ inc = 1 << s;
+ cntr = 1 << (stages-s-1);
+ i = j = index = 0;
+
+ do {
+ do {
+ qp[i] = index + inc;
+ pp[i++] = index++;
+ } while (++j < inc);
+ index = qp[i-1] + 1;
+ j = 0;
+ } while (--cntr);
+ }
+
+ /* compute actual address values using indicies as offsets */
+ for (s = 0; s < stages; s++) {
+ for (i = 0; i < bps; i++) {
+ *indexpr++ = regr + p[s][i];
+ *indexpi++ = regi + p[s][i];
+ *indexqr++ = regr + q[s][i];
+ *indexqi++ = regi + q[s][i];
+ }
+ }
+ }
+
+
+ /* calculate load indicies (bit reverse ordering) */
+ /* bit reverse ordering achieved by passing normal
+ order indicies backwards through the network */
+
+ /* init to normal order indicies */ {
+ int *load_index,*load_indexp;
+ int *temp_indexp, *temp_index;
+ temp_index=temp_indexp=(int *)malloc(n * INT_SIZE);
+
+ i = 0; j = n;
+ load_index = load_indexp = fft_net->load_index;
+
+ while (j--)
+ *load_indexp++ = i++;
+
+ /* pass indicies backwards through net */
+ for (s = stages - 1; s > 0; s--) {
+ pp = p[s];
+ qp = q[s];
+
+ for (i = 0; i < bps; i++) {
+ temp_index[pp[i]]=load_index[2*i];
+ temp_index[qp[i]]=load_index[2*i+1];
+ }
+ j = n;
+ load_indexp = load_index;
+ temp_indexp = temp_index;
+ while (j--)
+ *load_indexp++ = *temp_indexp++;
+ }
+
+ /* free all temporary arrays */
+ free((char *)temp_index);
+ for (s = 0; s < stages; s++) {
+ free((char *)p[s]);free((char *)q[s]);
+ }
+ free((char *)p);free((char *)q);
+ }
+}
+
+
+
+/*****************************************************************************/
+/* REGISTER LOAD AND STORE */
+/*****************************************************************************/
+
+void load_registers(FFT_NET *fft_net, float *buf, int buf_form,
+ int buf_scale, int trnsfrm_dir)
+
+/* effects: Multiplies the input buffer with the appropriate window and
+ stores the resulting values in the initial registers of the
+ network. Input buffer must contain values appropriate to form.
+ For RECT, the buffer contains real num. followed by imag num,
+ and for POLAR, it contains magnitude followed by phase. Pure
+ inputs are listed normally. Both LINEAR and DB scales are
+ interpreted.
+*/
+
+{
+ int *load_index = fft_net->load_index;
+ SAMPLE *window;
+ int index, i = 0, n = fft_net->n;
+
+ if (trnsfrm_dir==FORWARD) window = fft_net->window;
+ else if (trnsfrm_dir==INVERSE) window = fft_net->inv_window;
+ else {
+ fprintf(stderr, "load_registers:illegal transform direction\n");
+ exit(0);
+ }
+ fft_net->direction = trnsfrm_dir;
+
+ switch(buf_scale) {
+ case LINEAR: {
+
+ switch (buf_form) {
+ case REAL: { /* pure REAL */
+ while (i < fft_net->n) {
+ index = load_index[i];
+ fft_net->regr[i]=(SAMPLE)buf[index] * window[index];
+ fft_net->regi[i]=0.;
+ i++;
+ }
+ } break;
+
+ case IMAG: { /* pure IMAGinary */
+ while (i < fft_net->n) {
+ index = load_index[i];
+ fft_net->regr[i]=0;
+ fft_net->regi[i]=(SAMPLE)buf[index] * window[index];
+ i++;
+ }
+ } break;
+
+ case RECT: { /* both REAL and IMAGinary */
+ while (i < fft_net->n) {
+ index = load_index[i];
+ fft_net->regr[i]=(SAMPLE)buf[index*2] * window[index];
+ fft_net->regi[i]=(SAMPLE)buf[index*2+1] * window[index];
+ i++;
+ }
+ } break;
+
+ case POLAR: { /* magnitude followed by phase */
+ while (i < fft_net->n) {
+ index = load_index[i];
+ fft_net->regr[i]=(SAMPLE)(buf[index*2] * cos(buf[index*2+1]))
+ * window[index];
+ fft_net->regi[i]=(SAMPLE)(buf[index*2] * sin(buf[index*2+1]))
+ * window[index];
+ i++;
+ }
+ } break;
+
+ default: {
+ fprintf(stderr, "load_registers:illegal input form\n");
+ exit(0);
+ } break;
+ }
+ } break;
+
+ case DB: {
+
+ switch (buf_form) {
+ case REAL: { /* log pure REAL */
+ while (i < fft_net->n) {
+ index = load_index[i];
+ fft_net->regr[i]=(SAMPLE)pow(10., (1./20.)*buf[index])
+ * window[index]; /* window scaling after linearization */
+ fft_net->regi[i]=0.;
+ i++;
+ }
+ } break;
+
+ case IMAG: { /* log pure IMAGinary */
+ while (i < fft_net->n) {
+ index = load_index[i];
+ fft_net->regr[i]=0.;
+ fft_net->regi[i]=(SAMPLE)pow(10., (1./20.)*buf[index])
+ * window[index];
+ i++;
+ }
+ } break;
+
+ case RECT: { /* log REAL and log IMAGinary */
+ while (i < fft_net->n) {
+ index = load_index[i];
+ fft_net->regr[i]=(SAMPLE)pow(10., (1./20.)*buf[index*2])
+ * window[index];
+ fft_net->regi[i]=(SAMPLE)pow(10., (1./20.)*buf[index*2+1])
+ * window[index];
+ i++;
+ }
+ } break;
+
+ case POLAR: { /* log mag followed by phase */
+ while (i < fft_net->n) {
+ index = load_index[i];
+ fft_net->regr[i]=(SAMPLE)(pow(10., (1./20.)*buf[index*2])
+ * cos(buf[index*2+1])) * window[index];
+ fft_net->regi[i]=(SAMPLE)(pow(10., (1./20.)*buf[index*2])
+ * sin(buf[index*2+1])) * window[index];
+ i++;
+ }
+ } break;
+
+ default: {
+ fprintf(stderr, "load_registers:illegal input form\n");
+ exit(0);
+ } break;
+ }
+ } break;
+
+ default: {
+ fprintf(stderr, "load_registers:illegal input scale\n");
+ exit(0);
+ } break;
+ }
+}
+
+
+void store_registers(FFT_NET *fft_net, float *buf, int buf_form,
+ int buf_scale, int debug)
+
+/* modifies: buf
+ effects: Writes the final contents of the network registers into buf in
+ either linear or db scale, polar or rectangular form. If any of
+ the pure forms(REAL, IMAG, MAG, or PHASE) are used then only the
+ corresponding part of the registers is stored in buf.
+*/
+
+{
+ int i;
+ SAMPLE real, imag, mag, phase;
+ int n;
+
+ i = 0;
+ n = fft_net->n;
+
+ switch (buf_scale) {
+ case LINEAR: {
+
+ switch (buf_form) {
+ case REAL: { /* pure REAL */
+ do {
+ *buf++ = (float)fft_net->regr[i];
+ } while (++i < n);
+ } break;
+
+ case IMAG: { /* pure IMAGinary */
+ do {
+ *buf++ = (float)fft_net->regi[i];
+ } while (++i < n);
+ } break;
+
+ case RECT: { /* both REAL and IMAGinary */
+ do {
+ *buf++ = (float)fft_net->regr[i];
+ *buf++ = (float)fft_net->regi[i];
+ } while (++i < n);
+ } break;
+
+ case MAG: { /* magnitude only */
+ do {
+ real = fft_net->regr[i];
+ imag = fft_net->regi[i];
+ *buf++ = (float)sqrt(real*real+imag*imag);
+ } while (++i < n);
+ } break;
+
+ case PHASE: { /* phase only */
+ do {
+ real = fft_net->regr[i];
+ imag = fft_net->regi[i];
+ if (real > .00001)
+ *buf++ = (float)atan2(imag, real);
+ else { /* deal with bad case */
+ if (imag > 0){ *buf++ = PI / 2.;
+ if(debug) fprintf(stderr,"real=0 and imag > 0\n");}
+ else if (imag < 0){ *buf++ = -PI / 2.;
+ if(debug) fprintf(stderr,"real=0 and imag < 0\n");}
+ else { *buf++ = 0;
+ if(debug) fprintf(stderr,"real=0 and imag=0\n");}
+ }
+ } while (++i < n);
+ } break;
+
+ case POLAR: { /* magnitude and phase */
+ do {
+ real = fft_net->regr[i];
+ imag = fft_net->regi[i];
+ *buf++ = (float)sqrt(real*real+imag*imag);
+ if (real) /* a hack to avoid div by zero */
+ *buf++ = (float)atan2(imag, real);
+ else { /* deal with bad case */
+ if (imag > 0) *buf++ = PI / 2.;
+ else if (imag < 0) *buf++ = -PI / 2.;
+ else *buf++ = 0;
+ }
+ } while (++i < n);
+ } break;
+
+ default: {
+ fprintf(stderr, "store_registers:illegal output form\n");
+ exit(0);
+ } break;
+ }
+ } break;
+
+ case DB: {
+
+ switch (buf_form) {
+ case REAL: { /* real only */
+ do {
+ *buf++ = (float)20.*log10(fft_net->regr[i]);
+ } while (++i < n);
+ } break;
+
+ case IMAG: { /* imag only */
+ do {
+ *buf++ = (float)20.*log10(fft_net->regi[i]);
+ } while (++i < n);
+ } break;
+
+ case RECT: { /* real and imag */
+ do {
+ *buf++ = (float)20.*log10(fft_net->regr[i]);
+ *buf++ = (float)20.*log10(fft_net->regi[i]);
+ } while (++i < n);
+ } break;
+
+ case MAG: { /* magnitude only */
+ do {
+ real = fft_net->regr[i];
+ imag = fft_net->regi[i];
+ *buf++ = (float)20.*log10(sqrt(real*real+imag*imag));
+ } while (++i < n);
+ } break;
+
+ case PHASE: { /* phase only */
+ do {
+ real = fft_net->regr[i];
+ imag = fft_net->regi[i];
+ if (real)
+ *buf++ = (float)atan2(imag, real);
+ else { /* deal with bad case */
+ if (imag > 0) *buf++ = PI / 2.;
+ else if (imag < 0) *buf++ = -PI / 2.;
+ else *buf++ = 0;
+ }
+ } while (++i < n);
+ } break;
+
+ case POLAR: { /* magnitude and phase */
+ do {
+ real = fft_net->regr[i];
+ imag = fft_net->regi[i];
+ *buf++ = (float)20.*log10(sqrt(real*real+imag*imag));
+ if (real)
+ *buf++ = (float)atan2(imag, real);
+ else { /* deal with bad case */
+ if (imag > 0) *buf++ = PI / 2.;
+ else if (imag < 0) *buf++ = -PI / 2.;
+ else *buf++ = 0;
+ }
+ } while (++i < n);
+ } break;
+
+ default: {
+ fprintf(stderr, "store_registers:illegal output form\n");
+ exit(0);
+ } break;
+ }
+ } break;
+
+ default: {
+ fprintf(stderr, "store_registers:illegal output scale\n");
+ exit(0);
+ } break;
+ }
+}
+
+
+
+/*****************************************************************************/
+/* COMPUTE TRANSFORMATION */
+/*****************************************************************************/
+
+void compute_fft(FFT_NET *fft_net)
+
+
+/* modifies: fft_net
+ effects: Passes the values (already loaded) in the registers through
+ the network, multiplying with appropriate coefficients at each
+ stage. The fft result will be in the registers at the end of
+ the computation. The direction of the transformation is indicated
+ by the network flag 'direction'. The form of the computation is:
+
+ X(pn) = X(p) + C*X(q)
+ X(qn) = X(p) - C*X(q)
+
+ where X(pn,qn) represents the output of the registers at each stage.
+ The calculations are actually done in place. Register pointers are
+ used to speed up the calculations.
+
+ Register and coefficient addresses involved in the calculations
+ are stored sequentially and are accessed as such. fft_net->indexp,
+ indexq contain pointers to the relevant addresses, and fft_net->coeffs,
+ inv_coeffs points to the appropriate coefficients at each stage of the
+ computation.
+*/
+
+{
+ SAMPLE **xpr, **xpi, **xqr, **xqi, *cr, *ci;
+ int i;
+ SAMPLE tpr, tpi, tqr, tqi;
+ int bps = fft_net->bps;
+ int cnt = bps * (fft_net->stages - 1);
+
+ /* predetermined register addresses and coefficients */
+ xpr = fft_net->indexpr;
+ xpi = fft_net->indexpi;
+ xqr = fft_net->indexqr;
+ xqi = fft_net->indexqi;
+
+ if (fft_net->direction==FORWARD) { /* FORWARD FFT coefficients */
+ cr = fft_net->coeffr;
+ ci = fft_net->coeffi;
+ }
+ else { /* INVERSE FFT coefficients */
+ cr = fft_net->inv_coeffr;
+ ci = fft_net->inv_coeffi;
+ }
+
+ /* stage one coefficients are 1 + 0j so C*X(q)=X(q) */
+ /* bps mults can be avoided */
+
+ for (i = 0; i < bps; i++) {
+
+ /* add X(p) and X(q) */
+ tpr = **xpr + **xqr;
+ tpi = **xpi + **xqi;
+ tqr = **xpr - **xqr;
+ tqi = **xpi - **xqi;
+
+ /* exchange register with temp */
+ **xpr = tpr;
+ **xpi = tpi;
+ **xqr = tqr;
+ **xqi = tqi;
+
+ /* next set of register for calculations: */
+ xpr++; xpi++; xqr++; xqi++; cr++; ci++;
+
+ }
+
+ for (i = 0; i < cnt; i++) {
+
+ /* mult X(q) by coeff C */
+ tqr = **xqr * *cr - **xqi * *ci;
+ tqi = **xqr * *ci + **xqi * *cr;
+
+ /* exchange register with temp */
+ **xqr = tqr;
+ **xqi = tqi;
+
+ /* add X(p) and X(q) */
+ tpr = **xpr + **xqr;
+ tpi = **xpi + **xqi;
+ tqr = **xpr - **xqr;
+ tqi = **xpi - **xqi;
+
+ /* exchange register with temp */
+ **xpr = tpr;
+ **xpi = tpi;
+ **xqr = tqr;
+ **xqi = tqi;
+ /* next set of register for calculations: */
+ xpr++; xpi++; xqr++; xqi++; cr++; ci++;
+ }
+}
+
+
+/****************************************************************************/
+/* SUPPORT MODULES */
+/****************************************************************************/
+
+void net_alloc(FFT_NET *fft_net)
+
+
+/* effects: Allocates appropriate two dimensional arrays and assigns
+ correct internal pointers.
+*/
+
+{
+
+ int stages, bps, n;
+
+ n = fft_net->n;
+ stages = fft_net->stages;
+ bps = fft_net->bps;
+
+
+ /* two dimensional arrays with elements stored sequentially */
+
+ fft_net->load_index = (int *)malloc(n * INT_SIZE);
+ fft_net->regr = (SAMPLE *)malloc(n * SAMPLE_SIZE);
+ fft_net->regi = (SAMPLE *)malloc(n * SAMPLE_SIZE);
+ fft_net->coeffr = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
+ fft_net->coeffi = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
+ fft_net->inv_coeffr = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
+ fft_net->inv_coeffi = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
+ fft_net->indexpr = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
+ fft_net->indexpi = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
+ fft_net->indexqr = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
+ fft_net->indexqi = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
+
+ /* one dimensional load window */
+ fft_net->window = (SAMPLE *)malloc(n * SAMPLE_SIZE);
+ fft_net->inv_window = (SAMPLE *)malloc(n * SAMPLE_SIZE);
+}
+
+void net_dealloc(FFT_NET *fft_net)
+
+
+/* effects: Deallocates given FFT network.
+*/
+
+{
+
+ free((char *)fft_net->load_index);
+ free((char *)fft_net->regr);
+ free((char *)fft_net->regi);
+ free((char *)fft_net->coeffr);
+ free((char *)fft_net->coeffi);
+ free((char *)fft_net->inv_coeffr);
+ free((char *)fft_net->inv_coeffi);
+ free((char *)fft_net->indexpr);
+ free((char *)fft_net->indexpi);
+ free((char *)fft_net->indexqr);
+ free((char *)fft_net->indexqi);
+ free((char *)fft_net->window);
+ free((char *)fft_net->inv_window);
+}
+
+
+BOOL power_of_two(n)
+
+int n;
+
+/* effects: Returns TRUE if n is a power of two, otherwise FALSE.
+*/
+
+{
+ int i;
+
+ for (i = n; i > 1; i >>= 1)
+ if (i & 1) return FALSE; /* more than one bit high */
+ return TRUE;
+}
+
+
+void create_hanning(SAMPLE *window, int n, SAMPLE scale)
+
+/* effects: Fills the buffer window with a hanning window of the appropriate
+ size scaled by scale.
+*/
+
+{
+ SAMPLE a, pi_div_n = PI/n;
+ int k;
+
+ for (k=1; k <= n; k++) {
+ a = sin(k * pi_div_n);
+ *window++ = scale * a * a;
+ }
+}
+
+
+void create_rectangular(SAMPLE *window, int n, SAMPLE scale)
+
+/* effects: Fills the buffer window with a rectangular window of the
+ appropriate size of height scale.
+*/
+
+{
+ while (n--)
+ *window++ = scale;
+}
+
+
+void short_to_float(short *short_buf, float *float_buf, int n)
+
+/* effects; Converts short_buf to floats and stores them in float_buf.
+*/
+
+{
+ while (n--) {
+ *float_buf++ = (float)*short_buf++;
+ }
+}
+
+void ComplexFFT(float *buf, int npoints)
+{
+ float *fp, *fp2;
+ int i;
+ cfft(FORWARD, npoints, RECTANGULAR,
+ buf, RECT, LINEAR, buf, RECT, LINEAR, 0);
+}
+
+void ComplexIFFT(float *buf, int npoints)
+{
+ double renorm;
+ float *fp, *fp2;
+ int i;
+ renorm = (double) npoints;
+ cfft(INVERSE, npoints, RECTANGULAR,
+ buf, RECT, LINEAR, buf, RECT, LINEAR, 0);
+ for (i = npoints << 1, fp = buf; i--; fp++) *fp *= renorm;
+}
+
+
+
generated by cgit v1.2.3 (git 2.25.1) at 2025-04-05 10:54:33 +0000