fft4f2d.c

2维fft程序

                x3r = a[j2][i] - a[j3][i];
                x3i = a[j2][i + 1] - a[j3][i + 1];
                a[j][i] = x0r + x2r;
                a[j][i + 1] = x0i + x2i;
                a[j2][i] = x0r - x2r;
                a[j2][i + 1] = x0i - x2i;
                a[j1][i] = x1r + x3i;
                a[j1][i + 1] = x1i - x3r;
                a[j3][i] = x1r - x3i;
                a[j3][i + 1] = x1i + x3r;
            }
        }
        if (m < n) {
            wk1r = w[2];
            for (j = m; j <= l + m - 1; j++) {
                j1 = j + l;
                j2 = j1 + l;
                j3 = j2 + l;
                for (i = 0; i <= n2 - 2; i += 2) {
                    x0r = a[j][i] + a[j1][i];
                    x0i = a[j][i + 1] + a[j1][i + 1];
                    x1r = a[j][i] - a[j1][i];
                    x1i = a[j][i + 1] - a[j1][i + 1];
                    x2r = a[j2][i] + a[j3][i];
                    x2i = a[j2][i + 1] + a[j3][i + 1];
                    x3r = a[j2][i] - a[j3][i];
                    x3i = a[j2][i + 1] - a[j3][i + 1];
                    a[j][i] = x0r + x2r;
                    a[j][i + 1] = x0i + x2i;
                    a[j2][i] = x0i - x2i;
                    a[j2][i + 1] = x2r - x0r;
                    x0r = x1r + x3i;
                    x0i = x1i - x3r;
                    a[j1][i] = wk1r * (x0i + x0r);
                    a[j1][i + 1] = wk1r * (x0i - x0r);
                    x0r = x3i - x1r;
                    x0i = x3r + x1i;
                    a[j3][i] = wk1r * (x0r + x0i);
                    a[j3][i + 1] = wk1r * (x0r - x0i);
                }
            }
            k1 = 1;
            ks = -1;
            for (k = (m << 1); k <= n - m; k += m) {
                k1++;
                ks = -ks;
                wk1r = w[k1 << 1];
                wk1i = w[(k1 << 1) + 1];
                wk2r = ks * w[k1];
                wk2i = w[k1 + ks];
                wk3r = wk1r - 2 * wk2i * wk1i;
                wk3i = 2 * wk2i * wk1r - wk1i;
                for (j = k; j <= l + k - 1; j++) {
                    j1 = j + l;
                    j2 = j1 + l;
                    j3 = j2 + l;
                    for (i = 0; i <= n2 - 2; i += 2) {
                        x0r = a[j][i] + a[j1][i];
                        x0i = a[j][i + 1] + a[j1][i + 1];
                        x1r = a[j][i] - a[j1][i];
                        x1i = a[j][i + 1] - a[j1][i + 1];
                        x2r = a[j2][i] + a[j3][i];
                        x2i = a[j2][i + 1] + a[j3][i + 1];
                        x3r = a[j2][i] - a[j3][i];
                        x3i = a[j2][i + 1] - a[j3][i + 1];
                        a[j][i] = x0r + x2r;
                        a[j][i + 1] = x0i + x2i;
                        x0r -= x2r;
                        x0i -= x2i;
                        a[j2][i] = wk2r * x0r + wk2i * x0i;
                        a[j2][i + 1] = wk2r * x0i - wk2i * x0r;
                        x0r = x1r + x3i;
                        x0i = x1i - x3r;
                        a[j1][i] = wk1r * x0r + wk1i * x0i;
                        a[j1][i + 1] = wk1r * x0i - wk1i * x0r;
                        x0r = x1r - x3i;
                        x0i = x1i + x3r;
                        a[j3][i] = wk3r * x0r + wk3i * x0i;
                        a[j3][i + 1] = wk3r * x0i - wk3i * x0r;
                    }
                }
            }
        }
        l = m;
    }
    if (l < n) {
        for (j = 0; j <= l - 1; j++) {
            j1 = j + l;
            for (i = 0; i <= n2 - 2; i += 2) {
                x0r = a[j][i] - a[j1][i];
                x0i = a[j][i + 1] - a[j1][i + 1];
                a[j][i] += a[j1][i];
                a[j][i + 1] += a[j1][i + 1];
                a[j1][i] = x0r;
                a[j1][i + 1] = x0i;
            }
        }
    }
}


