fftsg_h.c

This is a package to calculate Discrete Fourier/Cosine/Sine Transforms of 1-dimensional sequences of

}


void bitrv208(double *a)
{
    double x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i;
    
    x1r = a[2];
    x1i = a[3];
    x3r = a[6];
    x3i = a[7];
    x4r = a[8];
    x4i = a[9];
    x6r = a[12];
    x6i = a[13];
    a[2] = x4r;
    a[3] = x4i;
    a[6] = x6r;
    a[7] = x6i;
    a[8] = x1r;
    a[9] = x1i;
    a[12] = x3r;
    a[13] = x3i;
}


void bitrv208neg(double *a)
{
    double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
        x5r, x5i, x6r, x6i, x7r, x7i;
    
    x1r = a[2];
    x1i = a[3];
    x2r = a[4];
    x2i = a[5];
    x3r = a[6];
    x3i = a[7];
    x4r = a[8];
    x4i = a[9];
    x5r = a[10];
    x5i = a[11];
    x6r = a[12];
    x6i = a[13];
    x7r = a[14];
    x7i = a[15];
    a[2] = x7r;
    a[3] = x7i;
    a[4] = x3r;
    a[5] = x3i;
    a[6] = x5r;
    a[7] = x5i;
    a[8] = x1r;
    a[9] = x1i;
    a[10] = x6r;
    a[11] = x6i;
    a[12] = x2r;
    a[13] = x2i;
    a[14] = x4r;
    a[15] = x4i;
}


void bitrv1(int n, double *a)
{
    int j0, k0, j1, k1, l, m, i, j, k, nh;
    double x;
    
    nh = n >> 1;
    x = a[1];
    a[1] = a[nh];
    a[nh] = x;
    m = 2;
    for (l = n >> 2; l > 2; l >>= 2) {
        m <<= 1;
    }
    if (l == 2) {
        j1 = m + 1;
        k1 = m + nh;
        x = a[j1];
        a[j1] = a[k1];
        a[k1] = x;
        j0 = 0;
        for (k0 = 2; k0 < m; k0 += 2) {
            for (i = nh >> 1; i > (j0 ^= i); i >>= 1);
            k = k0;
            for (j = j0; j < j0 + k0; j += 2) {
                x = a[j];
                a[j] = a[k];
                a[k] = x;
                j1 = j + m;
                k1 = k + m;
                x = a[j1];
                a[j1] = a[k1];
                a[k1] = x;
                j1 += nh;
                k1++;
                x = a[j1];
                a[j1] = a[k1];
                a[k1] = x;
                j1 -= m;
                k1 -= m;
                x = a[j1];
                a[j1] = a[k1];
                a[k1] = x;
                j1++;
                k1 += nh;
                x = a[j1];
                a[j1] = a[k1];
                a[k1] = x;
                j1 += m;
                k1 += m;
                x = a[j1];
                a[j1] = a[k1];
                a[k1] = x;
                j1 -= nh;
                k1--;
                x = a[j1];
                a[j1] = a[k1];
                a[k1] = x;
                j1 -= m;
                k1 -= m;
                x = a[j1];
                a[j1] = a[k1];
                a[k1] = x;
                for (i = nh >> 1; i > (k ^= i); i >>= 1);
            }
            k1 = j0 + k0;
            j1 = k1 + 1;
            k1 += nh;
            x = a[j1];
            a[j1] = a[k1];
            a[k1] = x;
            j1 += m;
            k1 += m;
            x = a[j1];
            a[j1] = a[k1];
            a[k1] = x;
        }
    } else {
        j0 = 0;
        for (k0 = 2; k0 < m; k0 += 2) {
            for (i = nh >> 1; i > (j0 ^= i); i >>= 1);
            k = k0;
            for (j = j0; j < j0 + k0; j += 2) {
                x = a[j];
                a[j] = a[k];
                a[k] = x;
                j1 = j + nh;
                k1 = k + 1;
                x = a[j1];
                a[j1] = a[k1];
                a[k1] = x;
                j1++;
                k1 += nh;
                x = a[j1];
                a[j1] = a[k1];
                a[k1] = x;
                j1 -= nh;
                k1--;
                x = a[j1];
                a[j1] = a[k1];
                a[k1] = x;
                for (i = nh >> 1; i > (k ^= i); i >>= 1);
            }
            k1 = j0 + k0;
            j1 = k1 + 1;
            k1 += nh;
            x = a[j1];
            a[j1] = a[k1];
            a[k1] = x;
        }
    }
}


