/* 
 * DSPFFT_RT - DSP Blockset 1-D FFT 
 *
 * Reference:
 * A COOLEY-TUKEY RADIX-2, DIF FFT PROGRAM
 *      COMPLEX INPUT DATA IN ARRAYS X AND Y
 *    C. S. BURRUS, RICE UNIVERSITY, SEPT 1983
 *
 *  Copyright (c) 1995-2000 The MathWorks, Inc.
 *  $Revision: 1.3 $  $Date: 2000/03/03 21:28:01 $
 */

#include "dspfft_rt.h"

void dspfft(const int_T n, creal_T *y)
{
    /* MAIN FFT LOOPS */
    
    /* Parameter adjustments */
    y--;
    
    /* Function Body */
    {
        int_T k;
        int_T n2 = n;
        int_T m;

        /*
         * Calculate m such that Nfft=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--;
        }

        for (k = 1; k <= m; ++k) {
            int_T  n1 = n2;
            int_T  j;
            const real_T e  = (real_T)6.283185307179586476925286766559005768394 / n1;
            real_T a  = 0.0;

	    n2 /= 2;

            for (j = 1; j <= n2; ++j) {
                int_T  i;
                const real_T c = cos(a);
	        const real_T s = sin(a);

	        a = j * e;

                for (i = j; (n1 < 0) ? i >= n : i <= n; i += n1) {
		    int_T   q = i + n2;
                    creal_T t;

                    t.re = y[i].re - y[q].re;
                    y[i].re += y[q].re;
        
                    t.im = y[i].im - y[q].im;
                    y[i].im += y[q].im;

                    y[q].re = c * t.re + s * t.im;
		    y[q].im = c * t.im - s * t.re;
	        }
	    }
        }
    }

    /* Bit Reverse */
    {
        int_T j = 1;
        int_T i;
    
        for (i = 1; i <= n-1; i++) {
	    if (i < j) {
                creal_T t = y[j]; y[j] = y[i]; y[i] = t;
	    }
            {
                int_T  k = n / 2;
	        while (k<j) {
	            j -= k;
	            k /= 2;
	        }
	        j += k;
            }
        }
    }
}

/* [EOF] dspfft_rt.c */