fftsg.cpp.svn-base

ffshow源码

                        delta = atanf(1.0f) / nwh;
                        wn4r = cosf(delta * nwh);
                        w[0] = 1;
                        w[1] = wn4r;
                        if (nwh >= 4) {
                                w[2] = 0.5f / cosf(delta * 2);
                                w[3] = 0.5f / cosf(delta * 6);
                        }
                        for (j = 4; j < nwh; j += 4) {
                                w[j] = cos(delta * j);
                                w[j + 1] = sinf(delta * j);
                                w[j + 2] = cosf(3 * delta * j);
                                w[j + 3] = sinf(3 * delta * j);
                        }
                        nw0 = 0;
                        while (nwh > 2) {
                                nw1 = nw0 + nwh;
                                nwh >>= 1;
                                w[nw1] = 1;
                                w[nw1 + 1] = wn4r;
                                if (nwh >= 4) {
                                        wk1r = w[nw0 + 4];
                                        wk3r = w[nw0 + 6];
                                        w[nw1 + 2] = 0.5f / wk1r;
                                        w[nw1 + 3] = 0.5f / wk3r;
                                }
                                for (j = 4; j < nwh; j += 4) {
                                        wk1r = w[nw0 + 2 * j];
                                        wk1i = w[nw0 + 2 * j + 1];
                                        wk3r = w[nw0 + 2 * j + 2];
                                        wk3i = w[nw0 + 2 * j + 3];
                                        w[nw1 + j] = wk1r;
                                        w[nw1 + j + 1] = wk1i;
                                        w[nw1 + j + 2] = wk3r;
                                        w[nw1 + j + 3] = wk3i;
                                }
                                nw0 = nw1;
                        }
                }
        }


        static void makect(int nc, int *ip, REAL *c)
        {
                int j, nch;
                REAL delta;
    
                ip[1] = nc;
                if (nc > 1) {
                        nch = nc >> 1;
                        delta = atanf(1.0f) / nch;
                        c[0] = cosf(delta * nch);
                        c[nch] = 0.5f * c[0];
                        for (j = 1; j < nch; j++) {
                                c[j] = 0.5f * cosf(delta * j);
                                c[nc - j] = 0.5f * sinf(delta * j);
                        }
                }
        }


        /* -------- child routines -------- */


        /* length of the recursive FFT mode */
        static const int CDFT_RECURSIVE_N=512;  /* <= (L1 cache size) / 16 */

        static void cftfsub(int n, REAL *a, int *ip, int nw, REAL *w)
        {
                int m;
    
                if (n > 32) {
                        m = n >> 2;
                        cftf1st(n, a, &w[nw - m]);
                        if (n > CDFT_RECURSIVE_N) {
                                cftrec1(m, a, nw, w);
                                cftrec2(m, &a[m], nw, w);
                                cftrec1(m, &a[2 * m], nw, w);
                                cftrec1(m, &a[3 * m], nw, w);
                        } else if (m > 32) {
                                cftexp1(n, a, nw, w);
                        } else {
                                cftfx41(n, a, nw, w);
                        }
                        bitrv2(n, ip, a);
                } else if (n > 8) {
                        if (n == 32) {
                                cftf161(a, &w[nw - 8]);
                                bitrv216(a);
                        } else {
                                cftf081(a, w);
                                bitrv208(a);
                        }
                } else if (n == 8) {
                        cftf040(a);
                } else if (n == 4) {
                        cftx020(a);
                }
        }


        static void cftbsub(int n, REAL *a, int *ip, int nw, REAL *w)
        {
                int m;
    
                if (n > 32) {
                        m = n >> 2;
                        cftb1st(n, a, &w[nw - m]);
                        if (n > CDFT_RECURSIVE_N) {
                                cftrec1(m, a, nw, w);
                                cftrec2(m, &a[m], nw, w);
                                cftrec1(m, &a[2 * m], nw, w);
                                cftrec1(m, &a[3 * m], nw, w);
                        } else if (m > 32) {
                                cftexp1(n, a, nw, w);
                        } else {
                                cftfx41(n, a, nw, w);
                        }
                        bitrv2conj(n, ip, a);
                } else if (n > 8) {
                        if (n == 32) {
                                cftf161(a, &w[nw - 8]);
                                bitrv216neg(a);
                        } else {
                                cftf081(a, w);
                                bitrv208neg(a);
                        }
                } else if (n == 8) {
                        cftb040(a);
                } else if (n == 4) {
                        cftx020(a);
                }
        }


        static void bitrv2(int n, int *ip, REAL *a)
        {
                int j, j1, k, k1, l, m, m2;
                REAL xr, xi, yr, yi;
    
                ip[0] = 0;
                l = n;
                m = 1;
                while ((m << 3) < l) {
                        l >>= 1;
                        for (j = 0; j < m; j++) {
                                ip[m + j] = ip[j] + l;
                        }
                        m <<= 1;
                }
                m2 = 2 * m;
                if ((m << 3) == l) {
                        for (k = 0; k < m; k++) {
                                for (j = 0; j < k; j++) {
                                        j1 = 2 * j + ip[k];
                                        k1 = 2 * k + ip[j];
                                        xr = a[j1];
                                        xi = a[j1 + 1];
                                        yr = a[k1];
                                        yi = a[k1 + 1];
                                        a[j1] = yr;
                                        a[j1 + 1] = yi;
                                        a[k1] = xr;
                                        a[k1 + 1] = xi;
                                        j1 += m2;
                                        k1 += 2 * m2;
                                        xr = a[j1];
                                        xi = a[j1 + 1];
                                        yr = a[k1];
                                        yi = a[k1 + 1];
                                        a[j1] = yr;
                                        a[j1 + 1] = yi;
                                        a[k1] = xr;
                                        a[k1 + 1] = xi;
                                        j1 += m2;
                                        k1 -= m2;
                                        xr = a[j1];
                                        xi = a[j1 + 1];
                                        yr = a[k1];
                                        yi = a[k1 + 1];
                                        a[j1] = yr;
                                        a[j1 + 1] = yi;
                                        a[k1] = xr;
                                        a[k1 + 1] = xi;
                                        j1 += m2;
                                        k1 += 2 * m2;
                                        xr = a[j1];
                                        xi = a[j1 + 1];
                                        yr = a[k1];
                                        yi = a[k1 + 1];
                                        a[j1] = yr;
                                        a[j1 + 1] = yi;
                                        a[k1] = xr;
                                        a[k1 + 1] = xi;
                                }
                                j1 = 2 * k + m2 + ip[k];
                                k1 = j1 + m2;
                                xr = a[j1];
                                xi = a[j1 + 1];
                                yr = a[k1];
                                yi = a[k1 + 1];
                                a[j1] = yr;
                                a[j1 + 1] = yi;
                                a[k1] = xr;
                                a[k1 + 1] = xi;
                        }
                } else {
                        for (k = 1; k < m; k++) {
                                for (j = 0; j < k; j++) {
                                        j1 = 2 * j + ip[k];
                                        k1 = 2 * k + ip[j];
                                        xr = a[j1];
                                        xi = a[j1 + 1];
                                        yr = a[k1];
                                        yi = a[k1 + 1];
                                        a[j1] = yr;
                                        a[j1 + 1] = yi;
                                        a[k1] = xr;
                                        a[k1 + 1] = xi;
                                        j1 += m2;
                                        k1 += m2;
                                        xr = a[j1];
                                        xi = a[j1 + 1];
                                        yr = a[k1];
                                        yi = a[k1 + 1];
                                        a[j1] = yr;
                                        a[j1 + 1] = yi;
                                        a[k1] = xr;
                                        a[k1 + 1] = xi;
                                }
                        }
                }
        }


        static void bitrv2conj(int n, int *ip, REAL *a)
        {
                int j, j1, k, k1, l, m, m2;
                REAL xr, xi, yr, yi;
    
                ip[0] = 0;
                l = n;
                m = 1;
                while ((m << 3) < l) {
                        l >>= 1;
                        for (j = 0; j < m; j++) {
                                ip[m + j] = ip[j] + l;
                        }
                        m <<= 1;
                }
                m2 = 2 * m;
                if ((m << 3) == l) {
                        for (k = 0; k < m; k++) {
                                for (j = 0; j < k; j++) {
                                        j1 = 2 * j + ip[k];
                                        k1 = 2 * k + ip[j];
                                        xr = a[j1];
                                        xi = -a[j1 + 1];
                                        yr = a[k1];
                                        yi = -a[k1 + 1];
                                        a[j1] = yr;
                                        a[j1 + 1] = yi;
                                        a[k1] = xr;
                                        a[k1 + 1] = xi;
                                        j1 += m2;
                                        k1 += 2 * m2;
                                        xr = a[j1];
                                        xi = -a[j1 + 1];
                                        yr = a[k1];
                                        yi = -a[k1 + 1];
                                        a[j1] = yr;
                                        a[j1 + 1] = yi;
                                        a[k1] = xr;
                                        a[k1 + 1] = xi;
                                        j1 += m2;
                                        k1 -= m2;
                                        xr = a[j1];
                                        xi = -a[j1 + 1];
                                        yr = a[k1];
                                        yi = -a[k1 + 1];
                                        a[j1] = yr;
                                        a[j1 + 1] = yi;
                                        a[k1] = xr;
                                        a[k1 + 1] = xi;
                                        j1 += m2;
                                        k1 += 2 * m2;
                                        xr = a[j1];
                                        xi = -a[j1 + 1];
                                        yr = a[k1];
                                        yi = -a[k1 + 1];
                                        a[j1] = yr;
                                        a[j1 + 1] = yi;
                                        a[k1] = xr;
                                        a[k1 + 1] = xi;
                                }
                                k1 = 2 * k + ip[k];
                                a[k1 + 1] = -a[k1 + 1];
                                j1 = k1 + m2;
                                k1 = j1 + m2;
                                xr = a[j1];
                                xi = -a[j1 + 1];
                                yr = a[k1];
                                yi = -a[k1 + 1];
                                a[j1] = yr;
                                a[j1 + 1] = yi;
                                a[k1] = xr;
                                a[k1 + 1] = xi;