fft4f2d.c

2维fft程序

    void dctbsub(int n1, int n2, double **a, int nc, double *c);
    void dctfsub(int n1, int n2, double **a, int nc, double *c);
    int n, nw, nc, n1h, n2h, i, ix, ic, j, jx, jc;
    double xi;
    
    n = n1 << 1;
    if (n < n2) {
        n = n2;
    }
    nw = ip[0];
    if (n > (nw << 2)) {
        nw = n >> 2;
        makewt(nw, ip, w);
    }
    nc = ip[1];
    if (n1 > nc || n2 > nc) {
        if (n1 > n2) {
            nc = n1;
        } else {
            nc = n2;
        }
        makect(nc, ip, w + nw);
    }
    n1h = n1 >> 1;
    n2h = n2 >> 1;
    if (isgn >= 0) {
        for (i = 0; i <= n1 - 1; i++) {
            for (j = 1; j <= n2h - 1; j++) {
                jx = j << 1;
                t[i][jx] = a[i][j];
                t[i][jx + 1] = a[i][n2 - j];
            }
        }
        t[0][0] = a[0][0];
        t[0][1] = a[0][n2h];
        t[n1h][0] = a[n1h][0];
        t[n1h][1] = a[n1h][n2h];
        for (i = 1; i <= n1h - 1; i++) {
            ic = n1 - i;
            t[i][0] = a[i][0];
            t[ic][1] = a[i][n2h];
            t[i][1] = a[ic][0];
            t[ic][0] = a[ic][n2h];
        }
        dctfsub(n1, n2, t, nc, w + nw);
        if (n1 > 2) {
            bitrv2row(n1, n2, ip + 2, t);
        }
        cftfrow(n1, n2, t, w);
        for (i = 0; i <= n1 - 1; i++) {
            t[i][1] = 0.5 * (t[i][0] - t[i][1]);
            t[i][0] -= t[i][1];
        }
        if (n2 > 4) {
            rftfcol(n1, n2, t, nc, w + nw);
            bitrv2col(n1, n2, ip + 2, t);
        }
        cftfcol(n1, n2, t, w);
        for (i = 0; i <= n1h - 1; i++) {
            ix = i << 1;
            ic = n1 - 1 - i;
            for (j = 0; j <= n2h - 1; j++) {
                jx = j << 1;
                jc = n2 - 1 - j;
                a[ix][jx] = t[i][j];
                a[ix][jx + 1] = t[i][jc];
                a[ix + 1][jx] = t[ic][j];
                a[ix + 1][jx + 1] = t[ic][jc];
            }
        }
    } else {
        for (i = 0; i <= n1h - 1; i++) {
            ix = i << 1;
            ic = n1 - 1 - i;
            for (j = 0; j <= n2h - 1; j++) {
                jx = j << 1;
                jc = n2 - 1 - j;
                t[i][j] = a[ix][jx];
                t[i][jc] = a[ix][jx + 1];
                t[ic][j] = a[ix + 1][jx];
                t[ic][jc] = a[ix + 1][jx + 1];
            }
        }
        if (n2 > 4) {
            bitrv2col(n1, n2, ip + 2, t);
        }
        cftbcol(n1, n2, t, w);
        if (n2 > 4) {
            rftbcol(n1, n2, t, nc, w + nw);
        }
        for (i = 0; i <= n1 - 1; i++) {
            xi = t[i][0] - t[i][1];
            t[i][0] += t[i][1];
            t[i][1] = xi;
        }
        if (n1 > 2) {
            bitrv2row(n1, n2, ip + 2, t);
        }
        cftbrow(n1, n2, t, w);
        dctbsub(n1, n2, t, nc, w + nw);
        for (i = 0; i <= n1 - 1; i++) {
            for (j = 1; j <= n2h - 1; j++) {
                jx = j << 1;
                a[i][j] = t[i][jx];
                a[i][n2 - j] = t[i][jx + 1];
            }
        }
        a[0][0] = t[0][0];
        a[0][n2h] = t[0][1];
        a[n1h][0] = t[n1h][0];
        a[n1h][n2h] = t[n1h][1];
        for (i = 1; i <= n1h - 1; i++) {
            ic = n1 - i;
            a[i][0] = t[i][0];
            a[i][n2h] = t[ic][1];
            a[ic][0] = t[i][1];
            a[ic][n2h] = t[ic][0];
        }
    }
}


