/* * 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 SCCSstatic 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;        }    }}