void rftbcol(int n1, int n, double **a, int nc, double *c)
{
    int i, j, k, kk, ks;
    double wkr, wki, xr, xi, yr, yi;
    
    ks = (nc << 2) / n;
    for (i = 0; i <= n1 - 1; i++) {
        kk = 0;
        for (k = (n >> 1) - 2; k >= 2; k -= 2) {
            j = n - k;
            kk += ks;
            wkr = 0.5 - c[kk];
            wki = c[nc - kk];
            xr = a[i][k] - a[i][j];
            xi = a[i][k + 1] + a[i][j + 1];
            yr = wkr * xr - wki * xi;
            yi = wkr * xi + wki * xr;
            a[i][k] -= yr;
            a[i][k + 1] -= yi;
            a[i][j] += yr;
            a[i][j + 1] -= yi;
        }
    }
}


void rftfcol(int n1, int n, double **a, int nc, double *c)
{
    int i, j, k, kk, ks;
    double wkr, wki, xr, xi, yr, yi;
    
    ks = (nc << 2) / n;
    for (i = 0; i <= n1 - 1; i++) {
        kk = 0;
        for (k = (n >> 1) - 2; k >= 2; k -= 2) {
            j = n - k;
            kk += ks;
            wkr = 0.5 - c[kk];
            wki = c[nc - kk];
            xr = a[i][k] - a[i][j];
            xi = a[i][k + 1] + a[i][j + 1];
            yr = wkr * xr + wki * xi;
            yi = wkr * xi - wki * xr;
            a[i][k] -= yr;
            a[i][k + 1] -= yi;
            a[i][j] += yr;
            a[i][j + 1] -= yi;
        }
    }
}


void dctbsub(int n1, int n2, double **a, int nc, double *c)
{
    int kk1, kk2, ks1, ks2, n1h, j1, k1, k2;
    double w1r, w1i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i;
    
    ks1 = nc / n1;
    ks2 = nc / n2;
    n1h = n1 >> 1;
    kk1 = ks1;
    for (k1 = 1; k1 <= n1h - 1; k1++) {
        j1 = n1 - k1;
        w1r = 2 * c[kk1];
        w1i = 2 * c[nc - kk1];
        kk1 += ks1;
        kk2 = ks2;
        for (k2 = 2; k2 <= n2 - 2; k2 += 2) {
            x0r = w1r * c[kk2];
            x0i = w1i * c[kk2];
            x1r = w1r * c[nc - kk2];
            x1i = w1i * c[nc - kk2];
            wkr = x0r - x1i;
            wki = x0i + x1r;
            wji = x0r + x1i;
            wjr = x0i - x1r;
            kk2 += ks2;
            x0r = wkr * a[k1][k2] - wki * a[k1][k2 + 1];
            x0i = wkr * a[k1][k2 + 1] + wki * a[k1][k2];
            x1r = wjr * a[j1][k2] - wji * a[j1][k2 + 1];
            x1i = wjr * a[j1][k2 + 1] + wji * a[j1][k2];
            a[k1][k2] = x0r + x1i;
            a[k1][k2 + 1] = x0i - x1r;
            a[j1][k2] = x1r + x0i;
            a[j1][k2 + 1] = x1i - x0r;
        }
        wkr = w1r * 0.5;
        wki = w1i * 0.5;
        wjr = w1r * c[kk2];
        wji = w1i * c[kk2];
        x0r = a[k1][0] + a[j1][0];
        x0i = a[k1][1] - a[j1][1];
        x1r = a[k1][0] - a[j1][0];
        x1i = a[k1][1] + a[j1][1];
        a[k1][0] = wkr * x0r - wki * x0i;
        a[k1][1] = wkr * x0i + wki * x0r;
        a[j1][0] = -wjr * x1r + wji * x1i;
        a[j1][1] = wjr * x1i + wji * x1r;
    }
    w1r = 2 * c[kk1];
    kk2 = ks2;
    for (k2 = 2; k2 <= n2 - 2; k2 += 2) {
        wkr = 2 * c[kk2];
        wki = 2 * c[nc - kk2];
        wjr = w1r * wkr;
        wji = w1r * wki;
        kk2 += ks2;
        x0i = wkr * a[0][k2 + 1] + wki * a[0][k2];
        a[0][k2] = wkr * a[0][k2] - wki * a[0][k2 + 1];
        a[0][k2 + 1] = x0i;
        x0i = wjr * a[n1h][k2 + 1] + wji * a[n1h][k2];
        a[n1h][k2] = wjr * a[n1h][k2] - wji * a[n1h][k2 + 1];
        a[n1h][k2 + 1] = x0i;
    }
    a[0][1] *= w1r;
    a[n1h][0] *= w1r;
    a[n1h][1] *= 0.5;
}


