fftsg_h.c

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

        cftleaf(m, isplt, &a[j - m]);
    }
    return (void *) 0;
}
#endif /* USE_CDFT_THREADS */


void cftrec4(int n, double *a)
{
    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;
    
    m = n;
    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]);
    }
}


int cfttree(int n, int j, int k, double *a)
{
    void cftmdl1(int n, double *a);
    void cftmdl2(int n, double *a);
    int i, isplt, m;
    
    if ((k & 3) != 0) {
        isplt = k & 1;
        if (isplt != 0) {
            cftmdl1(n, &a[j - n]);
        } else {
            cftmdl2(n, &a[j - n]);
        }
    } else {
        m = n;
        for (i = k; (i & 3) == 0; i >>= 2) {
            m <<= 2;
        }
        isplt = i & 1;
        if (isplt != 0) {
            while (m > 128) {
                cftmdl1(m, &a[j - m]);
                m >>= 2;
            }
        } else {
            while (m > 128) {
                cftmdl2(m, &a[j - m]);
                m >>= 2;
            }
        }
    }
    return isplt;
}


void cftleaf(int n, int isplt, double *a)
{
    void cftmdl1(int n, double *a);
    void cftmdl2(int n, double *a);
    void cftf161(double *a);
    void cftf162(double *a);
    void cftf081(double *a);
    void cftf082(double *a);
    
    if (n == 512) {
        cftmdl1(128, a);
        cftf161(a);
        cftf162(&a[32]);
        cftf161(&a[64]);
        cftf161(&a[96]);
        cftmdl2(128, &a[128]);
        cftf161(&a[128]);
        cftf162(&a[160]);
        cftf161(&a[192]);
        cftf162(&a[224]);
        cftmdl1(128, &a[256]);
        cftf161(&a[256]);
        cftf162(&a[288]);
        cftf161(&a[320]);
        cftf161(&a[352]);
        if (isplt != 0) {
            cftmdl1(128, &a[384]);
            cftf161(&a[480]);
        } else {
            cftmdl2(128, &a[384]);
            cftf162(&a[480]);
        }
        cftf161(&a[384]);
        cftf162(&a[416]);
        cftf161(&a[448]);
    } else {
        cftmdl1(64, a);
        cftf081(a);
        cftf082(&a[16]);
        cftf081(&a[32]);
        cftf081(&a[48]);
        cftmdl2(64, &a[64]);
        cftf081(&a[64]);
        cftf082(&a[80]);
        cftf081(&a[96]);
        cftf082(&a[112]);
        cftmdl1(64, &a[128]);
        cftf081(&a[128]);
        cftf082(&a[144]);
        cftf081(&a[160]);
        cftf081(&a[176]);
        if (isplt != 0) {
            cftmdl1(64, &a[192]);
            cftf081(&a[240]);
        } else {
            cftmdl2(64, &a[192]);
            cftf082(&a[240]);
        }
        cftf081(&a[192]);
        cftf082(&a[208]);
        cftf081(&a[224]);
    }
}


void cftmdl1(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;
}


