dsp_fft32x32.c

TI c64x的FFT程序

        n10_n30 = _addsub(_lo(x_54), _lo(x_76));
        n11_n31 = _addsub(_hi(x_54), _hi(x_76));

        /*-------------------------------------------------------------*/
        /*  Store out the four outputs of 2 radix2 butterflies.        */
        /*-------------------------------------------------------------*/
        n0_dword0 = _itod(_hill(n01_n21), _hill(n00_n20));
        n0_dword1 = _itod(_hill(n11_n31), _hill(n10_n30));
        n0_dword2 = _itod(_loll(n01_n21), _loll(n00_n20));
        n0_dword3 = _itod(_loll(n11_n31), _loll(n10_n30));

        _amemd8(&y0[2*h2]) = n0_dword0;
        _amemd8(&y1[2*h2]) = n0_dword1;
        _amemd8(&y2[2*h2]) = n0_dword2;
        _amemd8(&y3[2*h2]) = n0_dword3;

        /*-------------------------------------------------------------*/
        /* Read in the next set of inputs from pointer "x2". These will*/
        /* produce outputs that are contiguous to the previous outputs.*/
        /*-------------------------------------------------------------*/
        x_98 = _amemd8(&x2[0]);
        x_BA = _amemd8(&x2[2]);
        x_DC = _amemd8(&x2[4]);
        x_FE = _amemd8(&x2[6]);

        x2 += 8;

        /*-------------------------------------------------------------*/
        /* Perform radix2 style decompositions.                        */
        /*    n02 = x_8 + x_A;    n03 = x_9 + x_B;                     */
        /*    n22 = x_8 - x_A;    n23 = x_9 - x_B;                     */
        /*    n12 = x_C + x_E;    n13 = x_D + x_F;                     */
        /*    n32 = x_C - x_E;    n33 = x_D - x_F;                     */
        /*-------------------------------------------------------------*/
        n02_n22 = _addsub(_lo(x_98), _lo(x_BA));
        n03_n23 = _addsub(_hi(x_98), _hi(x_BA));
        n12_n32 = _addsub(_lo(x_DC), _lo(x_FE));
        n13_n33 = _addsub(_hi(x_DC), _hi(x_FE));

        /*----------------------------------------------------------------*/
        /* Points that are read from succesive locations map              */
        /* to y, y[N/8], y[N/2],y[5N/8] in a radix4 scheme                */
        /*----------------------------------------------------------------*/
        n1_dword0 = _itod(_hill(n03_n23), _hill(n02_n22));
        n1_dword1 = _itod(_hill(n13_n33), _hill(n12_n32));
        n1_dword2 = _itod(_loll(n03_n23), _loll(n02_n22));
        n1_dword3 = _itod(_loll(n13_n33), _loll(n12_n32));

        _amemd8(&y0[2*h2 + 2]) = n1_dword0;
        _amemd8(&y1[2*h2 + 2]) = n1_dword1;
        _amemd8(&y2[2*h2 + 2]) = n1_dword2;
        _amemd8(&y3[2*h2 + 2]) = n1_dword3;

        /*---------------------------------------------------------------*/
        /* Increment j by "j0", if j is equal to n0, increment j by n0,  */
        /* that double reversal is avoided.                              */
        /*---------------------------------------------------------------*/
        j += j0;

        if (j == n0) {
            j  += n0;
            x0 += (int)npoints >> 1;
            x2 += (int)npoints >> 1;
        }
    }
}


