This is fft_net.c in view mode; [Download] [Up]
/* Princeton version 1.1, 10/89 */
/*****************************************************************************/
/* */
/* Fast Fourier Transform */
/* FFT Network Interaction and Support Modules */
/* Kevin Peterson, MIT Media Lab, EMS */
/* UROP - Fall '86 */
/* REV: 6/12/87(KHP) - Generalized to one procedure call with typed I/O */
/* */
/*****************************************************************************/
/*****************************************************************************/
/* INCLUDE FILES */
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#include "fft.h"
#define BOOL int
#define TRUE 1
#define FALSE 0
/*****************************************************************************/
/* GLOBAL DECLARATIONS */
/*****************************************************************************/
static FFT_NET *firstnet;
/*****************************************************************************/
/* GENERALIZED FAST FOURIER TRANSFORM MODULE */
/*****************************************************************************/
fft(trnsfrm_dir, npnt, window,
source_buf, source_form, source_scale,
result_buf, result_form, result_scale)
int trnsfrm_dir, npnt, window;
float *source_buf, *result_buf;
int source_form, source_scale;
int result_form, result_scale;
/* 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
transformed 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.
*/
{
register FFT_NET *thisnet, *prevnet = NULL;
/* 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 != NULL && /* search linkedlist */
(thisnet->n != npnt || thisnet->wtype != window)) {
prevnet = thisnet;
thisnet = thisnet->next;
}
if (thisnet == NULL) { /* if didn't find right size & wtype */
thisnet = (FFT_NET *)malloc(sizeof(FFT_NET)); /* alloc new */
build_fft_network(thisnet, npnt, window); /* build net */
if (prevnet != NULL)
prevnet->next = thisnet; /* & add to list */
else firstnet = thisnet;
}
load_registers(thisnet, source_buf, source_form, source_scale,
trnsfrm_dir);
compute_fft(thisnet);
store_registers(thisnet, result_buf, result_form, result_scale);
}
fft_clear() /* Deallocate all preserved FFT networks.
Should use when finished with all computations. */
{
register FFT_NET *thisnet, *nextnet;
if (firstnet) {
thisnet=firstnet;
do {
nextnet = thisnet->next;
net_dealloc(thisnet);
free((char *)thisnet);
} while (thisnet = nextnet);
}
}
/*****************************************************************************/
/* NETWORK CONSTRUCTION */
/*****************************************************************************/
build_fft_network(fft_net, n, window_type)
FFT_NET *fft_net;
int n;
BOOL 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.
*/
{
register int cntr, i, j, s;
int stages, bps;
int **p, **q, *pp, *qp;
double two_pi_div_n = TWO_PI / n;
/* network physical constraints */
fft_net->n = n;
fft_net->bps = bps = n/2;
fft_net->wtype = window_type;
for (i = 0, j = n; j > 1; j >>= 1, i++);
fft_net->stages = stages = i;
fft_net->direction = FORWARD;
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;
register double *coeffpr = fft_net->coeffr;
register double *coeffpi = fft_net->coeffi;
register double *inv_coeffpr = fft_net->inv_coeffr;
register double *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);
*coeffpi = -sin(two_pi_div_n*power);
*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;
register double **indexpr = fft_net->indexpr;
register double **indexpi = fft_net->indexpi;
register double **indexqr = fft_net->indexqr;
register double **indexqi = fft_net->indexqi;
register double *regr = fft_net->regr;
register double *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 */ {
register int *load_index,*load_indexp;
register 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(temp_index);
for (s = 0; s < stages; s++) {
free(p[s]);free(q[s]);
}
free(p);free(q);
}
}
/*****************************************************************************/
/* REGISTER LOAD AND STORE */
/*****************************************************************************/
load_registers(fft_net, buf, buf_form, buf_scale, trnsfrm_dir)
FFT_NET *fft_net;
float *buf;
int buf_form, buf_scale, 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;
register double *window;
register 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]=(double)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]=(double)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]=(double)buf[index*2] * window[index];
fft_net->regi[i]=(double)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]=(double)(buf[index*2] * cos(buf[index*2+1]))
* window[index];
fft_net->regi[i]=(double)(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]=(double)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]=(double)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]=(double)pow(10., (1./20.)*buf[index*2])
* window[index];
fft_net->regi[i]=(double)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]=(double)(pow(10., (1./20.)*buf[index*2])
* cos(buf[index*2+1])) * window[index];
fft_net->regi[i]=(double)(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;
}
}
store_registers(fft_net, buf, buf_form, buf_scale)
FFT_NET *fft_net;
register float *buf;
int buf_form, buf_scale;
/* 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.
*/
{
register int i;
double 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)
*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)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 */
/*****************************************************************************/
compute_fft(fft_net)
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.
*/
{
register double **xpr, **xpi, **xqr, **xqi, *cr, *ci;
register int i;
double 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 if (fft_net->direction==INVERSE) { /* 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 */
/****************************************************************************/
net_alloc(fft_net)
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 = (double *)malloc(n * DOUBLE_SIZE);
fft_net->regi = (double *)malloc(n * DOUBLE_SIZE);
fft_net->coeffr = (double *)malloc(stages*bps*DOUBLE_SIZE);
fft_net->coeffi = (double *)malloc(stages*bps*DOUBLE_SIZE);
fft_net->inv_coeffr = (double *)malloc(stages*bps*DOUBLE_SIZE);
fft_net->inv_coeffi = (double *)malloc(stages*bps*DOUBLE_SIZE);
fft_net->indexpr = (double **)malloc(stages * bps * PNTR_SIZE);
fft_net->indexpi = (double **)malloc(stages * bps * PNTR_SIZE);
fft_net->indexqr = (double **)malloc(stages * bps * PNTR_SIZE);
fft_net->indexqi = (double **)malloc(stages * bps * PNTR_SIZE);
/* one dimensional load window */
fft_net->window = (double *)malloc(n * DOUBLE_SIZE);
fft_net->inv_window = (double *)malloc(n * DOUBLE_SIZE);
}
net_dealloc(fft_net)
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;
}
create_hanning(window, n, scale)
register double *window;
int n;
double scale;
/* effects: Fills the buffer window with a hanning window of the appropriate
size scaled by scale.
*/
{
double a, pi_div_n = PI/n;
register int k;
for (k=1; k <= n; k++) {
a = sin(k * pi_div_n);
*window++ = scale * a * a;
}
}
create_rectangular(window, n, scale)
register double *window;
register int n;
double scale;
/* effects: Fills the buffer window with a rectangular window of the
appropriate size of height scale.
*/
{
while (n--)
*window++ = scale;
}
short_to_float(short_buf, float_buf, n)
register short *short_buf;
register float *float_buf;
register int n;
/* effects; Converts short_buf to floats and stores them in float_buf.
*/
{
while (n--) {
*float_buf++ = (float)*short_buf++;
}
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.