fftsg2d.c

2维fft程序

        for (j = 0; j < n2; j += 4) {
            for (i = 0; i < n1; i++) {
                t[i] = a[i][j];
                t[n1 + i] = a[i][j + 1];
                t[2 * n1 + i] = a[i][j + 2];
                t[3 * n1 + i] = a[i][j + 3];
            }
            if (ics == 0) {
                ddct(n1, isgn, t, ip, w);
                ddct(n1, isgn, &t[n1], ip, w);
                ddct(n1, isgn, &t[2 * n1], ip, w);
                ddct(n1, isgn, &t[3 * n1], ip, w);
            } else {
                ddst(n1, isgn, t, ip, w);
                ddst(n1, isgn, &t[n1], ip, w);
                ddst(n1, isgn, &t[2 * n1], ip, w);
                ddst(n1, isgn, &t[3 * n1], ip, w);
            }
            for (i = 0; i < n1; i++) {
                a[i][j] = t[i];
                a[i][j + 1] = t[n1 + i];
                a[i][j + 2] = t[2 * n1 + i];
                a[i][j + 3] = t[3 * n1 + i];
            }
        }
    } else if (n2 == 2) {
        for (i = 0; i < n1; i++) {
            t[i] = a[i][0];
            t[n1 + i] = a[i][1];
        }
        if (ics == 0) {
            ddct(n1, isgn, t, ip, w);
            ddct(n1, isgn, &t[n1], ip, w);
        } else {
            ddst(n1, isgn, t, ip, w);
            ddst(n1, isgn, &t[n1], ip, w);
        }
        for (i = 0; i < n1; i++) {
            a[i][0] = t[i];
            a[i][1] = t[n1 + i];
        }
    }
}


#ifdef USE_FFT2D_THREADS
struct fft2d_arg_st {
    int nthread;
    int n0;
    int n1;
    int n2;
    int ic;
    int isgn;
    double **a;
    double *t;
    int *ip;
    double *w;
};
typedef struct fft2d_arg_st fft2d_arg_t;


void xdft2d0_subth(int n1, int n2, int icr, int isgn, double **a, 
    int *ip, double *w)
{
    void *xdft2d0_th(void *p);
    fft2d_thread_t th[FFT2D_MAX_THREADS];
    fft2d_arg_t ag[FFT2D_MAX_THREADS];
    int nthread, i;
    
    nthread = FFT2D_MAX_THREADS;
    if (nthread > n1) {
        nthread = n1;
    }
    for (i = 0; i < nthread; i++) {
        ag[i].nthread = nthread;
        ag[i].n0 = i;
        ag[i].n1 = n1;
        ag[i].n2 = n2;
        ag[i].ic = icr;
        ag[i].isgn = isgn;
        ag[i].a = a;
        ag[i].ip = ip;
        ag[i].w = w;
        fft2d_thread_create(&th[i], xdft2d0_th, &ag[i]);
    }
    for (i = 0; i < nthread; i++) {
        fft2d_thread_wait(th[i]);
    }
}


void cdft2d_subth(int n1, int n2, int isgn, double **a, double *t, 
    int *ip, double *w)
{
    void *cdft2d_th(void *p);
    fft2d_thread_t th[FFT2D_MAX_THREADS];
    fft2d_arg_t ag[FFT2D_MAX_THREADS];
    int nthread, nt, i;
    
    nthread = FFT2D_MAX_THREADS;
    nt = 8 * n1;
    if (n2 == 4 * FFT2D_MAX_THREADS) {
        nt >>= 1;
    } else if (n2 < 4 * FFT2D_MAX_THREADS) {
        nthread = n2 >> 1;
        nt >>= 2;
    }
    for (i = 0; i < nthread; i++) {
        ag[i].nthread = nthread;
        ag[i].n0 = i;
        ag[i].n1 = n1;
        ag[i].n2 = n2;
        ag[i].isgn = isgn;
        ag[i].a = a;
        ag[i].t = &t[nt * i];
        ag[i].ip = ip;
        ag[i].w = w;
        fft2d_thread_create(&th[i], cdft2d_th, &ag[i]);
    }
    for (i = 0; i < nthread; i++) {
        fft2d_thread_wait(th[i]);
    }
}


