fftsg3d.c

2维fft程序

                        a[j1][j2][j3] *= 8.0 / n1 / n2 / n3;
                    }
                }
            }
        .
*/


#include <stdio.h>
#include <stdlib.h>
#define fft3d_alloc_error_check(p) { \
    if ((p) == NULL) { \
        fprintf(stderr, "fft3d memory allocation error\n"); \
        exit(1); \
    } \
}


#ifdef USE_FFT3D_PTHREADS
#define USE_FFT3D_THREADS
#ifndef FFT3D_MAX_THREADS
#define FFT3D_MAX_THREADS 4
#endif
#ifndef FFT3D_THREADS_BEGIN_N
#define FFT3D_THREADS_BEGIN_N 65536
#endif
#include <pthread.h>
#define fft3d_thread_t pthread_t
#define fft3d_thread_create(thp,func,argp) { \
    if (pthread_create(thp, NULL, func, (void *) (argp)) != 0) { \
        fprintf(stderr, "fft3d thread error\n"); \
        exit(1); \
    } \
}
#define fft3d_thread_wait(th) { \
    if (pthread_join(th, NULL) != 0) { \
        fprintf(stderr, "fft3d thread error\n"); \
        exit(1); \
    } \
}
#endif /* USE_FFT3D_PTHREADS */


#ifdef USE_FFT3D_WINTHREADS
#define USE_FFT3D_THREADS
#ifndef FFT3D_MAX_THREADS
#define FFT3D_MAX_THREADS 4
#endif
#ifndef FFT3D_THREADS_BEGIN_N
#define FFT3D_THREADS_BEGIN_N 131072
#endif
#include <windows.h>
#define fft3d_thread_t HANDLE
#define fft3d_thread_create(thp,func,argp) { \
    DWORD thid; \
    *(thp) = CreateThread(NULL, 0, (LPTHREAD_START_ROUTINE) (func), (LPVOID) (argp), 0, &thid); \
    if (*(thp) == 0) { \
        fprintf(stderr, "fft3d thread error\n"); \
        exit(1); \
    } \
}
#define fft3d_thread_wait(th) { \
    WaitForSingleObject(th, INFINITE); \
    CloseHandle(th); \
}
#endif /* USE_FFT3D_WINTHREADS */


void cdft3d(int n1, int n2, int n3, int isgn, double ***a, 
    double *t, int *ip, double *w)
{
    void makewt(int nw, int *ip, double *w);
    void xdft3da_sub(int n1, int n2, int n3, int icr, int isgn, 
        double ***a, double *t, int *ip, double *w);
    void cdft3db_sub(int n1, int n2, int n3, int isgn, double ***a, 
        double *t, int *ip, double *w);
#ifdef USE_FFT3D_THREADS
    void xdft3da_subth(int n1, int n2, int n3, int icr, int isgn, 
        double ***a, double *t, int *ip, double *w);
    void cdft3db_subth(int n1, int n2, int n3, int isgn, double ***a, 
        double *t, int *ip, double *w);
#endif /* USE_FFT3D_THREADS */
    int n, itnull, nt;
    
    n = n1;
    if (n < n2) {
        n = n2;
    }
    n <<= 1;
    if (n < n3) {
        n = n3;
    }
    if (n > (ip[0] << 2)) {
        makewt(n >> 2, ip, w);
    }
    itnull = 0;
    if (t == NULL) {
        itnull = 1;
        nt = n1;
        if (nt < n2) {
            nt = n2;
        }
        nt *= 8;
#ifdef USE_FFT3D_THREADS
        nt *= FFT3D_MAX_THREADS;
#endif /* USE_FFT3D_THREADS */
        if (n3 == 4) {
            nt >>= 1;
        } else if (n3 < 4) {
            nt >>= 2;
        }
        t = (double *) malloc(sizeof(double) * nt);
        fft3d_alloc_error_check(t);
    }
#ifdef USE_FFT3D_THREADS
    if ((double) n1 * n2 * n3 >= (double) FFT3D_THREADS_BEGIN_N) {
        xdft3da_subth(n1, n2, n3, 0, isgn, a, t, ip, w);
        cdft3db_subth(n1, n2, n3, isgn, a, t, ip, w);
    } else 
#endif /* USE_FFT3D_THREADS */
    {
        xdft3da_sub(n1, n2, n3, 0, isgn, a, t, ip, w);
        cdft3db_sub(n1, n2, n3, isgn, a, t, ip, w);
    }
    if (itnull != 0) {
        free(t);
    }
}