void cftb1st(int n, double *a)
{
    int i, i0, j, j0, j1, j2, j3, m, mh;
    double ew, w1r, w1i, wk1r, wk1i, wk3r, wk3i, 
        wd1r, wd1i, wd3r, wd3i, ss1, ss3;
    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
    
    mh = n >> 3;
    m = 2 * mh;
    j1 = m;
    j2 = j1 + m;
    j3 = j2 + m;
    x0r = a[0] + a[j2];
    x0i = -a[1] - a[j2 + 1];
    x1r = a[0] - a[j2];
    x1i = -a[1] + a[j2 + 1];
    x2r = a[j1] + a[j3];
    x2i = a[j1 + 1] + a[j3 + 1];
    x3r = a[j1] - a[j3];
    x3i = a[j1 + 1] - a[j3 + 1];
    a[0] = x0r + x2r;
    a[1] = x0i - x2i;
    a[j1] = x0r - x2r;
    a[j1 + 1] = x0i + x2i;
    a[j2] = x1r + x3i;
    a[j2 + 1] = x1i + x3r;
    a[j3] = x1r - x3i;
    a[j3 + 1] = x1i - x3r;
    wd1r = 1;
    wd1i = 0;
    wd3r = 1;
    wd3i = 0;
    ew = M_PI_2 / m;
    w1r = cos(2 * ew);
    w1i = sin(2 * ew);
    wk1r = w1r;
    wk1i = w1i;
    ss1 = 2 * w1i;
    wk3i = 2 * ss1 * wk1r;
    wk3r = wk1r - wk3i * wk1i;
    wk3i = wk1i - wk3i * wk1r;
    ss3 = 2 * wk3i;
    i = 0;
    for (;;) {
        i0 = i + 4 * CDFT_LOOP_DIV;
        if (i0 > mh - 4) {
            i0 = mh - 4;
        }
        for (j = i + 2; j < i0; j += 4) {
            wd1r -= ss1 * wk1i;
            wd1i += ss1 * wk1r;
            wd3r -= ss3 * wk3i;
            wd3i += ss3 * wk3r;
            j1 = j + m;
            j2 = j1 + m;
            j3 = j2 + m;
            x0r = a[j] + a[j2];
            x0i = -a[j + 1] - a[j2 + 1];
            x1r = a[j] - a[j2];
            x1i = -a[j + 1] + a[j2 + 1];
            x2r = a[j1] + a[j3];
            x2i = a[j1 + 1] + a[j3 + 1];
            x3r = a[j1] - a[j3];
            x3i = a[j1 + 1] - a[j3 + 1];
            a[j] = x0r + x2r;
            a[j + 1] = x0i - x2i;
            a[j1] = x0r - x2r;
            a[j1 + 1] = x0i + x2i;
            x0r = x1r + x3i;
            x0i = x1i + x3r;
            a[j2] = wk1r * x0r - wk1i * x0i;
            a[j2 + 1] = wk1r * x0i + wk1i * x0r;
            x0r = x1r - x3i;
            x0i = x1i - x3r;
            a[j3] = wk3r * x0r + wk3i * x0i;
            a[j3 + 1] = wk3r * x0i - wk3i * x0r;
            x0r = a[j + 2] + a[j2 + 2];
            x0i = -a[j + 3] - a[j2 + 3];
            x1r = a[j + 2] - a[j2 + 2];
            x1i = -a[j + 3] + a[j2 + 3];
            x2r = a[j1 + 2] + a[j3 + 2];
            x2i = a[j1 + 3] + a[j3 + 3];
            x3r = a[j1 + 2] - a[j3 + 2];
            x3i = a[j1 + 3] - a[j3 + 3];
            a[j + 2] = x0r + x2r;
            a[j + 3] = x0i - x2i;
            a[j1 + 2] = x0r - x2r;
            a[j1 + 3] = x0i + x2i;
            x0r = x1r + x3i;
            x0i = x1i + x3r;
            a[j2 + 2] = wd1r * x0r - wd1i * x0i;
            a[j2 + 3] = wd1r * x0i + wd1i * x0r;
            x0r = x1r - x3i;
            x0i = x1i - x3r;
            a[j3 + 2] = wd3r * x0r + wd3i * x0i;
            a[j3 + 3] = wd3r * x0i - wd3i * x0r;
            j0 = m - j;
            j1 = j0 + m;
            j2 = j1 + m;
            j3 = j2 + m;
            x0r = a[j0] + a[j2];
            x0i = -a[j0 + 1] - a[j2 + 1];
            x1r = a[j0] - a[j2];
            x1i = -a[j0 + 1] + a[j2 + 1];
            x2r = a[j1] + a[j3];
            x2i = a[j1 + 1] + a[j3 + 1];
            x3r = a[j1] - a[j3];
            x3i = a[j1 + 1] - a[j3 + 1];
            a[j0] = x0r + x2r;
            a[j0 + 1] = x0i - x2i;
            a[j1] = x0r - x2r;
            a[j1 + 1] = x0i + x2i;
            x0r = x1r + x3i;
            x0i = x1i + x3r;
            a[j2] = wk1i * x0r - wk1r * x0i;
            a[j2 + 1] = wk1i * x0i + wk1r * x0r;
            x0r = x1r - x3i;
            x0i = x1i - x3r;
            a[j3] = wk3i * x0r + wk3r * x0i;
            a[j3 + 1] = wk3i * x0i - wk3r * x0r;
            x0r = a[j0 - 2] + a[j2 - 2];
            x0i = -a[j0 - 1] - a[j2 - 1];
            x1r = a[j0 - 2] - a[j2 - 2];
            x1i = -a[j0 - 1] + a[j2 - 1];
            x2r = a[j1 - 2] + a[j3 - 2];
            x2i = a[j1 - 1] + a[j3 - 1];
            x3r = a[j1 - 2] - a[j3 - 2];
            x3i = a[j1 - 1] - a[j3 - 1];
            a[j0 - 2] = x0r + x2r;
            a[j0 - 1] = x0i - x2i;
            a[j1 - 2] = x0r - x2r;
            a[j1 - 1] = x0i + x2i;
            x0r = x1r + x3i;
            x0i = x1i + x3r;
            a[j2 - 2] = wd1i * x0r - wd1r * x0i;
            a[j2 - 1] = wd1i * x0i + wd1r * x0r;
            x0r = x1r - x3i;
            x0i = x1i - x3r;
            a[j3 - 2] = wd3i * x0r + wd3r * x0i;
            a[j3 - 1] = wd3i * x0i - wd3r * x0r;
            wk1r -= ss1 * wd1i;
            wk1i += ss1 * wd1r;
            wk3r -= ss3 * wd3i;
            wk3i += ss3 * wd3r;
        }
        if (i0 == mh - 4) {
            break;
        }
        wd1r = cos(ew * i0);
        wd1i = sin(ew * i0);
        wd3i = 4 * wd1i * wd1r;
        wd3r = wd1r - wd3i * wd1i;
        wd3i = wd1i - wd3i * wd1r;
        wk1r = w1r * wd1r - w1i * wd1i;
        wk1i = w1r * wd1i + w1i * wd1r;
        wk3i = 4 * wk1i * wk1r;
        wk3r = wk1r - wk3i * wk1i;
        wk3i = wk1i - wk3i * wk1r;
        i = i0;
    }
    wd1r = WR5000;
    j0 = mh;
    j1 = j0 + m;
    j2 = j1 + m;
    j3 = j2 + m;
    x0r = a[j0 - 2] + a[j2 - 2];
    x0i = -a[j0 - 1] - a[j2 - 1];
    x1r = a[j0 - 2] - a[j2 - 2];
    x1i = -a[j0 - 1] + a[j2 - 1];
    x2r = a[j1 - 2] + a[j3 - 2];
    x2i = a[j1 - 1] + a[j3 - 1];
    x3r = a[j1 - 2] - a[j3 - 2];
    x3i = a[j1 - 1] - a[j3 - 1];
    a[j0 - 2] = x0r + x2r;
    a[j0 - 1] = x0i - x2i;
    a[j1 - 2] = x0r - x2r;
    a[j1 - 1] = x0i + x2i;
    x0r = x1r + x3i;
    x0i = x1i + x3r;
    a[j2 - 2] = wk1r * x0r - wk1i * x0i;
    a[j2 - 1] = wk1r * x0i + wk1i * x0r;
    x0r = x1r - x3i;
    x0i = x1i - x3r;
    a[j3 - 2] = wk3r * x0r + wk3i * x0i;
    a[j3 - 1] = wk3r * x0i - wk3i * x0r;
    x0r = a[j0] + a[j2];
    x0i = -a[j0 + 1] - a[j2 + 1];
    x1r = a[j0] - a[j2];
    x1i = -a[j0 + 1] + a[j2 + 1];
    x2r = a[j1] + a[j3];
    x2i = a[j1 + 1] + a[j3 + 1];
    x3r = a[j1] - a[j3];
    x3i = a[j1 + 1] - a[j3 + 1];
    a[j0] = x0r + x2r;
    a[j0 + 1] = x0i - x2i;
    a[j1] = x0r - x2r;
    a[j1 + 1] = x0i + x2i;
    x0r = x1r + x3i;
    x0i = x1i + x3r;
    a[j2] = wd1r * (x0r - x0i);
    a[j2 + 1] = wd1r * (x0i + x0r);
    x0r = x1r - x3i;
    x0i = x1i - x3r;
    a[j3] = -wd1r * (x0r + x0i);
    a[j3 + 1] = -wd1r * (x0i - x0r);
    x0r = a[j0 + 2] + a[j2 + 2];
    x0i = -a[j0 + 3] - a[j2 + 3];
    x1r = a[j0 + 2] - a[j2 + 2];
    x1i = -a[j0 + 3] + a[j2 + 3];
    x2r = a[j1 + 2] + a[j3 + 2];
    x2i = a[j1 + 3] + a[j3 + 3];
    x3r = a[j1 + 2] - a[j3 + 2];
    x3i = a[j1 + 3] - a[j3 + 3];
    a[j0 + 2] = x0r + x2r;
    a[j0 + 3] = x0i - x2i;
    a[j1 + 2] = x0r - x2r;
    a[j1 + 3] = x0i + x2i;
    x0r = x1r + x3i;
    x0i = x1i + x3r;
    a[j2 + 2] = wk1i * x0r - wk1r * x0i;
    a[j2 + 3] = wk1i * x0i + wk1r * x0r;
    x0r = x1r - x3i;
    x0i = x1i - x3r;
    a[j3 + 2] = wk3i * x0r + wk3r * x0i;
    a[j3 + 3] = wk3i * x0i - wk3r * x0r;
}