void ddxt2d0_subth(int n1, int n2, int ics, int isgn, double **a, 
    int *ip, double *w)
{
    void *ddxt2d0_th(void *p);
    fft2d_thread_t th[FFT2D_MAX_THREADS];
    fft2d_arg_t ag[FFT2D_MAX_THREADS];
    int nthread, i;
    
    nthread = FFT2D_MAX_THREADS;
    if (nthread > n1) {
        nthread = n1;
    }
    for (i = 0; i < nthread; i++) {
        ag[i].nthread = nthread;
        ag[i].n0 = i;
        ag[i].n1 = n1;
        ag[i].n2 = n2;
        ag[i].ic = ics;
        ag[i].isgn = isgn;
        ag[i].a = a;
        ag[i].ip = ip;
        ag[i].w = w;
        fft2d_thread_create(&th[i], ddxt2d0_th, &ag[i]);
    }
    for (i = 0; i < nthread; i++) {
        fft2d_thread_wait(th[i]);
    }
}


void ddxt2d_subth(int n1, int n2, int ics, int isgn, double **a, 
    double *t, int *ip, double *w)
{
    void *ddxt2d_th(void *p);
    fft2d_thread_t th[FFT2D_MAX_THREADS];
    fft2d_arg_t ag[FFT2D_MAX_THREADS];
    int nthread, nt, i;
    
    nthread = FFT2D_MAX_THREADS;
    nt = 4 * n1;
    if (n2 == 2 * FFT2D_MAX_THREADS) {
        nt >>= 1;
    } else if (n2 < 2 * FFT2D_MAX_THREADS) {
        nthread = n2;
        nt >>= 2;
    }
    for (i = 0; i < nthread; i++) {
        ag[i].nthread = nthread;
        ag[i].n0 = i;
        ag[i].n1 = n1;
        ag[i].n2 = n2;
        ag[i].ic = ics;
        ag[i].isgn = isgn;
        ag[i].a = a;
        ag[i].t = &t[nt * i];
        ag[i].ip = ip;
        ag[i].w = w;
        fft2d_thread_create(&th[i], ddxt2d_th, &ag[i]);
    }
    for (i = 0; i < nthread; i++) {
        fft2d_thread_wait(th[i]);
    }
}


void *xdft2d0_th(void *p)
{
    void cdft(int n, int isgn, double *a, int *ip, double *w);
    void rdft(int n, int isgn, double *a, int *ip, double *w);
    int nthread, n0, n1, n2, icr, isgn, *ip, i;
    double **a, *w;
    
    nthread = ((fft2d_arg_t *) p)->nthread;
    n0 = ((fft2d_arg_t *) p)->n0;
    n1 = ((fft2d_arg_t *) p)->n1;
    n2 = ((fft2d_arg_t *) p)->n2;
    icr = ((fft2d_arg_t *) p)->ic;
    isgn = ((fft2d_arg_t *) p)->isgn;
    a = ((fft2d_arg_t *) p)->a;
    ip = ((fft2d_arg_t *) p)->ip;
    w = ((fft2d_arg_t *) p)->w;
    if (icr == 0) {
        for (i = n0; i < n1; i += nthread) {
            cdft(n2, isgn, a[i], ip, w);
        }
    } else {
        for (i = n0; i < n1; i += nthread) {
            rdft(n2, isgn, a[i], ip, w);
        }
    }
    return (void *) 0;
}