void rdft3d(int n1, int n2, int n3, 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 xdft3da_sub(int n1, int n2, int n3, int icr, int isgn, 
        double ***a, double *t, int *ip, double *w);
    void cdft3db_sub(int n1, int n2, int n3, int isgn, double ***a, 
        double *t, int *ip, double *w);
    void rdft3d_sub(int n1, int n2, int n3, int isgn, double ***a);
#ifdef USE_FFT3D_THREADS
    void xdft3da_subth(int n1, int n2, int n3, int icr, int isgn, 
        double ***a, double *t, int *ip, double *w);
    void cdft3db_subth(int n1, int n2, int n3, int isgn, double ***a, 
        double *t, int *ip, double *w);
#endif /* USE_FFT3D_THREADS */
    int n, nw, nc, itnull, nt;
    
    n = n1;
    if (n < n2) {
        n = n2;
    }
    n <<= 1;
    if (n < n3) {
        n = n3;
    }
    nw = ip[0];
    if (n > (nw << 2)) {
        nw = n >> 2;
        makewt(nw, ip, w);
    }
    nc = ip[1];
    if (n3 > (nc << 2)) {
        nc = n3 >> 2;
        makect(nc, ip, w + nw);
    }
    itnull = 0;
    if (t == NULL) {
        itnull = 1;
        nt = n1;
        if (nt < n2) {
            nt = n2;
        }
        nt *= 8;
#ifdef USE_FFT3D_THREADS
        nt *= FFT3D_MAX_THREADS;
#endif /* USE_FFT3D_THREADS */
        if (n3 == 4) {
            nt >>= 1;
        } else if (n3 < 4) {
            nt >>= 2;
        }
        t = (double *) malloc(sizeof(double) * nt);
        fft3d_alloc_error_check(t);
    }
#ifdef USE_FFT3D_THREADS
    if ((double) n1 * n2 * n3 >= (double) FFT3D_THREADS_BEGIN_N) {
        if (isgn < 0) {
            rdft3d_sub(n1, n2, n3, isgn, a);
            cdft3db_subth(n1, n2, n3, isgn, a, t, ip, w);
        }
        xdft3da_subth(n1, n2, n3, 1, isgn, a, t, ip, w);
        if (isgn >= 0) {
            cdft3db_subth(n1, n2, n3, isgn, a, t, ip, w);
            rdft3d_sub(n1, n2, n3, isgn, a);
        }
    } else 
#endif /* USE_FFT3D_THREADS */
    {
        if (isgn < 0) {
            rdft3d_sub(n1, n2, n3, isgn, a);
            cdft3db_sub(n1, n2, n3, isgn, a, t, ip, w);
        }
        xdft3da_sub(n1, n2, n3, 1, isgn, a, t, ip, w);
        if (isgn >= 0) {
            cdft3db_sub(n1, n2, n3, isgn, a, t, ip, w);
            rdft3d_sub(n1, n2, n3, isgn, a);
        }
    }
    if (itnull != 0) {
        free(t);
    }
}


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


void ddct3d(int n1, int n2, int n3, 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 ddxt3da_sub(int n1, int n2, int n3, int ics, int isgn, 
        double ***a, double *t, int *ip, double *w);
    void ddxt3db_sub(int n1, int n2, int n3, int ics, int isgn, 
        double ***a, double *t, int *ip, double *w);
#ifdef USE_FFT3D_THREADS
    void ddxt3da_subth(int n1, int n2, int n3, int ics, int isgn, 
        double ***a, double *t, int *ip, double *w);
    void ddxt3db_subth(int n1, int n2, int n3, int ics, int isgn, 
        double ***a, double *t, int *ip, double *w);
#endif /* USE_FFT3D_THREADS */
    int n, nw, nc, itnull, nt;
    
    n = n1;
    if (n < n2) {
        n = n2;
    }
    if (n < n3) {
        n = n3;
    }
    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;
        nt = n1;
        if (nt < n2) {
            nt = n2;
        }
        nt *= 4;
#ifdef USE_FFT3D_THREADS
        nt *= FFT3D_MAX_THREADS;
#endif /* USE_FFT3D_THREADS */
        if (n3 == 2) {
            nt >>= 1;
        }
        t = (double *) malloc(sizeof(double) * nt);
        fft3d_alloc_error_check(t);
    }
#ifdef USE_FFT3D_THREADS
    if ((double) n1 * n2 * n3 >= (double) FFT3D_THREADS_BEGIN_N) {
        ddxt3da_subth(n1, n2, n3, 0, isgn, a, t, ip, w);
        ddxt3db_subth(n1, n2, n3, 0, isgn, a, t, ip, w);
    } else 
#endif /* USE_FFT3D_THREADS */
    {
        ddxt3da_sub(n1, n2, n3, 0, isgn, a, t, ip, w);
        ddxt3db_sub(n1, n2, n3, 0, isgn, a, t, ip, w);
    }
    if (itnull != 0) {
        free(t);
    }
}


void ddst3d(int n1, int n2, int n3, 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 ddxt3da_sub(int n1, int n2, int n3, int ics, int isgn, 
        double ***a, double *t, int *ip, double *w);
    void ddxt3db_sub(int n1, int n2, int n3, int ics, int isgn, 
        double ***a, double *t, int *ip, double *w);
#ifdef USE_FFT3D_THREADS
    void ddxt3da_subth(int n1, int n2, int n3, int ics, int isgn, 
        double ***a, double *t, int *ip, double *w);
    void ddxt3db_subth(int n1, int n2, int n3, int ics, int isgn, 
        double ***a, double *t, int *ip, double *w);
#endif /* USE_FFT3D_THREADS */
    int n, nw, nc, itnull, nt;
    
    n = n1;
    if (n < n2) {
        n = n2;
    }
    if (n < n3) {
        n = n3;
    }
    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;
        nt = n1;
        if (nt < n2) {
            nt = n2;
        }
        nt *= 4;
#ifdef USE_FFT3D_THREADS
        nt *= FFT3D_MAX_THREADS;
#endif /* USE_FFT3D_THREADS */
        if (n3 == 2) {
            nt >>= 1;
        }
        t = (double *) malloc(sizeof(double) * nt);
        fft3d_alloc_error_check(t);
    }
#ifdef USE_FFT3D_THREADS
    if ((double) n1 * n2 * n3 >= (double) FFT3D_THREADS_BEGIN_N) {
        ddxt3da_subth(n1, n2, n3, 1, isgn, a, t, ip, w);
        ddxt3db_subth(n1, n2, n3, 1, isgn, a, t, ip, w);
    } else 
#endif /* USE_FFT3D_THREADS */
    {
        ddxt3da_sub(n1, n2, n3, 1, isgn, a, t, ip, w);
        ddxt3db_sub(n1, n2, n3, 1, isgn, a, t, ip, w);
    }