#ifdef USE_CDFT_THREADS
struct cdft_arg_st {
    int n0;
    int n;
    double *a;
};
typedef struct cdft_arg_st cdft_arg_t;


void cftrec4_th(int n, double *a)
{
    void *cftrec1_th(void *p);
    void *cftrec2_th(void *p);
    int i, idiv4, m, nthread;
    cdft_thread_t th[4];
    cdft_arg_t ag[4];
    
    nthread = 2;
    idiv4 = 0;
    m = n >> 1;
    if (n > CDFT_4THREADS_BEGIN_N) {
        nthread = 4;
        idiv4 = 1;
        m >>= 1;
    }
    for (i = 0; i < nthread; i++) {
        ag[i].n0 = n;
        ag[i].n = m;
        ag[i].a = &a[i * m];
        if (i != idiv4) {
            cdft_thread_create(&th[i], cftrec1_th, &ag[i]);
        } else {
            cdft_thread_create(&th[i], cftrec2_th, &ag[i]);
        }
    }
    for (i = 0; i < nthread; i++) {
        cdft_thread_wait(th[i]);
    }
}


void *cftrec1_th(void *p)
{
    int cfttree(int n, int j, int k, double *a);
    void cftleaf(int n, int isplt, double *a);
    void cftmdl1(int n, double *a);
    int isplt, j, k, m, n, n0;
    double *a;
    
    n0 = ((cdft_arg_t *) p)->n0;
    n = ((cdft_arg_t *) p)->n;
    a = ((cdft_arg_t *) p)->a;
    m = n0;
    while (m > 512) {
        m >>= 2;
        cftmdl1(m, &a[n - m]);
    }
    cftleaf(m, 1, &a[n - m]);
    k = 0;
    for (j = n - m; j > 0; j -= m) {
        k++;
        isplt = cfttree(m, j, k, a);
        cftleaf(m, isplt, &a[j - m]);
    }
    return (void *) 0;
}


void *cftrec2_th(void *p)
{
    int cfttree(int n, int j, int k, double *a);
    void cftleaf(int n, int isplt, double *a);
    void cftmdl2(int n, double *a);
    int isplt, j, k, m, n, n0;
    double *a;
    
    n0 = ((cdft_arg_t *) p)->n0;
    n = ((cdft_arg_t *) p)->n;
    a = ((cdft_arg_t *) p)->a;
    k = 1;
    m = n0;
    while (m > 512) {
        m >>= 2;
        k <<= 2;
        cftmdl2(m, &a[n - m]);
    }
    cftleaf(m, 0, &a[n - m]);
    k >>= 1;
    for (j = n - m; j > 0; j -= m) {
        k++;
        isplt = cfttree(m, j, k, a);