/*
 *  fft_rt.c - Cooley-Tukey radix-2 DIF FFT
 *  Complex input data in arrays x and y
 *  C.S. Burrus, Rice University, Sept. 1983
 *
 *  Copyright 1995-2000 The MathWorks, Inc.
 *  $Revision: 1.10 $  $Date: 2000/03/22 20:43:24 $
 */

#include "fft_rt.h"
#include <math.h>

void fft(int_T n, real_T *x, real_T *y)
{
    real_T a, e, xt, yt;
    int_T  n1, n2, i, j, k, m;
    
    /*
     * Calculate m such that n=2^m
     *
     * NOTE: If frexp() == 1, then frexp does not conform
     * to the ANSI C spec of [0.5, 1)
     */
    if ( frexp((real_T)n, &m) != 1.0 ) {
        m--;
    }
    
    /* --------------MAIN FFT LOOPS----------------------------- */
    
    /* Parameter adjustments */
    --y;
    --x;
    
    /* Function Body */
    n2 = n;
    for (k = 1; k <= m; ++k) {
        n1 = n2;
        n2 /= 2;
        e = (real_T)6.283185307179586476925286766559005768394 / n1;
        a = 0.0;
        for (j = 1; j <= n2; ++j) {
            real_T c = cos(a);
            real_T s = sin(a);

            a = j * e;
            for (i = j; n1 < 0 ? i >= n : i <= n; i += n1) {
                int_T q = i + n2;

                xt = x[i] - x[q];
                x[i] += x[q];
                yt = y[i] - y[q];
                y[i] += y[q];
                x[q] = c * xt + s * yt;
                y[q] = c * yt - s * xt;
            }
        }
    }
    
    /* ------------DIGIT REVERSE COUNTER----------------- */
    
    j = 1;
    n1 = n - 1;
    for (i = 1; i <= n1; ++i) {
        if (i < j) {
            xt = x[j]; x[j] = x[i];     x[i] = xt;
            xt = y[j]; y[j] = y[i];     y[i] = xt;
        }
        k = n / 2;
        while (k<j) {
            j -= k;
            k /= 2;
        }
        j += k;
    }
}

/* [EOF] fft_rt.c */