ftp.ee.lbl.gov/2014.05.ftp.ee.lbl.gov.tar/ftp.ee.lbl.gov/sst.tar.Z/sst.tar/sst/libfft.c

C/C++ Source or Header | 1990-01-10 | 4KB | 230 lines

/* libfft.c - fast Fourier transform library
**
** Copyright (C) 1989 by Jef Poskanzer.
**
** Permission to use, copy, modify, and distribute this software and its
** documentation for any purpose and without fee is hereby granted, provided
** that the above copyright notice appear in all copies and that both that
** copyright notice and this permission notice appear in supporting
** documentation.  This software is provided "as is" without express or
** implied warranty.
*/

#include <stdio.h>
#include <math.h>
#include "libfft.h"

#define MAXFFTSIZE 32768
#define LOG2_MAXFFTSIZE 15

static int bitreverse[MAXFFTSIZE], bits, n;
static complex w[MAXFFTSIZE];

/* initfft - initialize for fast Fourier transform
**
** b    power of two such that 2**nu = number of samples
*/
void
initfft( b )
int b;
    {
    register int n2, i, i0, j, k;
    float ws;
    complex u[32];
    register complex *wP, *v1;

    bits = b;
    if ( bits > LOG2_MAXFFTSIZE )
    {
    fprintf(
        stderr, "%d is too many bits, max is %d\n", bits, LOG2_MAXFFTSIZE );
    exit( 1 );
    }
    n = 1 << bits;

    /* Precompute the bitreverse array. */
    for ( i0 = ( 1 << bits ) - 1; i0 >= 0; --i0 )
    {
    k = 0;
    for ( j = 0; j < bits; ++j )
        {
        k *= 2;
        if ( i0 & ( 1 << j ) )
        k += 1;
        }
    bitreverse[i0] = k;
    }

    /* Precompute the u and w arrays - the "twiddle factors". */
    u[0].r = 0.;
    u[0].i = 1.;
    for ( i0 = 1; i0 < 32; ++i0 )
    {
    u[i0].r = sqrt( ( 1.0 + u[i0 - 1].r ) / 2.0 );
    u[i0].i = u[i0 - 1].i / ( 2.0 * u[i0].r );
    }
    n2 = n / 2;
    wP = w;
    for ( i0 = bits; --i0 >= 0; )
    {
    for ( k = 0; k < n2; ++k )
        {
        wP->r = 1.0;
        wP->i = 0.0;
        for ( i = k, v1 = &u[i0 - 1]; i; i >>= 1, --v1 )
        {
        if ( i & 1 )
            {
            ws    = wP->r * v1->r - wP->i * v1->i;
            wP->i = wP->r * v1->i + wP->i * v1->r;
            wP->r = ws;
            }
        }
        ++wP;
        }
    n2 /= 2;
    }
    }

/* fft - a fast Fourier transform routine
**
** x    data to be transformed
** inv  flag for inverse
*/

void
fft( x, inv )
complex x[];
int inv;
    {
    register complex *xP, *wP, *v1, *v2, *v3;
    register int i, i0, n2, j, k;

    n2 = n / 2;
    wP = w;
    for ( i0 = bits; --i0 >= 0; )
    {
    for ( k = 0; k < n2; ++k )
        {
        /* Transform the kth subseries. */
        v1 = &x[k];
        v3 = &x[n];
        while ( v1 < v3 )
        {
        register float tr, ti, z;

        v2 = v1 + n2;
        tr = v1->r + v2->r;
        ti = v1->i + v2->i;
        if ( inv )
            {
            z =     wP->r * (v1->r - v2->r) - wP->i * (v1->i - v2->i);
            v2->i = wP->r * (v1->i - v2->i) + wP->i * (v1->r - v2->r);
            }
        else
            {
            z =     wP->r * (v1->r - v2->r) + wP->i * (v1->i - v2->i);
            v2->i = wP->r * (v1->i - v2->i) - wP->i * (v1->r - v2->r);
            }
        v2->r = z;
        v1->r = tr;
        v1->i = ti;
        v1 = v2 + n2;
        }
        ++wP;
        }
    n2 /= 2;
    }

    /* Unscramble the data. */
    xP = x;
    for ( j = 0; j < n ; ++j )
    {
    k = bitreverse[j];
    if ( k > j )
        {
        complex t;

        t = *xP;
        *xP = x[k];
        x[k] = t;
        }
    ++xP;
    }

    /* Finally, divide each value by n, if this is the forward transform. */
    if ( ! inv )
    {
    register float f;

    f = 1.0 / n;
    xP = x;
    for ( j = 0; j < n ; ++j )
        {
        xP->r *= f;
        xP->i *= f;
        ++xP;
        }
    }
    }

/* fft_real_inverse - a fast Fourier transform routine
**
** x    data to be transformed
**
** Does an inverse transform without bothering to compute the imaginary
** values.  This is somewhat faster
*/

void
fft_real_inverse( x )
complex x[];
    {
    register complex *xP, *wP, *v1, *v2, *v3;
    register int i, i0, n2, j, k, wind;

    n2 = n / 2;
    wP = w;
    for ( i0 = bits; --i0 >= 0; )
    {
    for ( k = 0; k < n2; ++k )
        {
        /* Transform the kth subseries. */
        v1 = &x[k];
        v3 = &x[n];
        while ( v1 < v3 )
        {
        register float tr, ti, z;

        v2 = v1 + n2;
        tr = v1->r + v2->r;
        ti = v1->i + v2->i;
        z =     wP->r * ( v1->r - v2->r ) - wP->i * ( v1->i - v2->i );
        v2->i = wP->r * ( v1->i - v2->i ) + wP->i * ( v1->r - v2->r );
        v2->r = z;
        v1->r = tr;
        v1->i = ti;
        v1 = v2 + n2;
        }
        ++wP;
        }
    n2 /= 2;
    }

    /* Unscramble the data. */
    xP = x;
    for ( j = 0; j < n ; ++j )
    {
    k = bitreverse[j];
    if ( k > j )
        {
        register float tr;

        tr = xP->r;
        xP->r = x[k].r;
        x[k].r = tr;
        }
    ++xP;
    }
    }