fftsg2d.c

2维fft程序

    int *ip, double *w)
{
    void makewt(int nw, int *ip, double *w);
    void cdft(int n, int isgn, double *a, int *ip, double *w);
    void cdft2d_sub(int n1, int n2, int isgn, double **a, double *t, 
        int *ip, double *w);
#ifdef USE_FFT2D_THREADS
    void xdft2d0_subth(int n1, int n2, int icr, int isgn, double **a, 
        int *ip, double *w);
    void cdft2d_subth(int n1, int n2, int isgn, double **a, double *t, 
        int *ip, double *w);
#endif /* USE_FFT2D_THREADS */
    int n, itnull, nthread, nt, i;
    
    n = n1 << 1;
    if (n < n2) {
        n = n2;
    }
    if (n > (ip[0] << 2)) {
        makewt(n >> 2, ip, w);
    }
    itnull = 0;
    if (t == NULL) {
        itnull = 1;
        nthread = 1;
#ifdef USE_FFT2D_THREADS
        nthread = FFT2D_MAX_THREADS;
#endif /* USE_FFT2D_THREADS */
        nt = 8 * nthread * n1;
        if (n2 == 4 * nthread) {
            nt >>= 1;
        } else if (n2 < 4 * nthread) {
            nt >>= 2;
        }
        t = (double *) malloc(sizeof(double) * nt);
        fft2d_alloc_error_check(t);
    }
#ifdef USE_FFT2D_THREADS
    if ((double) n1 * n2 >= (double) FFT2D_THREADS_BEGIN_N) {
        xdft2d0_subth(n1, n2, 0, isgn, a, ip, w);
        cdft2d_subth(n1, n2, isgn, a, t, ip, w);
    } else 
#endif /* USE_FFT2D_THREADS */
    {
        for (i = 0; i < n1; i++) {
            cdft(n2, isgn, a[i], ip, w);
        }
        cdft2d_sub(n1, n2, isgn, a, t, ip, w);
    }
    if (itnull != 0) {
        free(t);
    }
}


void rdft2d(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 rdft(int n, int isgn, double *a, int *ip, double *w);
    void cdft2d_sub(int n1, int n2, int isgn, double **a, double *t, 
        int *ip, double *w);
    void rdft2d_sub(int n1, int n2, int isgn, double **a);
#ifdef USE_FFT2D_THREADS
    void xdft2d0_subth(int n1, int n2, int icr, int isgn, double **a, 
        int *ip, double *w);
    void cdft2d_subth(int n1, int n2, int isgn, double **a, double *t, 
        int *ip, double *w);
#endif /* USE_FFT2D_THREADS */
    int n, nw, nc, itnull, nthread, nt, i;
    
    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 (n2 > (nc << 2)) {
        nc = n2 >> 2;
        makect(nc, ip, w + nw);
    }
    itnull = 0;
    if (t == NULL) {
        itnull = 1;
        nthread = 1;
#ifdef USE_FFT2D_THREADS
        nthread = FFT2D_MAX_THREADS;
#endif /* USE_FFT2D_THREADS */
        nt = 8 * nthread * n1;
        if (n2 == 4 * nthread) {
            nt >>= 1;
        } else if (n2 < 4 * nthread) {
            nt >>= 2;
        }
        t = (double *) malloc(sizeof(double) * nt);
        fft2d_alloc_error_check(t);
    }
#ifdef USE_FFT2D_THREADS
    if ((double) n1 * n2 >= (double) FFT2D_THREADS_BEGIN_N) {
        if (isgn < 0) {
            rdft2d_sub(n1, n2, isgn, a);
            cdft2d_subth(n1, n2, isgn, a, t, ip, w);
        }
        xdft2d0_subth(n1, n2, 1, isgn, a, ip, w);
        if (isgn >= 0) {
            cdft2d_subth(n1, n2, isgn, a, t, ip, w);
            rdft2d_sub(n1, n2, isgn, a);
        }
    } else 
#endif /* USE_FFT2D_THREADS */
    {
        if (isgn < 0) {
            rdft2d_sub(n1, n2, isgn, a);
            cdft2d_sub(n1, n2, isgn, a, t, ip, w);
        }
        for (i = 0; i < n1; i++) {
            rdft(n2, isgn, a[i], ip, w);
        }
        if (isgn >= 0) {
            cdft2d_sub(n1, n2, isgn, a, t, ip, w);
            rdft2d_sub(n1, n2, isgn, a);
        }
    }
    if (itnull != 0) {
        free(t);
    }
}