void dctfsub(int n1, int n2, double **a, int nc, double *c)
{
    int kk1, kk2, ks1, ks2, n1h, j1, k1, k2;
    double w1r, w1i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i;
    
    ks1 = nc / n1;
    ks2 = nc / n2;
    n1h = n1 >> 1;
    kk1 = ks1;
    for (k1 = 1; k1 <= n1h - 1; k1++) {
        j1 = n1 - k1;
        w1r = 2 * c[kk1];
        w1i = 2 * c[nc - kk1];
        kk1 += ks1;
        kk2 = ks2;
        for (k2 = 2; k2 <= n2 - 2; k2 += 2) {
            x0r = w1r * c[kk2];
            x0i = w1i * c[kk2];
            x1r = w1r * c[nc - kk2];
            x1i = w1i * c[nc - kk2];
            wkr = x0r - x1i;
            wki = x0i + x1r;
            wji = x0r + x1i;
            wjr = x0i - x1r;
            kk2 += ks2;
            x0r = a[k1][k2] - a[j1][k2 + 1];
            x0i = a[j1][k2] + a[k1][k2 + 1];
            x1r = a[j1][k2] - a[k1][k2 + 1];
            x1i = a[k1][k2] + a[j1][k2 + 1];
            a[k1][k2] = wkr * x0r + wki * x0i;
            a[k1][k2 + 1] = wkr * x0i - wki * x0r;
            a[j1][k2] = wjr * x1r + wji * x1i;
            a[j1][k2 + 1] = wjr * x1i - wji * x1r;
        }
        x0r = 2 * c[kk2];
        wjr = x0r * w1r;
        wji = x0r * w1i;
        x0r = w1r * a[k1][0] + w1i * a[k1][1];
        x0i = w1r * a[k1][1] - w1i * a[k1][0];
        x1r = -wjr * a[j1][0] + wji * a[j1][1];
        x1i = wjr * a[j1][1] + wji * a[j1][0];
        a[k1][0] = x0r + x1r;
        a[k1][1] = x1i + x0i;
        a[j1][0] = x0r - x1r;
        a[j1][1] = x1i - x0i;
    }
    w1r = 2 * c[kk1];
    kk2 = ks2;
    for (k2 = 2; k2 <= n2 - 2; k2 += 2) {
        wkr = 2 * c[kk2];
        wki = 2 * c[nc - kk2];
        wjr = w1r * wkr;
        wji = w1r * wki;
        kk2 += ks2;
        x0i = wkr * a[0][k2 + 1] - wki * a[0][k2];
        a[0][k2] = wkr * a[0][k2] + wki * a[0][k2 + 1];
        a[0][k2 + 1] = x0i;
        x0i = wjr * a[n1h][k2 + 1] - wji * a[n1h][k2];
        a[n1h][k2] = wjr * a[n1h][k2] + wji * a[n1h][k2 + 1];
        a[n1h][k2 + 1] = x0i;
    }
    w1r *= 2;
    a[0][0] *= 2;
    a[0][1] *= w1r;
    a[n1h][0] *= w1r;
}


void dstbsub(int n1, int n2, double **a, int nc, double *c)
{
    int kk1, kk2, ks1, ks2, n1h, j1, k1, k2;
    double w1r, w1i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i;
    
    ks1 = nc / n1;
    ks2 = nc / n2;
    n1h = n1 >> 1;
    kk1 = ks1;
    for (k1 = 1; k1 <= n1h - 1; k1++) {
        j1 = n1 - k1;
        w1r = 2 * c[kk1];
        w1i = 2 * c[nc - kk1];
        kk1 += ks1;
        kk2 = ks2;
        for (k2 = 2; k2 <= n2 - 2; k2 += 2) {
            x0r = w1r * c[kk2];
            x0i = w1i * c[kk2];
            x1r = w1r * c[nc - kk2];
            x1i = w1i * c[nc - kk2];
            wkr = x0r - x1i;
            wki = x0i + x1r;
            wji = x0r + x1i;
            wjr = x0i - x1r;
            kk2 += ks2;
            x0r = wkr * a[k1][k2] - wki * a[k1][k2 + 1];
            x0i = wkr * a[k1][k2 + 1] + wki * a[k1][k2];
            x1r = wjr * a[j1][k2] - wji * a[j1][k2 + 1];
            x1i = wjr * a[j1][k2 + 1] + wji * a[j1][k2];
            a[k1][k2] = x1i - x0r;
            a[k1][k2 + 1] = x1r + x0i;
            a[j1][k2] = x0i - x1r;
            a[j1][k2 + 1] = x0r + x1i;
        }
        wkr = w1r * 0.5;
        wki = w1i * 0.5;
        wjr = w1r * c[kk2];
        wji = w1i * c[kk2];
        x0r = a[k1][0] + a[j1][0];
        x0i = a[k1][1] - a[j1][1];
        x1r = a[k1][0] - a[j1][0];
        x1i = a[k1][1] + a[j1][1];
        a[k1][1] = wkr * x0r - wki * x0i;
        a[k1][0] = wkr * x0i + wki * x0r;
        a[j1][1] = -wjr * x1r + wji * x1i;
        a[j1][0] = wjr * x1i + wji * x1r;
    }
    w1r = 2 * c[kk1];
    kk2 = ks2;
    for (k2 = 2; k2 <= n2 - 2; k2 += 2) {
        wkr = 2 * c[kk2];
        wki = 2 * c[nc - kk2];
        wjr = w1r * wkr;
        wji = w1r * wki;
        kk2 += ks2;
        x0i = wkr * a[0][k2 + 1] + wki * a[0][k2];
        a[0][k2 + 1] = wkr * a[0][k2] - wki * a[0][k2 + 1];
        a[0][k2] = x0i;
        x0i = wjr * a[n1h][k2 + 1] + wji * a[n1h][k2];
        a[n1h][k2 + 1] = wjr * a[n1h][k2] - wji * a[n1h][k2 + 1];
        a[n1h][k2] = x0i;
    }
    a[0][1] *= w1r;
    a[n1h][0] *= w1r;
    a[n1h][1] *= 0.5;
}


