hz_dit.c

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

#define PI 3.1415926535
/**************************************************************
*
*
*
*
* 
***************************************************************/
void fft_2dit(int n, double* x, double* y)
{
    double    a, e, xt, yt, c, s, xtt,ytt;    
    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;

    /* ------------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 = j - k;
            k = k / 2;
        }
   
        j = j + k;
    }

    //calculate FFT
    for(k = 1; k <= m; k++)
    {
        n2 = 1<<(k-1);
        n1 = 1<<k;
   
        e = -2*PI/n1;  
        for(j = 1; j <= n2; j++)
        {
            a = (j-1)*e;
            c = cos(a);
            s = sin(a);
            for(i = j; i <= n; i += n1)
            {
                q = i + n2;
         
                xtt = x[i];
                ytt = y[i];
         
                xt = x[q]*c - y[q]*s;
                yt = x[q]*s + y[q]*c;
         
                x[i] = xtt + xt;
                y[i] = ytt + yt;
         
                x[q] = xtt - xt;
                y[q] = ytt - yt;
            }
        }
    }    
}

/**************************************************************
*
*
*
*
* 
***************************************************************/
void ifft_2dit(int n, double* x, double* y)
{
    double    a, e, xt, yt, c, s, xtt,ytt;    
    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;

    /* ------------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 = j - k;
            k = k / 2;
        }
   
        j = j + k;
    }

    //calculate IFFT
    for(k = 1; k <= m; k++)
    {
        n2 = 1<<(k-1);
        n1 = 1<<k;
   
        e = 2*PI/n1;  
        for(j = 1; j <= n2; j++)
        {
            a = (j-1)*e;
            c = cos(a);
            s = sin(a);
            for(i = j; i <= n; i += n1)
            {
                q = i + n2;
         
                xtt = x[i];
                ytt = y[i];
         
                xt = x[q]*c - y[q]*s;
                yt = x[q]*s + y[q]*c;
         
                x[i] = xtt + xt;
                y[i] = ytt + yt;
         
                x[q] = xtt - xt;
                y[q] = ytt - yt;
            }
        }
    }    
    //each IFFT-ed data must be divided by 8
    for(i = 1; i <= n; i++)
    {
        x[i] = x[i]/n; 
        y[i] = y[i]/n;
    }
}