void rdft2dsort(int n1, int n2, int isgn, double **a)
{
    int n1h, i;
    double x, y;
    
    n1h = n1 >> 1;
    if (isgn < 0) {
        for (i = n1h + 1; i < n1; i++) {
            a[i][0] = a[i][n2 + 1];
            a[i][1] = a[i][n2];
        }
        a[0][1] = a[0][n2];
        a[n1h][1] = a[n1h][n2];
    } else {
        for (i = n1h + 1; i < n1; i++) {
            y = a[i][0];
            x = a[i][1];
            a[i][n2] = x;
            a[i][n2 + 1] = y;
            a[n1 - i][n2] = x;
            a[n1 - i][n2 + 1] = -y;
            a[i][0] = a[n1 - i][0];
            a[i][1] = -a[n1 - i][1];
        }
        a[0][n2] = a[0][1];
        a[0][n2 + 1] = 0;
        a[0][1] = 0;
        a[n1h][n2] = a[n1h][1];
        a[n1h][n2 + 1] = 0;
        a[n1h][1] = 0;
    }
}


void ddct2d(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 ddct(int n, int isgn, double *a, int *ip, double *w);
    void ddxt2d_sub(int n1, int n2, int ics, int isgn, double **a, 
        double *t, int *ip, double *w);
#ifdef USE_FFT2D_THREADS
    void ddxt2d0_subth(int n1, int n2, int ics, int isgn, double **a, 
        int *ip, double *w);
    void ddxt2d_subth(int n1, int n2, int ics, int isgn, double **a, 
        double *t, int *ip, double *w);
#endif /* USE_FFT2D_THREADS */
    int n, nw, nc, itnull, nthread, nt, i;
    
    n = n1;
    if (n < n2) {
        n = n2;
    }
    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);
    }
    itnull = 0;
    if (t == NULL) {
        itnull = 1;
        nthread = 1;
#ifdef USE_FFT2D_THREADS
        nthread = FFT2D_MAX_THREADS;
#endif /* USE_FFT2D_THREADS */
        nt = 4 * nthread * n1;
        if (n2 == 2 * nthread) {
            nt >>= 1;
        } else if (n2 < 2 * nthread) {
            nt >>= 2;
        }
        t = (double *) malloc(sizeof(double) * nt);
        fft2d_alloc_error_check(t);
    }
#ifdef USE_FFT2D_THREADS
    if ((double) n1 * n2 >= (double) FFT2D_THREADS_BEGIN_N) {
        ddxt2d0_subth(n1, n2, 0, isgn, a, ip, w);
        ddxt2d_subth(n1, n2, 0, isgn, a, t, ip, w);
    } else 
#endif /* USE_FFT2D_THREADS */
    {
        for (i = 0; i < n1; i++) {
            ddct(n2, isgn, a[i], ip, w);
        }
        ddxt2d_sub(n1, n2, 0, isgn, a, t, ip, w);
    }
    if (itnull != 0) {
        free(t);
    }
}


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 ddst(int n, int isgn, double *a, int *ip, double *w);
    void ddxt2d_sub(int n1, int n2, int ics, int isgn, double **a, 
        double *t, int *ip, double *w);
#ifdef USE_FFT2D_THREADS
    void ddxt2d0_subth(int n1, int n2, int ics, int isgn, double **a, 
        int *ip, double *w);
    void ddxt2d_subth(int n1, int n2, int ics, int isgn, double **a, 
        double *t, int *ip, double *w);
#endif /* USE_FFT2D_THREADS */
    int n, nw, nc, itnull, nthread, nt, i;
    
    n = n1;
    if (n < n2) {
        n = n2;
    }
    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);
    }
    itnull = 0;
    if (t == NULL) {
        itnull = 1;
        nthread = 1;
#ifdef USE_FFT2D_THREADS
        nthread = FFT2D_MAX_THREADS;
#endif /* USE_FFT2D_THREADS */
        nt = 4 * nthread * n1;
        if (n2 == 2 * nthread) {
            nt >>= 1;
        } else if (n2 < 2 * nthread) {
            nt >>= 2;
        }
        t = (double *) malloc(sizeof(double) * nt);
        fft2d_alloc_error_check(t);
    }
#ifdef USE_FFT2D_THREADS
    if ((double) n1 * n2 >= (double) FFT2D_THREADS_BEGIN_N) {
        ddxt2d0_subth(n1, n2, 1, isgn, a, ip, w);
        ddxt2d_subth(n1, n2, 1, isgn, a, t, ip, w);
    } else 
