fftsg.cpp.svn-base

ffshow源码

                if (n > (nw << 2)) {
                        nw = n >> 2;
                        makewt(nw, ip, w);
                }
                nc = ip[1];
                if (n > (nc << 2)) {
                        nc = n >> 2;
                        makect(nc, ip, w + nw);
                }
                if (isgn >= 0) {
                        if (n > 4) {
                                cftfsub(n, a, ip + 2, nw, w);
                                rftfsub(n, a, nc, w + nw);
                        } else if (n == 4) {
                                cftfsub(n, a, ip + 2, nw, w);
                        }
                        xi = a[0] - a[1];
                        a[0] += a[1];
                        a[1] = xi;
                } else {
                        a[1] = 0.5f * (a[0] - a[1]);
                        a[0] -= a[1];
                        if (n > 4) {
                                rftbsub(n, a, nc, w + nw);
                                cftbsub(n, a, ip + 2, nw, w);
                        } else if (n == 4) {
                                cftbsub(n, a, ip + 2, nw, w);
                        }
                }
        }


        static void ddct(int n, int isgn, REAL *a, int *ip, REAL *w)
        {
                int j, nw, nc;
                REAL xr;
    
                nw = ip[0];
                if (n > (nw << 2)) {
                        nw = n >> 2;
                        makewt(nw, ip, w);
                }
                nc = ip[1];
                if (n > nc) {
                        nc = n;
                        makect(nc, ip, w + nw);
                }
                if (isgn < 0) {
                        xr = a[n - 1];
                        for (j = n - 2; j >= 2; j -= 2) {
                                a[j + 1] = a[j] - a[j - 1];
                                a[j] += a[j - 1];
                        }
                        a[1] = a[0] - xr;
                        a[0] += xr;
                        if (n > 4) {
                                rftbsub(n, a, nc, w + nw);
                                cftbsub(n, a, ip + 2, nw, w);
                        } else if (n == 4) {
                                cftbsub(n, a, ip + 2, nw, w);
                        }
                }
                dctsub(n, a, nc, w + nw);
                if (isgn >= 0) {
                        if (n > 4) {
                                cftfsub(n, a, ip + 2, nw, w);
                                rftfsub(n, a, nc, w + nw);
                        } else if (n == 4) {
                                cftfsub(n, a, ip + 2, nw, w);
                        }
                        xr = a[0] - a[1];
                        a[0] += a[1];
                        for (j = 2; j < n; j += 2) {
                                a[j - 1] = a[j] - a[j + 1];
                                a[j] += a[j + 1];
                        }
                        a[n - 1] = xr;
                }
        }


        static void ddst(int n, int isgn, REAL *a, int *ip, REAL *w)
        {
                int j, nw, nc;
                REAL xr;
    
                nw = ip[0];
                if (n > (nw << 2)) {
                        nw = n >> 2;
                        makewt(nw, ip, w);
                }
                nc = ip[1];
                if (n > nc) {
                        nc = n;
                        makect(nc, ip, w + nw);
                }
                if (isgn < 0) {
                        xr = a[n - 1];
                        for (j = n - 2; j >= 2; j -= 2) {
                                a[j + 1] = -a[j] - a[j - 1];
                                a[j] -= a[j - 1];
                        }
                        a[1] = a[0] + xr;
                        a[0] -= xr;
                        if (n > 4) {
                                rftbsub(n, a, nc, w + nw);
                                cftbsub(n, a, ip + 2, nw, w);
                        } else if (n == 4) {
                                cftbsub(n, a, ip + 2, nw, w);
                        }
                }
                dstsub(n, a, nc, w + nw);
                if (isgn >= 0) {
                        if (n > 4) {
                                cftfsub(n, a, ip + 2, nw, w);
                                rftfsub(n, a, nc, w + nw);
                        } else if (n == 4) {
                                cftfsub(n, a, ip + 2, nw, w);
                        }
                        xr = a[0] - a[1];
                        a[0] += a[1];
                        for (j = 2; j < n; j += 2) {
                                a[j - 1] = -a[j] - a[j + 1];
                                a[j] -= a[j + 1];
                        }
                        a[n - 1] = -xr;
                }
        }


        static void dfct(int n, REAL *a, REAL *t, int *ip, REAL *w)
        {
                int j, k, l, m, mh, nw, nc;
                REAL xr, xi, yr, yi;
    
                nw = ip[0];
                if (n > (nw << 3)) {
                        nw = n >> 3;
                        makewt(nw, ip, w);
                }
                nc = ip[1];
                if (n > (nc << 1)) {
                        nc = n >> 1;
                        makect(nc, ip, w + nw);
                }
                m = n >> 1;
                yi = a[m];
                xi = a[0] + a[n];
                a[0] -= a[n];
                t[0] = xi - yi;
                t[m] = xi + yi;
                if (n > 2) {
                        mh = m >> 1;
                        for (j = 1; j < mh; j++) {
                                k = m - j;
                                xr = a[j] - a[n - j];
                                xi = a[j] + a[n - j];
                                yr = a[k] - a[n - k];
                                yi = a[k] + a[n - k];
                                a[j] = xr;
                                a[k] = yr;
                                t[j] = xi - yi;
                                t[k] = xi + yi;
                        }
                        t[mh] = a[mh] + a[n - mh];
                        a[mh] -= a[n - mh];
                        dctsub(m, a, nc, w + nw);
                        if (m > 4) {
                                cftfsub(m, a, ip + 2, nw, w);
                                rftfsub(m, a, nc, w + nw);
                        } else if (m == 4) {
                                cftfsub(m, a, ip + 2, nw, w);
                        }
                        a[n - 1] = a[0] - a[1];
                        a[1] = a[0] + a[1];
                        for (j = m - 2; j >= 2; j -= 2) {
                                a[2 * j + 1] = a[j] + a[j + 1];
                                a[2 * j - 1] = a[j] - a[j + 1];
                        }
                        l = 2;
                        m = mh;
                        while (m >= 2) {
                                dctsub(m, t, nc, w + nw);
                                if (m > 4) {
                                        cftfsub(m, t, ip + 2, nw, w);
                                        rftfsub(m, t, nc, w + nw);
                                } else if (m == 4) {
                                        cftfsub(m, t, ip + 2, nw, w);
                                }
                                a[n - l] = t[0] - t[1];
                                a[l] = t[0] + t[1];
                                k = 0;
                                for (j = 2; j < m; j += 2) {
                                        k += l << 2;
                                        a[k - l] = t[j] - t[j + 1];
                                        a[k + l] = t[j] + t[j + 1];
                                }
                                l <<= 1;
                                mh = m >> 1;
                                for (j = 0; j < mh; j++) {
                                        k = m - j;
                                        t[j] = t[m + k] - t[m + j];
                                        t[k] = t[m + k] + t[m + j];
                                }
                                t[mh] = t[m + mh];
                                m = mh;
                        }
                        a[l] = t[0];
                        a[n] = t[2] - t[1];
                        a[0] = t[2] + t[1];
                } else {
                        a[1] = a[0];
                        a[2] = t[0];
                        a[0] = t[1];
                }
        }


        static void dfst(int n, REAL *a, REAL *t, int *ip, REAL *w)
        {
                int j, k, l, m, mh, nw, nc;
                REAL xr, xi, yr, yi;
    
                nw = ip[0];
                if (n > (nw << 3)) {
                        nw = n >> 3;
                        makewt(nw, ip, w);
                }
                nc = ip[1];
                if (n > (nc << 1)) {
                        nc = n >> 1;
                        makect(nc, ip, w + nw);
                }
                if (n > 2) {
                        m = n >> 1;
                        mh = m >> 1;
                        for (j = 1; j < mh; j++) {
                                k = m - j;
                                xr = a[j] + a[n - j];
                                xi = a[j] - a[n - j];
                                yr = a[k] + a[n - k];
                                yi = a[k] - a[n - k];
                                a[j] = xr;
                                a[k] = yr;
                                t[j] = xi + yi;
                                t[k] = xi - yi;
                        }
                        t[0] = a[mh] - a[n - mh];
                        a[mh] += a[n - mh];
                        a[0] = a[m];
                        dstsub(m, a, nc, w + nw);
                        if (m > 4) {
                                cftfsub(m, a, ip + 2, nw, w);
                                rftfsub(m, a, nc, w + nw);
                        } else if (m == 4) {
                                cftfsub(m, a, ip + 2, nw, w);
                        }
                        a[n - 1] = a[1] - a[0];
                        a[1] = a[0] + a[1];
                        for (j = m - 2; j >= 2; j -= 2) {
                                a[2 * j + 1] = a[j] - a[j + 1];
                                a[2 * j - 1] = -a[j] - a[j + 1];
                        }
                        l = 2;
                        m = mh;
                        while (m >= 2) {
                                dstsub(m, t, nc, w + nw);
                                if (m > 4) {
                                        cftfsub(m, t, ip + 2, nw, w);
                                        rftfsub(m, t, nc, w + nw);
                                } else if (m == 4) {
                                        cftfsub(m, t, ip + 2, nw, w);
                                }
                                a[n - l] = t[1] - t[0];
                                a[l] = t[0] + t[1];
                                k = 0;
                                for (j = 2; j < m; j += 2) {
                                        k += l << 2;
                                        a[k - l] = -t[j] - t[j + 1];
                                        a[k + l] = t[j] - t[j + 1];
                                }
                                l <<= 1;
                                mh = m >> 1;
                                for (j = 1; j < mh; j++) {
                                        k = m - j;
                                        t[j] = t[m + k] + t[m + j];
                                        t[k] = t[m + k] - t[m + j];
                                }
                                t[0] = t[m + mh];
                                m = mh;
                        }
                        a[l] = t[0];
                }
                a[0] = 0;
        }



        static void makewt(int nw, int *ip, REAL *w)
        {
                int j, nwh, nw0, nw1;
                REAL delta, wn4r, wk1r, wk1i, wk3r, wk3i;
    
                ip[0] = nw;
                ip[1] = 1;
                if (nw > 2) {
                        nwh = nw >> 1;