void *cdft2d_th(void *p)
{
    void cdft(int n, int isgn, double *a, int *ip, double *w);
    int nthread, n0, n1, n2, isgn, *ip, i, j;
    double **a, *t, *w;
    
    nthread = ((fft2d_arg_t *) p)->nthread;
    n0 = ((fft2d_arg_t *) p)->n0;
    n1 = ((fft2d_arg_t *) p)->n1;
    n2 = ((fft2d_arg_t *) p)->n2;
    isgn = ((fft2d_arg_t *) p)->isgn;
    a = ((fft2d_arg_t *) p)->a;
    t = ((fft2d_arg_t *) p)->t;
    ip = ((fft2d_arg_t *) p)->ip;
    w = ((fft2d_arg_t *) p)->w;
    if (n2 > 4 * nthread) {
        for (j = 8 * n0; j < n2; j += 8 * nthread) {
            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 * nthread) {
        for (i = 0; i < n1; i++) {
            t[2 * i] = a[i][4 * n0];
            t[2 * i + 1] = a[i][4 * n0 + 1];
            t[2 * n1 + 2 * i] = a[i][4 * n0 + 2];
            t[2 * n1 + 2 * i + 1] = a[i][4 * n0 + 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][4 * n0] = t[2 * i];
            a[i][4 * n0 + 1] = t[2 * i + 1];
            a[i][4 * n0 + 2] = t[2 * n1 + 2 * i];
            a[i][4 * n0 + 3] = t[2 * n1 + 2 * i + 1];
        }
    } else if (n2 == 2 * nthread) {
        for (i = 0; i < n1; i++) {
            t[2 * i] = a[i][2 * n0];
            t[2 * i + 1] = a[i][2 * n0 + 1];
        }
        cdft(2 * n1, isgn, t, ip, w);
        for (i = 0; i < n1; i++) {
            a[i][2 * n0] = t[2 * i];
            a[i][2 * n0 + 1] = t[2 * i + 1];
        }
    }
    return (void *) 0;
}


void *ddxt2d0_th(void *p)
{
    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 nthread, n0, n1, n2, ics, isgn, *ip, i;
    double **a, *w;
    
    nthread = ((fft2d_arg_t *) p)->nthread;
    n0 = ((fft2d_arg_t *) p)->n0;
    n1 = ((fft2d_arg_t *) p)->n1;
    n2 = ((fft2d_arg_t *) p)->n2;
    ics = ((fft2d_arg_t *) p)->ic;
    isgn = ((fft2d_arg_t *) p)->isgn;
    a = ((fft2d_arg_t *) p)->a;
    ip = ((fft2d_arg_t *) p)->ip;
    w = ((fft2d_arg_t *) p)->w;
    if (ics == 0) {
        for (i = n0; i < n1; i += nthread) {
            ddct(n2, isgn, a[i], ip, w);
        }
    } else {
        for (i = n0; i < n1; i += nthread) {
            ddst(n2, isgn, a[i], ip, w);
        }
    }
    return (void *) 0;
}


void *ddxt2d_th(void *p)
{
    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 nthread, n0, n1, n2, ics, isgn, *ip, i, j;
    double **a, *t, *w;
    
    nthread = ((fft2d_arg_t *) p)->nthread;
    n0 = ((fft2d_arg_t *) p)->n0;
    n1 = ((fft2d_arg_t *) p)->n1;
    n2 = ((fft2d_arg_t *) p)->n2;
    ics = ((fft2d_arg_t *) p)->ic;
    isgn = ((fft2d_arg_t *) p)->isgn;
    a = ((fft2d_arg_t *) p)->a;
    t = ((fft2d_arg_t *) p)->t;
    ip = ((fft2d_arg_t *) p)->ip;
    w = ((fft2d_arg_t *) p)->w;
    if (n2 > 2 * nthread) {
        for (j = 4 * n0; j < n2; j += 4 * nthread) {
            for (i = 0; i < n1; i++) {
                t[i] = a[i][j];
                t[n1 + i] = a[i][j + 1];
                t[2 * n1 + i] = a[i][j + 2];
                t[3 * n1 + i] = a[i][j + 3];
            }
            if (ics == 0) {
                ddct(n1, isgn, t, ip, w);
                ddct(n1, isgn, &t[n1], ip, w);
                ddct(n1, isgn, &t[2 * n1], ip, w);
                ddct(n1, isgn, &t[3 * n1], ip, w);
            } else {
                ddst(n1, isgn, t, ip, w);
                ddst(n1, isgn, &t[n1], ip, w);
                ddst(n1, isgn, &t[2 * n1], ip, w);
                ddst(n1, isgn, &t[3 * n1], ip, w);
            }
            for (i = 0; i < n1; i++) {
                a[i][j] = t[i];
                a[i][j + 1] = t[n1 + i];
                a[i][j + 2] = t[2 * n1 + i];
                a[i][j + 3] = t[3 * n1 + i];
            }
        }
    } else if (n2 == 2 * nthread) {
        for (i = 0; i < n1; i++) {
            t[i] = a[i][2 * n0];
            t[n1 + i] = a[i][2 * n0 + 1];
        }
        if (ics == 0) {
            ddct(n1, isgn, t, ip, w);
            ddct(n1, isgn, &t[n1], ip, w);
        } else {
            ddst(n1, isgn, t, ip, w);
            ddst(n1, isgn, &t[n1], ip, w);
        }
        for (i = 0; i < n1; i++) {
            a[i][2 * n0] = t[i];
            a[i][2 * n0 + 1] = t[n1 + i];
        }
    } else if (n2 == nthread) {
        for (i = 0; i < n1; i++) {
            t[i] = a[i][n0];
        }
        if (ics == 0) {
            ddct(n1, isgn, t, ip, w);
        } else {
            ddst(n1, isgn, t, ip, w);
        }
        for (i = 0; i < n1; i++) {
            a[i][n0] = t[i];
        }
    }
    return (void *) 0;
}
#endif /* USE_FFT2D_THREADS */