hz_dif.c

胡广书 数字信号处理 理论算法与实现 第二版

/*****************************************************************************
*
*    fft:        FFT implementation with Cooley-Tukey radix-2 DIF
*
*    Arguments:    n -- data number to be FFT-ed
*                x -- pointer to real      section of original data
*                y -- pointer to imaginary section of original data
*
*    Notes:        (original comments for the implementation)
*                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-2002 The MathWorks, Inc.
*                $Revision: 1.15 $  $Date: 2002/04/12 22:18:20 $
*****************************************************************************
*/
void fft_2dif(int n, double *x, double *y)
{
    double    a, e, xt, yt, c, s;    
    int        n1, n2, i, j, k, m, q;
    
    /*
     * 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)
     */
    
    m = (int)(log10(n) / log10(2));
    if( pow(2, m) != n)
    {
        printf("m must be 2's pow!\n");
    }
    
    /* --------------MAIN FFT LOOPS----------------------------- */
    
    /* Parameter adjustments */
    --y;
    --x;
    
    /* Function Body */
    n2 = n;
    for (k = 1; k <= m; ++k)                 
    {
        n1 = n2;
        n2 /= 2;
        e = (double)-6.283185307179586476925286766559005768394 / n1;
        a = 0.0;

        for (j = 1; j <= n2; ++j)             
        {
            c = cos(a);
            s = sin(a);
            a = j * e;                        
            
            for (i = j;                       
                 i <= n; 
                 i += n1
                )                             
            {
                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;
    }
}

/*****************************************************************************
*
*    ifft:        IFFT implementation with Cooley-Tukey radix-2 DIF
*
*    Arguments:    n -- data number to be IFFT-ed
*                x -- poiter to the real      section of FFT-ed data
*                y -- poiter to the imaginary section of FFT-ed data
*
*    Notes:        the only two difference between fft() and ifft() are:
*        (1) e = (double) 6.283185307179586476925286766559005768394 / n1, in 
fft()
*        changed to        
*            e = (double)-6.283185307179586476925286766559005768394 / n1, in 
iff();
*        (2)each IFFT-ed data must be divied by n; 
*****************************************************************************
*/
void ifft_2dif(int n, double *x, double *y)
{
    double    a, e, xt, yt, c, s;    
    int        n1, n2, i, j, k, m, q;
    
    /*
     * 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)
     */

    m = (int)(log10(n) / log10(2));
    if( pow(2, m) != n)
    {
        printf("m must be 2's pow!\n");
    }
    
    /* --------------MAIN FFT LOOPS----------------------------- */
    
    /* Parameter adjustments */
    --y;
    --x;
    
    /* Function Body */
    n2 = n;
    for (k = 1; k <= m; ++k)                 //step loops
    {
        n1 = n2;
        n2 /= 2;
        e = (double)6.283185307179586476925286766559005768394 / n1;
        a = 0.0;

        for (j = 1; j <= n2; ++j)             //butterfly loops within each gr
oup
        {
            c = cos(a);
            s = sin(a);
            a = j * e;                        //theta for calculating scale fa
ctor
            
            for (i = j;                        //group loops                
                i <= n; 
                 i += n1
                )                             //top    input of a butterfly
            {
                q = i + n2;                    //bottom input of a butterfly

                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;
    }

    //each IFFT-ed data must be divided by n
    for(i = 1; i <= n; i++)
    {
        x[i] = x[i]/n;    
        y[i] = y[i]/n;
    }    
}