#endif /* USE_FFT2D_THREADS */
    {
        for (i = 0; i < n1; i++) {
            ddst(n2, isgn, a[i], ip, w);
        }
        ddxt2d_sub(n1, n2, 1, isgn, a, t, ip, w);
    }
    if (itnull != 0) {
        free(t);
    }
}


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


void cdft2d_sub(int n1, int n2, int isgn, double **a, double *t, 
    int *ip, double *w)
{
    void cdft(int n, int isgn, double *a, int *ip, double *w);
    int i, j;
    
    if (n2 > 4) {
        for (j = 0; j < n2; j += 8) {
            for (i = 0; i < n1; i++) {
                t[2 * i] = a[i][j];
                t[2 * i + 1] = a[i][j + 1];
                t[2 * n1 + 2 * i] = a[i][j + 2];
                t[2 * n1 + 2 * i + 1] = a[i][j + 3];
                t[4 * n1 + 2 * i] = a[i][j + 4];
                t[4 * n1 + 2 * i + 1] = a[i][j + 5];
                t[6 * n1 + 2 * i] = a[i][j + 6];
                t[6 * n1 + 2 * i + 1] = a[i][j + 7];
            }
            cdft(2 * n1, isgn, t, ip, w);
            cdft(2 * n1, isgn, &t[2 * n1], ip, w);
            cdft(2 * n1, isgn, &t[4 * n1], ip, w);
            cdft(2 * n1, isgn, &t[6 * n1], ip, w);
            for (i = 0; i < n1; i++) {
                a[i][j] = t[2 * i];
                a[i][j + 1] = t[2 * i + 1];
                a[i][j + 2] = t[2 * n1 + 2 * i];
                a[i][j + 3] = t[2 * n1 + 2 * i + 1];
                a[i][j + 4] = t[4 * n1 + 2 * i];
                a[i][j + 5] = t[4 * n1 + 2 * i + 1];
                a[i][j + 6] = t[6 * n1 + 2 * i];
                a[i][j + 7] = t[6 * n1 + 2 * i + 1];
            }
        }
    } else if (n2 == 4) {
        for (i = 0; i < n1; i++) {
            t[2 * i] = a[i][0];
            t[2 * i + 1] = a[i][1];
            t[2 * n1 + 2 * i] = a[i][2];
            t[2 * n1 + 2 * i + 1] = a[i][3];
        }
        cdft(2 * n1, isgn, t, ip, w);
        cdft(2 * n1, isgn, &t[2 * n1], ip, w);
        for (i = 0; i < n1; i++) {
            a[i][0] = t[2 * i];
            a[i][1] = t[2 * i + 1];
            a[i][2] = t[2 * n1 + 2 * i];
            a[i][3] = t[2 * n1 + 2 * i + 1];
        }
    } else if (n2 == 2) {
        for (i = 0; i < n1; i++) {
            t[2 * i] = a[i][0];
            t[2 * i + 1] = a[i][1];
        }
        cdft(2 * n1, isgn, t, ip, w);
        for (i = 0; i < n1; i++) {
            a[i][0] = t[2 * i];
            a[i][1] = t[2 * i + 1];
        }
    }
}


void rdft2d_sub(int n1, int n2, int isgn, double **a)
{
    int n1h, i, j;
    double xi;
    
    n1h = n1 >> 1;
    if (isgn < 0) {
        for (i = 1; i < n1h; i++) {
            j = n1 - i;
            xi = a[i][0] - a[j][0];
            a[i][0] += a[j][0];
            a[j][0] = xi;
            xi = a[j][1] - a[i][1];
            a[i][1] += a[j][1];
            a[j][1] = xi;
        }
    } else {
        for (i = 1; i < n1h; i++) {
            j = n1 - i;
            a[j][0] = 0.5 * (a[i][0] - a[j][0]);
            a[i][0] -= a[j][0];
            a[j][1] = 0.5 * (a[i][1] + a[j][1]);
            a[i][1] -= a[j][1];
        }
    }
}


void ddxt2d_sub(int n1, int n2, int ics, int isgn, double **a, 
    double *t, int *ip, double *w)
{
    void ddct(int n, int isgn, double *a, int *ip, double *w);
    void ddst(int n, int isgn, double *a, int *ip, double *w);
    int i, j;
    
    if (n2 > 2) {