void ddst2d(int n1, int n2, int isgn, double **a, double **t, 
    int *ip, double *w)
{
    void makewt(int nw, int *ip, double *w);
    void makect(int nc, int *ip, double *c);
    void bitrv2col(int n1, int n, int *ip, double **a);
    void bitrv2row(int n, int n2, int *ip, double **a);
    void cftbcol(int n1, int n, double **a, double *w);
    void cftbrow(int n, int n2, double **a, double *w);
    void cftfcol(int n1, int n, double **a, double *w);
    void cftfrow(int n, int n2, double **a, double *w);
    void rftbcol(int n1, int n, double **a, int nc, double *c);
    void rftfcol(int n1, int n, double **a, int nc, double *c);
    void dstbsub(int n1, int n2, double **a, int nc, double *c);
    void dstfsub(int n1, int n2, double **a, int nc, double *c);
    int n, nw, nc, n1h, n2h, i, ix, ic, j, jx, jc;
    double xi;
    
    n = n1 << 1;
    if (n < n2) {
        n = n2;
    }
    nw = ip[0];
    if (n > (nw << 2)) {
        nw = n >> 2;
        makewt(nw, ip, w);
    }
    nc = ip[1];
    if (n1 > nc || n2 > nc) {
        if (n1 > n2) {
            nc = n1;
        } else {
            nc = n2;
        }
        makect(nc, ip, w + nw);
    }
    n1h = n1 >> 1;
    n2h = n2 >> 1;
    if (isgn >= 0) {
        for (i = 0; i <= n1 - 1; i++) {
            for (j = 1; j <= n2h - 1; j++) {
                jx = j << 1;
                t[i][jx] = a[i][j];
                t[i][jx + 1] = a[i][n2 - j];
            }
        }
        t[0][0] = a[0][0];
        t[0][1] = a[0][n2h];
        t[n1h][0] = a[n1h][0];
        t[n1h][1] = a[n1h][n2h];
        for (i = 1; i <= n1h - 1; i++) {
            ic = n1 - i;
            t[i][0] = a[i][0];
            t[ic][1] = a[i][n2h];
            t[i][1] = a[ic][0];
            t[ic][0] = a[ic][n2h];
        }
        dstfsub(n1, n2, t, nc, w + nw);
        if (n1 > 2) {
            bitrv2row(n1, n2, ip + 2, t);
        }
        cftfrow(n1, n2, t, w);
        for (i = 0; i <= n1 - 1; i++) {
            t[i][1] = 0.5 * (t[i][0] - t[i][1]);
            t[i][0] -= t[i][1];
        }
        if (n2 > 4) {
            rftfcol(n1, n2, t, nc, w + nw);
            bitrv2col(n1, n2, ip + 2, t);
        }
        cftfcol(n1, n2, t, w);
        for (i = 0; i <= n1h - 1; i++) {
            ix = i << 1;
            ic = n1 - 1 - i;
            for (j = 0; j <= n2h - 1; j++) {
                jx = j << 1;
                jc = n2 - 1 - j;
                a[ix][jx] = t[i][j];
                a[ix][jx + 1] = -t[i][jc];
                a[ix + 1][jx] = -t[ic][j];
                a[ix + 1][jx + 1] = t[ic][jc];
            }
        }
    } else {
        for (i = 0; i <= n1h - 1; i++) {
            ix = i << 1;
            ic = n1 - 1 - i;
            for (j = 0; j <= n2h - 1; j++) {
                jx = j << 1;
                jc = n2 - 1 - j;
                t[i][j] = a[ix][jx];
                t[i][jc] = -a[ix][jx + 1];
                t[ic][j] = -a[ix + 1][jx];
                t[ic][jc] = a[ix + 1][jx + 1];
            }
        }
        if (n2 > 4) {
            bitrv2col(n1, n2, ip + 2, t);
        }
        cftbcol(n1, n2, t, w);
        if (n2 > 4) {
            rftbcol(n1, n2, t, nc, w + nw);
        }
        for (i = 0; i <= n1 - 1; i++) {
            xi = t[i][0] - t[i][1];
            t[i][0] += t[i][1];
            t[i][1] = xi;
        }
        if (n1 > 2) {
            bitrv2row(n1, n2, ip + 2, t);
        }
        cftbrow(n1, n2, t, w);
        dstbsub(n1, n2, t, nc, w + nw);
        for (i = 0; i <= n1 - 1; i++) {
            for (j = 1; j <= n2h - 1; j++) {
                jx = j << 1;
                a[i][j] = t[i][jx];
                a[i][n2 - j] = t[i][jx + 1];
            }
        }
        a[0][0] = t[0][0];
        a[0][n2h] = t[0][1];
        a[n1h][0] = t[n1h][0];
        a[n1h][n2h] = t[n1h][1];
        for (i = 1; i <= n1h - 1; i++) {
            ic = n1 - i;
            a[i][0] = t[i][0];
            a[i][n2h] = t[ic][1];
            a[ic][0] = t[i][1];
            a[ic][n2h] = t[ic][0];
        }
    }
}