void radix_4 (
    int * restrict ptr_x,
    int * restrict ptr_y,
    int npoints
)
{
    int n0, j0;
    int * restrict x2, * restrict x0;
    int * restrict y0, * restrict y1, * restrict y2, * restrict y3;
    int i, j, l1, h2;
    int xh0_2, xh1_2, xl0_2, xl1_2, xh0_3, xh1_3, xl0_3, xl1_3;
    int n02, n12, n22, n32, n03, n13, n23, n33;
    int xh0_0, xh1_0, xh0_1, xh1_1, xl0_0, xl1_0, xl0_1, xl1_1;
    double n0_dword0, n0_dword1, n0_dword2, n0_dword3;
    double n1_dword0, n1_dword1, n1_dword2, n1_dword3;
    double x_10, x_32, x_54, x_76, x_98, x_BA, x_DC, x_FE;
    long long xh0_0_xl0_0, xh1_0_xl1_0, xh0_1_xl0_1, xh1_1_xl1_1;
    long long xh0_2_xl0_2, xh1_2_xl1_2, xh0_3_xl0_3, xh1_3_xl1_3;
    long long n00_n20, n10_n30, n01_n21, n31_n11;
    long long n02_n22, n12_n32, n03_n23, n33_n13;

    /*-----------------------------------------------------------------*/
    /* The following code performs either a standard radix4 pass or a  */
    /* radix2 pass. Two pointers are used to access the input data.    */
    /* The input data is read "N/4" complex samples apart or "N/2"     */
    /* words apart using pointers "x0" and "x2". This produces out-    */
    /* puts that are 0, N/4, N/2, 3N/4 for a radix4 FFT, and 0, N/8    */
    /* N/2, 3N/8 for radix 2.                                          */
    /*-----------------------------------------------------------------*/
    y0 = ptr_y;
    y2 = ptr_y + (int)npoints;
    x0 = ptr_x;
    x2 = ptr_x + (int)(npoints >> 1);

    /*------------------------------------------------------------*/
    /* The pointers are set at the following locations.           */
    /*------------------------------------------------------------*/
    y1 = y0 + (int)(npoints >> 1);
    y3 = y2 + (int)(npoints >> 1);
    l1 = _norm(npoints) + 2;
    j0 = 4;
    n0 = npoints >> 2;

    /*-----------------------------------------------------------------*/
    /* The following code reads data indentically for either a radix 4 */
    /* or a radix 2 style decomposition. It writes out at different    */
    /* locations though. It checks if either half the points, or a     */
    /* quarter of the complex points have been exhausted to jump to    */
    /* pervent double reversal.                                        */
    /*-----------------------------------------------------------------*/
    j = 0;

    _nassert((int)(n0) % 4  == 0);
    _nassert((int)(x0) % 8 == 0);
    _nassert((int)(x2) % 8 == 0);
    _nassert((int)(y0) % 8 == 0);
    _nassert(npoints >= 16);
    #pragma MUST_ITERATE(2,,2);
    for (i = 0; i < npoints; i += 8) {
        /*-------------------------------------------------------------*/
        /* Digit reverse the index starting from 0. The increment to   */
        /* "j" is either by 4, or 8.                                   */
        /*-------------------------------------------------------------*/
        DIG_REV(j, l1, h2);

        /*-------------------------------------------------------------*/
        /* Read in the input data, from the first eight locations.     */
        /* These are transformed either as a radix4 or as a radix 2.   */
        /*-------------------------------------------------------------*/
        x_10 = _amemd8(&x0[0]);
        x_32 = _amemd8(&x0[2]);
        x_54 = _amemd8(&x0[4]);
        x_76 = _amemd8(&x0[6]);
        x0 += 8;

        xh0_0_xl0_0 = _addsub(_lo(x_10), _lo(x_54));
        xh1_0_xl1_0 = _addsub(_hi(x_10), _hi(x_54));
        xh0_1_xl0_1 = _addsub(_lo(x_32), _lo(x_76));
        xh1_1_xl1_1 = _addsub(_hi(x_32), _hi(x_76));

        xh0_0 = _hill(xh0_0_xl0_0);
        xh1_0 = _hill(xh1_0_xl1_0);
        xh0_1 = _hill(xh0_1_xl0_1);
        xh1_1 = _hill(xh1_1_xl1_1);
        xl0_0 = _loll(xh0_0_xl0_0);
        xl1_1 = _loll(xh1_1_xl1_1);
        xl1_0 = _loll(xh1_0_xl1_0);
        xl0_1 = _loll(xh0_1_xl0_1);

        n00_n20 = _addsub(xh0_0, xh0_1);
        n10_n30 = _addsub(xl0_0, xl1_1);
        n01_n21 = _addsub(xh1_0, xh1_1);
        n31_n11 = _addsub(xl1_0, xl0_1);

        /*-------------------------------------------------------------*/
        /*  Store out the four outputs of 1 radix4 butterfly or 2      */
        /*  radix2 butterflies.                                        */
        /*-------------------------------------------------------------*/
        n0_dword0 = _itod(_hill(n01_n21), _hill(n00_n20));
        n0_dword1 = _itod(_loll(n31_n11), _hill(n10_n30));
        n0_dword2 = _itod(_loll(n01_n21), _loll(n00_n20));
        n0_dword3 = _itod(_hill(n31_n11), _loll(n10_n30));

        _amemd8(&y0[2*h2]) = n0_dword0;
        _amemd8(&y1[2*h2]) = n0_dword1;
        _amemd8(&y2[2*h2]) = n0_dword2;
        _amemd8(&y3[2*h2]) = n0_dword3;

        /*-------------------------------------------------------------*/
        /* Read in the next set of inputs from pointer "x2". These will*/
        /* produce outputs that are contiguous to the previous outputs.*/
        /*-------------------------------------------------------------*/
        x_98 = _amemd8(&x2[0]);
        x_BA = _amemd8(&x2[2]);
        x_DC = _amemd8(&x2[4]);
        x_FE = _amemd8(&x2[6]);

        x2 += 8;

        /*-------------------------------------------------------------*/
        /* Perform radix4 style decompositions and overwrite results   */
        /* if it is dtermined that the radix to be used is radix 2.    */
        /*-------------------------------------------------------------*/
        xh0_2_xl0_2 = _addsub(_lo(x_98), _lo(x_DC));
        xh1_2_xl1_2 = _addsub(_hi(x_98), _hi(x_DC));
        xh0_3_xl0_3 = _addsub(_lo(x_BA), _lo(x_FE));
        xh1_3_xl1_3 = _addsub(_hi(x_BA), _hi(x_FE));

        xh0_2 = _hill(xh0_2_xl0_2);
        xh1_2 = _hill(xh1_2_xl1_2);
        xh0_3 = _hill(xh0_3_xl0_3);
        xh1_3 = _hill(xh1_3_xl1_3);
        xl0_2 = _loll(xh0_2_xl0_2);
        xl1_2 = _loll(xh1_2_xl1_2);
        xl0_3 = _loll(xh0_3_xl0_3);
        xl1_3 = _loll(xh1_3_xl1_3);

        n02_n22 = _addsub(xh0_2, xh0_3);
        n12_n32 = _addsub(xl0_2, xl1_3);
        n03_n23 = _addsub(xh1_2, xh1_3);
        n33_n13 = _addsub(xl1_2, xl0_3);

        n02 = _hill(n02_n22);
        n12 = _hill(n12_n32);
        n03 = _hill(n03_n23);
        n33 = _hill(n33_n13);
        n22 = _loll(n02_n22);
        n32 = _loll(n12_n32);
        n23 = _loll(n03_n23);
        n13 = _loll(n33_n13);

        /*----------------------------------------------------------------*/
        /* Points that are read from succesive locations map to y, y[N/4] */
        /* y[N/2], y[3N/4] in a radix4 scheme, y, y[N/8], y[N/2],y[5N/8]  */
        /*----------------------------------------------------------------*/
        n1_dword0 = _itod(n03, n02);
        n1_dword1 = _itod(n13, n12);
        n1_dword2 = _itod(n23, n22);
        n1_dword3 = _itod(n33, n32);

        _amemd8(&y0[2*h2 + 2]) = n1_dword0;
        _amemd8(&y1[2*h2 + 2]) = n1_dword1;
        _amemd8(&y2[2*h2 + 2]) = n1_dword2;
        _amemd8(&y3[2*h2 + 2]) = n1_dword3;

        /*---------------------------------------------------------------*/
        /* Increment j by "j0", if j is equal to n0, increment j by n0,  */
        /* that double reversal is avoided.                              */
        /*---------------------------------------------------------------*/
        j += j0;

        if (j == n0) {
            j += n0;
            x0 += (int)npoints >> 1;
            x2 += (int)npoints >> 1;
        }
    }
}

/* ======================================================================== */
/*  End of file: DSP_fft32x32.c                                             */
/* ------------------------------------------------------------------------ */
/*          Copyright (C) 2007 Texas Instruments, Incorporated.             */
/*                          All Rights Reserved.                            */
/* ======================================================================== */