void cftmdl2(int n, double *a)
{
    int i, i0, j, j0, j1, j2, j3, m, mh;
    double ew, w1r, w1i, wn4r, wk1r, wk1i, wk3r, wk3i, 
        wl1r, wl1i, wl3r, wl3i, wd1r, wd1i, wd3r, wd3i, 
        we1r, we1i, we3r, we3i, ss1, ss3;
    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y2r, y2i;
    
    mh = n >> 3;
    m = 2 * mh;
    wn4r = WR5000;
    j1 = m;
    j2 = j1 + m;
    j3 = j2 + m;
    x0r = a[0] - a[j2 + 1];
    x0i = a[1] + a[j2];
    x1r = a[0] + a[j2 + 1];
    x1i = a[1] - a[j2];
    x2r = a[j1] - a[j3 + 1];
    x2i = a[j1 + 1] + a[j3];
    x3r = a[j1] + a[j3 + 1];
    x3i = a[j1 + 1] - a[j3];
    y0r = wn4r * (x2r - x2i);
    y0i = wn4r * (x2i + x2r);
    a[0] = x0r + y0r;
    a[1] = x0i + y0i;
    a[j1] = x0r - y0r;
    a[j1 + 1] = x0i - y0i;
    y0r = wn4r * (x3r - x3i);
    y0i = wn4r * (x3i + x3r);
    a[j2] = x1r - y0i;
    a[j2 + 1] = x1i + y0r;
    a[j3] = x1r + y0i;
    a[j3 + 1] = x1i - y0r;
    wl1r = 1;
    wl1i = 0;
    wl3r = 1;
    wl3i = 0;
    we1r = wn4r;
    we1i = wn4r;
    we3r = -wn4r;
    we3i = -wn4r;
    ew = M_PI_2 / (2 * m);
    w1r = cos(2 * ew);
    w1i = sin(2 * ew);
    wk1r = w1r;
    wk1i = w1i;
    wd1r = wn4r * (w1r - w1i);
    wd1i = wn4r * (w1i + w1r);
    ss1 = 2 * w1i;
    wk3i = 2 * ss1 * wk1r;
    wk3r = wk1r - wk3i * wk1i;
    wk3i = wk1i - wk3i * wk1r;
    ss3 = 2 * wk3i;
    wd3r = -wn4r * (wk3r - wk3i);
    wd3i = -wn4r * (wk3i + wk3r);
    i = 0;
    for (;;) {
        i0 = i + 4 * CDFT_LOOP_DIV;
        if (i0 > mh - 4) {
            i0 = mh - 4;
        }
        for (j = i + 2; j < i0; j += 4) {
            wl1r -= ss1 * wk1i;
            wl1i += ss1 * wk1r;
            wl3r -= ss3 * wk3i;
            wl3i += ss3 * wk3r;
            we1r -= ss1 * wd1i;
            we1i += ss1 * wd1r;
            we3r -= ss3 * wd3i;
            we3i += ss3 * wd3r;
            j1 = j + m;
            j2 = j1 + m;
            j3 = j2 + m;
            x0r = a[j] - a[j2 + 1];
            x0i = a[j + 1] + a[j2];
            x1r = a[j] + a[j2 + 1];
            x1i = a[j + 1] - a[j2];
            x2r = a[j1] - a[j3 + 1];
            x2i = a[j1 + 1] + a[j3];
            x3r = a[j1] + a[j3 + 1];
            x3i = a[j1 + 1] - a[j3];
            y0r = wk1r * x0r - wk1i * x0i;
            y0i = wk1r * x0i + wk1i * x0r;
            y2r = wd1r * x2r - wd1i * x2i;
            y2i = wd1r * x2i + wd1i * x2r;
            a[j] = y0r + y2r;
            a[j + 1] = y0i + y2i;
            a[j1] = y0r - y2r;
            a[j1 + 1] = y0i - y2i;
            y0r = wk3r * x1r + wk3i * x1i;
            y0i = wk3r * x1i - wk3i * x1r;
            y2r = wd3r * x3r + wd3i * x3i;
            y2i = wd3r * x3i - wd3i * x3r;
            a[j2] = y0r + y2r;
            a[j2 + 1] = y0i + y2i;
            a[j3] = y0r - y2r;
            a[j3 + 1] = y0i - y2i;
            x0r = a[j + 2] - a[j2 + 3];
            x0i = a[j + 3] + a[j2 + 2];
            x1r = a[j + 2] + a[j2 + 3];
            x1i = a[j + 3] - a[j2 + 2];
            x2r = a[j1 + 2] - a[j3 + 3];
            x2i = a[j1 + 3] + a[j3 + 2];
            x3r = a[j1 + 2] + a[j3 + 3];
            x3i = a[j1 + 3] - a[j3 + 2];
            y0r = wl1r * x0r - wl1i * x0i;
            y0i = wl1r * x0i + wl1i * x0r;
            y2r = we1r * x2r - we1i * x2i;
            y2i = we1r * x2i + we1i * x2r;
            a[j + 2] = y0r + y2r;
            a[j + 3] = y0i + y2i;
            a[j1 + 2] = y0r - y2r;
            a[j1 + 3] = y0i - y2i;
            y0r = wl3r * x1r + wl3i * x1i;
            y0i = wl3r * x1i - wl3i * x1r;
            y2r = we3r * x3r + we3i * x3i;
            y2i = we3r * x3i - we3i * x3r;
            a[j2 + 2] = y0r + y2r;
            a[j2 + 3] = y0i + y2i;
            a[j3 + 2] = y0r - y2r;
            a[j3 + 3] = y0i - y2i;
            j0 = m - j;
            j1 = j0 + m;
            j2 = j1 + m;
            j3 = j2 + m;
            x0r = a[j0] - a[j2 + 1];
            x0i = a[j0 + 1] + a[j2];
            x1r = a[j0] + a[j2 + 1];
            x1i = a[j0 + 1] - a[j2];
            x2r = a[j1] - a[j3 + 1];
            x2i = a[j1 + 1] + a[j3];
            x3r = a[j1] + a[j3 + 1];
            x3i = a[j1 + 1] - a[j3];