/* -------- initializing routines -------- */


#include <math.h>

void makewt(int nw, int *ip, double *w)
{
    void bitrv2(int n, int *ip, double *a);
    int nwh, j;
    double delta, x, y;
    
    ip[0] = nw;
    ip[1] = 1;
    if (nw > 2) {
        nwh = nw >> 1;
        delta = atan(1.0) / nwh;
        w[0] = 1;
        w[1] = 0;
        w[nwh] = cos(delta * nwh);
        w[nwh + 1] = w[nwh];
        for (j = 2; j <= nwh - 2; j += 2) {
            x = cos(delta * j);
            y = sin(delta * j);
            w[j] = x;
            w[j + 1] = y;
            w[nw - j] = y;
            w[nw - j + 1] = x;
        }
        bitrv2(nw, ip + 2, w);
    }
}


void makect(int nc, int *ip, double *c)
{
    int nch, j;
    double delta;
    
    ip[1] = nc;
    if (nc > 1) {
        nch = nc >> 1;
        delta = atan(1.0) / nch;
        c[0] = 0.5;
        c[nch] = 0.5 * cos(delta * nch);
        for (j = 1; j <= nch - 1; j++) {
            c[j] = 0.5 * cos(delta * j);
            c[nc - j] = 0.5 * sin(delta * j);
        }
    }
}


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


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


void bitrv2col(int n1, int n, int *ip, double **a)
{
    int i, j, j1, k, k1, l, m, m2;
    double xr, xi;
    
    ip[0] = 0;
    l = n;
    m = 1;
    while ((m << 2) < l) {
        l >>= 1;
        for (j = 0; j <= m - 1; j++) {
            ip[m + j] = ip[j] + l;
        }
        m <<= 1;
    }
    if ((m << 2) > l) {
        for (i = 0; i <= n1 - 1; i++) {
            for (k = 1; k <= m - 1; k++) {
                for (j = 0; j <= k - 1; j++) {
                    j1 = (j << 1) + ip[k];
                    k1 = (k << 1) + ip[j];
                    xr = a[i][j1];
                    xi = a[i][j1 + 1];
                    a[i][j1] = a[i][k1];
                    a[i][j1 + 1] = a[i][k1 + 1];
                    a[i][k1] = xr;
                    a[i][k1 + 1] = xi;
                }
            }
        }
    } else {
        m2 = m << 1;
        for (i = 0; i <= n1 - 1; i++) {
            for (k = 1; k <= m - 1; k++) {
                for (j = 0; j <= k - 1; j++) {
                    j1 = (j << 1) + ip[k];
                    k1 = (k << 1) + ip[j];
                    xr = a[i][j1];
                    xi = a[i][j1 + 1];
                    a[i][j1] = a[i][k1];
                    a[i][j1 + 1] = a[i][k1 + 1];
                    a[i][k1] = xr;
                    a[i][k1 + 1] = xi;
                    j1 += m2;
                    k1 += m2;
                    xr = a[i][j1];
                    xi = a[i][j1 + 1];
                    a[i][j1] = a[i][k1];
                    a[i][j1 + 1] = a[i][k1 + 1];
                    a[i][k1] = xr;
                    a[i][k1 + 1] = xi;
                }
            }
        }
    }
}


void bitrv2row(int n, int n2, int *ip, double **a)
{
    int i, j, j1, k, k1, l, m;