getfft.c

speech signal process tools

/*
 * This material contains proprietary software of Entropic Speech, Inc.
 * Any reproduction, distribution, or publication without the prior
 * written permission of Entropic Speech, Inc. is strictly prohibited.
 * Any public distribution of copies of this work authorized in writing by
 * Entropic Speech, Inc. must bear the notice
 *
 *    "Copyright (c) 1986, 1987 Entropic Speech, Inc.; All rights reserved"
 *
 */
/*--------------------------------------------------------------+
|								|
|  get_fft -- compute discrete Fourier transform		|
|								|
|  by Shankar Narayan, EPI; some modifications by Rod Johnson.	|
|								|
|  Module:  getfft.c						|
|								|
+--------------------------------------------------------------*/
#ifdef SCCS
static char *sccs_id = "@(#)getfft.c	1.1 6/11/87";
#endif

#include <stdio.h>
typedef struct complex {
    double  real;
    double  imag;
}       COMPLEX;

COMPLEX cadd(), csub(), cmult();
double  sin(), cos(), log();

extern  debug_level;

#define ROUND(X) ((X)<0 ? -(int)(0.5-(X)) : (int)(0.5+(X)))

	/* Fourier transform of array x[n] */ 

get_fft_comp(x, n)
    COMPLEX x[];
    int     n;
{
    COMPLEX t, u[15], v, w;
    double  twopi = 6.283185307179586477;
    double  ang;
    int     i, l, m, p, q, s, sby2, stage;

    m = ROUND(log((double) n) / log(2.0));
    ang = twopi/n; w.real = cos(ang); w.imag = -sin(ang);

    u[1] = w;
    for (i=2; i<= m; i++)
        u[i] = cmult(u[i-1], u[i-1]);

    sby2 = n;
    for (stage = 1; stage <= m; stage++)
    {
        v = u[stage];
        w.real = 1.0; w.imag = 0.0;

        s = sby2;
        sby2 = s/2;
        for (l = 0; l < sby2; l++)
        {
            for (i = 0; i< n; i += s)
            {
                p = i + l;
                q = p + sby2;
                t = cadd(x[p], x[q]);
                x[q] = csub(x[p], x[q]);
                x[q] = cmult(x[q], w);
                x[p] = t;
            }
            w = cmult(w, v);
        }
    }
    /* perform bit-reversal operation */
    for (i = 0; i < n; i++)
    {
        p = i;
        q = 0;
        for (l = 0; l < m; l++)
        {
            q = 2*q + p%2;
            p = p/2;
        }
        if (i < q)
        {
            t = x[i];
            x[i] = x[q];
            x[q] = t;
        }
    }
}