void dstfsub(int n1, int n2, double **a, int nc, double *c)
{
    int kk1, kk2, ks1, ks2, n1h, j1, k1, k2;
    double w1r, w1i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i;
    
    ks1 = nc / n1;
    ks2 = nc / n2;
    n1h = n1 >> 1;
    kk1 = ks1;
    for (k1 = 1; k1 <= n1h - 1; k1++) {
        j1 = n1 - k1;
        w1r = 2 * c[kk1];
        w1i = 2 * c[nc - kk1];
        kk1 += ks1;
        kk2 = ks2;
        for (k2 = 2; k2 <= n2 - 2; k2 += 2) {
            x0r = w1r * c[kk2];
            x0i = w1i * c[kk2];
            x1r = w1r * c[nc - kk2];
            x1i = w1i * c[nc - kk2];
            wkr = x0r - x1i;
            wki = x0i + x1r;
            wji = x0r + x1i;
            wjr = x0i - x1r;
            kk2 += ks2;
            x0r = a[j1][k2 + 1] - a[k1][k2];
            x0i = a[k1][k2 + 1] + a[j1][k2];
            x1r = a[k1][k2 + 1] - a[j1][k2];
            x1i = a[j1][k2 + 1] + a[k1][k2];
            a[k1][k2] = wkr * x0r + wki * x0i;
            a[k1][k2 + 1] = wkr * x0i - wki * x0r;
            a[j1][k2] = wjr * x1r + wji * x1i;
            a[j1][k2 + 1] = wjr * x1i - wji * x1r;
        }
        x0r = 2 * c[kk2];
        wjr = x0r * w1r;
        wji = x0r * w1i;
        x0r = w1r * a[k1][1] + w1i * a[k1][0];
        x0i = w1r * a[k1][0] - w1i * a[k1][1];
        x1r = -wjr * a[j1][1] + wji * a[j1][0];
        x1i = wjr * a[j1][0] + wji * a[j1][1];
        a[k1][0] = x0r + x1r;
        a[k1][1] = x1i + x0i;
        a[j1][0] = x0r - x1r;
        a[j1][1] = x1i - x0i;
    }
    w1r = 2 * c[kk1];
    kk2 = ks2;
    for (k2 = 2; k2 <= n2 - 2; k2 += 2) {
        wkr = 2 * c[kk2];
        wki = 2 * c[nc - kk2];
        wjr = w1r * wkr;
        wji = w1r * wki;
        kk2 += ks2;
        x0i = wkr * a[0][k2] - wki * a[0][k2 + 1];
        a[0][k2] = wkr * a[0][k2 + 1] + wki * a[0][k2];
        a[0][k2 + 1] = x0i;
        x0i = wjr * a[n1h][k2] - wji * a[n1h][k2 + 1];
        a[n1h][k2] = wjr * a[n1h][k2 + 1] + wji * a[n1h][k2];
        a[n1h][k2 + 1] = x0i;
    }
    w1r *= 2;
    a[0][0] *= 2;
    a[0][1] *= w1r;
    a[n1h][0] *= w1r;
}