fft8g_h.c

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

            x0i = y1i + y5i;
            a[j + 2] = wk1r * x0r - wk1i * x0i;
            a[j + 3] = wk1r * x0i + wk1i * x0r;
            x0r = y1r - y5r;
            x0i = y1i - y5i;
            a[j + 10] = wk5r * x0r - wk5i * x0i;
            a[j + 11] = wk5r * x0i + wk5i * x0r;
            x0r = y3r - y7i;
            x0i = y3i + y7r;
            a[j + 6] = wk3r * x0r - wk3i * x0i;
            a[j + 7] = wk3r * x0i + wk3i * x0r;
            x0r = y3r + y7i;
            x0i = y3i - y7r;
            a[j + 14] = wk7r * x0r - wk7i * x0i;
            a[j + 15] = wk7r * x0i + wk7i * x0r;
            a[j] = y0r + y4r;
            a[j + 1] = y0i + y4i;
            x0r = y0r - y4r;
            x0i = y0i - y4i;
            a[j + 8] = wk4r * x0r - wk4i * x0i;
            a[j + 9] = wk4r * x0i + wk4i * x0r;
            x0r = y2r - y6i;
            x0i = y2i + y6r;
            a[j + 4] = wk2r * x0r - wk2i * x0i;
            a[j + 5] = wk2r * x0i + wk2i * x0r;
            x0r = y2r + y6i;
            x0i = y2i - y6r;
            a[j + 12] = wk6r * x0r - wk6i * x0i;
            a[j + 13] = wk6r * x0i + wk6i * x0r;
        }
    }
}


void cftmdl(int n, int l, double *a)
{
    int j, j1, j2, j3, j4, j5, j6, j7, k, kj, kr, m;
    double ew, wn4r, wtmp, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i, 
        wk4r, wk4i, wk5r, wk5i, wk6r, wk6i, wk7r, wk7i;
    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
        y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
        y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
    
    m = l << 3;
    wn4r = WR5000;
    for (j = 0; j < l; j += 2) {
        j1 = j + l;
        j2 = j1 + l;
        j3 = j2 + l;
        j4 = j3 + l;
        j5 = j4 + l;
        j6 = j5 + l;
        j7 = j6 + l;
        x0r = a[j] + a[j1];
        x0i = a[j + 1] + a[j1 + 1];
        x1r = a[j] - a[j1];
        x1i = a[j + 1] - a[j1 + 1];
        x2r = a[j2] + a[j3];
        x2i = a[j2 + 1] + a[j3 + 1];
        x3r = a[j2] - a[j3];
        x3i = a[j2 + 1] - a[j3 + 1];
        y0r = x0r + x2r;
        y0i = x0i + x2i;
        y2r = x0r - x2r;
        y2i = x0i - x2i;
        y1r = x1r - x3i;
        y1i = x1i + x3r;
        y3r = x1r + x3i;
        y3i = x1i - x3r;
        x0r = a[j4] + a[j5];
        x0i = a[j4 + 1] + a[j5 + 1];
        x1r = a[j4] - a[j5];
        x1i = a[j4 + 1] - a[j5 + 1];
        x2r = a[j6] + a[j7];
        x2i = a[j6 + 1] + a[j7 + 1];
        x3r = a[j6] - a[j7];
        x3i = a[j6 + 1] - a[j7 + 1];
        y4r = x0r + x2r;
        y4i = x0i + x2i;
        y6r = x0r - x2r;
        y6i = x0i - x2i;
        x0r = x1r - x3i;
        x0i = x1i + x3r;
        x2r = x1r + x3i;
        x2i = x1i - x3r;
        y5r = wn4r * (x0r - x0i);
        y5i = wn4r * (x0r + x0i);
        y7r = wn4r * (x2r - x2i);
        y7i = wn4r * (x2r + x2i);
        a[j1] = y1r + y5r;
        a[j1 + 1] = y1i + y5i;
        a[j5] = y1r - y5r;
        a[j5 + 1] = y1i - y5i;
        a[j3] = y3r - y7i;
        a[j3 + 1] = y3i + y7r;
        a[j7] = y3r + y7i;
        a[j7 + 1] = y3i - y7r;
        a[j] = y0r + y4r;
        a[j + 1] = y0i + y4i;
        a[j4] = y0r - y4r;
        a[j4 + 1] = y0i - y4i;
        a[j2] = y2r - y6i;
        a[j2 + 1] = y2i + y6r;
        a[j6] = y2r + y6i;
        a[j6 + 1] = y2i - y6r;
    }
    if (m < n) {
        wk1r = WR2500;
        wk1i = WI2500;
        for (j = m; j < l + m; j += 2) {
            j1 = j + l;
            j2 = j1 + l;
            j3 = j2 + l;
            j4 = j3 + l;
            j5 = j4 + l;
            j6 = j5 + l;
            j7 = j6 + l;
            x0r = a[j] + a[j1];
            x0i = a[j + 1] + a[j1 + 1];
            x1r = a[j] - a[j1];
            x1i = a[j + 1] - a[j1 + 1];
            x2r = a[j2] + a[j3];
            x2i = a[j2 + 1] + a[j3 + 1];
            x3r = a[j2] - a[j3];
            x3i = a[j2 + 1] - a[j3 + 1];
            y0r = x0r + x2r;
            y0i = x0i + x2i;
            y2r = x0r - x2r;
            y2i = x0i - x2i;
            y1r = x1r - x3i;
            y1i = x1i + x3r;
            y3r = x1r + x3i;
            y3i = x1i - x3r;
            x0r = a[j4] + a[j5];
            x0i = a[j4 + 1] + a[j5 + 1];
            x1r = a[j4] - a[j5];
            x1i = a[j4 + 1] - a[j5 + 1];
            x2r = a[j6] + a[j7];
            x2i = a[j6 + 1] + a[j7 + 1];
            x3r = a[j6] - a[j7];
            x3i = a[j6 + 1] - a[j7 + 1];
            y4r = x0r + x2r;
            y4i = x0i + x2i;
            y6r = x0r - x2r;
            y6i = x0i - x2i;
            x0r = x1r - x3i;
            x0i = x1i + x3r;
            x2r = x1r + x3i;
            x2i = x3r - x1i;
            y5r = wk1i * x0r - wk1r * x0i;
            y5i = wk1i * x0i + wk1r * x0r;
            y7r = wk1r * x2r + wk1i * x2i;
            y7i = wk1r * x2i - wk1i * x2r;
            x0r = wk1r * y1r - wk1i * y1i;
            x0i = wk1r * y1i + wk1i * y1r;
            a[j1] = x0r + y5r;
            a[j1 + 1] = x0i + y5i;
            a[j5] = y5i - x0i;
            a[j5 + 1] = x0r - y5r;
            x0r = wk1i * y3r - wk1r * y3i;
            x0i = wk1i * y3i + wk1r * y3r;
            a[j3] = x0r - y7r;
            a[j3 + 1] = x0i + y7i;
            a[j7] = y7i - x0i;
            a[j7 + 1] = x0r + y7r;
            a[j] = y0r + y4r;
            a[j + 1] = y0i + y4i;
            a[j4] = y4i - y0i;
            a[j4 + 1] = y0r - y4r;
            x0r = y2r - y6i;
            x0i = y2i + y6r;
            a[j2] = wn4r * (x0r - x0i);
            a[j2 + 1] = wn4r * (x0i + x0r);
            x0r = y6r - y2i;
            x0i = y2r + y6i;
            a[j6] = wn4r * (x0r - x0i);
            a[j6 + 1] = wn4r * (x0i + x0r);
        }
        ew = M_PI_2 / n;
        kr = n >> 2;
        for (k = 2 * m; k < n; k += m) {
            for (kj = n >> 2; kj > (kr ^= kj); kj >>= 1);
            wk1r = cos(ew * kr);
            wk1i = sin(ew * kr);
            wk2r = 1 - 2 * wk1i * wk1i;
            wk2i = 2 * wk1i * wk1r;
            wtmp = 2 * wk2i;
            wk3r = wk1r - wtmp * wk1i;
            wk3i = wtmp * wk1r - wk1i;
            wk4r = 1 - wtmp * wk2i;
            wk4i = wtmp * wk2r;
            wtmp = 2 * wk4i;
            wk5r = wk3r - wtmp * wk1i;
            wk5i = wtmp * wk1r - wk3i;
            wk6r = wk2r - wtmp * wk2i;
            wk6i = wtmp * wk2r - wk2i;
            wk7r = wk1r - wtmp * wk3i;
            wk7i = wtmp * wk3r - wk1i;
            for (j = k; j < l + k; j += 2) {
                j1 = j + l;
                j2 = j1 + l;
                j3 = j2 + l;
                j4 = j3 + l;
                j5 = j4 + l;
                j6 = j5 + l;
                j7 = j6 + l;
                x0r = a[j] + a[j1];
                x0i = a[j + 1] + a[j1 + 1];
                x1r = a[j] - a[j1];
                x1i = a[j + 1] - a[j1 + 1];
                x2r = a[j2] + a[j3];
                x2i = a[j2 + 1] + a[j3 + 1];
                x3r = a[j2] - a[j3];
                x3i = a[j2 + 1] - a[j3 + 1];
                y0r = x0r + x2r;
                y0i = x0i + x2i;
                y2r = x0r - x2r;
                y2i = x0i - x2i;
                y1r = x1r - x3i;
                y1i = x1i + x3r;
                y3r = x1r + x3i;
                y3i = x1i - x3r;
                x0r = a[j4] + a[j5];
                x0i = a[j4 + 1] + a[j5 + 1];
                x1r = a[j4] - a[j5];
                x1i = a[j4 + 1] - a[j5 + 1];
                x2r = a[j6] + a[j7];
                x2i = a[j6 + 1] + a[j7 + 1];
                x3r = a[j6] - a[j7];
                x3i = a[j6 + 1] - a[j7 + 1];
                y4r = x0r + x2r;
                y4i = x0i + x2i;
                y6r = x0r - x2r;
                y6i = x0i - x2i;
                x0r = x1r - x3i;
                x0i = x1i + x3r;
                x2r = x1r + x3i;
                x2i = x1i - x3r;
                y5r = wn4r * (x0r - x0i);
                y5i = wn4r * (x0r + x0i);
                y7r = wn4r * (x2r - x2i);
                y7i = wn4r * (x2r + x2i);
                x0r = y1r + y5r;
                x0i = y1i + y5i;
                a[j1] = wk1r * x0r - wk1i * x0i;
                a[j1 + 1] = wk1r * x0i + wk1i * x0r;
                x0r = y1r - y5r;
                x0i = y1i - y5i;
                a[j5] = wk5r * x0r - wk5i * x0i;
                a[j5 + 1] = wk5r * x0i + wk5i * x0r;
                x0r = y3r - y7i;
                x0i = y3i + y7r;
                a[j3] = wk3r * x0r - wk3i * x0i;
                a[j3 + 1] = wk3r * x0i + wk3i * x0r;
                x0r = y3r + y7i;
                x0i = y3i - y7r;
                a[j7] = wk7r * x0r - wk7i * x0i;
                a[j7 + 1] = wk7r * x0i + wk7i * x0r;
                a[j] = y0r + y4r;
                a[j + 1] = y0i + y4i;
                x0r = y0r - y4r;
                x0i = y0i - y4i;
                a[j4] = wk4r * x0r - wk4i * x0i;
                a[j4 + 1] = wk4r * x0i + wk4i * x0r;
                x0r = y2r - y6i;
                x0i = y2i + y6r;
                a[j2] = wk2r * x0r - wk2i * x0i;
                a[j2 + 1] = wk2r * x0i + wk2i * x0r;
                x0r = y2r + y6i;
                x0i = y2i - y6r;
                a[j6] = wk6r * x0r - wk6i * x0i;
                a[j6 + 1] = wk6r * x0i + wk6i * x0r;
            }
        }
    }
}


void rftfsub(int n, double *a)
{
    int i, i0, j, k;
    double ec, w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi;
    
    ec = 2 * M_PI_2 / n;
    wkr = 0;
    wki = 0;
    wdi = cos(ec);
    wdr = sin(ec);
    wdi *= wdr;
    wdr *= wdr;
    w1r = 1 - 2 * wdr;
    w1i = 2 * wdi;
    ss = 2 * w1i;
    i = n >> 1;
    for (;;) {
        i0 = i - 4 * RDFT_LOOP_DIV;
        if (i0 < 4) {
            i0 = 4;
        }
        for (j = i - 4; j >= i0; j -= 4) {
            k = n - j;
            xr = a[j + 2] - a[k - 2];
            xi = a[j + 3] + a[k - 1];
            yr = wdr * xr - wdi * xi;
            yi = wdr * xi + wdi * xr;
            a[j + 2] -= yr;
            a[j + 3] -= yi;
            a[k - 2] += yr;
            a[k - 1] -= yi;
            wkr += ss * wdi;
            wki += ss * (0.5 - wdr);
            xr = a[j] - a[k];
            xi = a[j + 1] + a[k + 1];
            yr = wkr * xr - wki * xi;
            yi = wkr * xi + wki * xr;
            a[j] -= yr;
            a[j + 1] -= yi;
            a[k] += yr;
            a[k + 1] -= yi;
            wdr += ss * wki;
            wdi += ss * (0.5 - wkr);
        }
        if (i0 == 4) {
            break;
        }
        wkr = 0.5 * sin(ec * i0);
        wki = 0.5 * cos(ec * i0);
        wdr = 0.5 - (wkr * w1r - wki * w1i);
        wdi = wkr * w1i + wki * w1r;
        wkr = 0.5 - wkr;
        i = i0;
    }
    xr = a[2] - a[n - 2];
    xi = a[3] + a[n - 1];
    yr = wdr * xr - wdi * xi;
    yi = wdr * xi + wdi * xr;
    a[2] -= yr;
    a[3] -= yi;
    a[n - 2] += yr;
    a[n - 1] -= yi;
}


void rftbsub(int n, double *a)
{
    int i, i0, j, k;
    double ec, w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi;
    
    ec = 2 * M_PI_2 / n;
    wkr = 0;
    wki = 0;
    wdi = cos(ec);
    wdr = sin(ec);
    wdi *= wdr;
    wdr *= wdr;
    w1r = 1 - 2 * wdr;
    w1i = 2 * wdi;
    ss = 2 * w1i;
    i = n >> 1;
    a[i + 1] = -a[i + 1];
    for (;;) {
        i0 = i - 4 * RDFT_LOOP_DIV;
        if (i0 < 4) {
            i0 = 4;
        }
        for (j = i - 4; j >= i0; j -= 4) {
            k = n - j;
            xr = a[j + 2] - a[k - 2];
            xi = a[j + 3] + a[k - 1];
            yr = wdr * xr + wdi * xi;
            yi = wdr * xi - wdi * xr;
            a[j + 2] -= yr;
            a[j + 3] = yi - a[j + 3];
            a[k - 2] += yr;
            a[k - 1] = yi - a[k - 1];
            wkr += ss * wdi;
            wki += ss * (0.5 - wdr);
            xr = a[j] - a[k];
            xi = a[j + 1] + a[k + 1];
            yr = wkr * xr + wki * xi;
            yi = wkr * xi - wki * xr;
            a[j] -= yr;
            a[j + 1] = yi - a[j + 1];
            a[k] += yr;
            a[k + 1] = yi - a[k + 1];
            wdr += ss * wki;
            wdi += ss * (0.5 - wkr);
        }
        if (i0 == 4) {
            break;
        }
        wkr = 0.5 * sin(ec * i0);
        wki = 0.5 * cos(ec * i0);
        wdr = 0.5 - (wkr * w1r - wki * w1i);
        wdi = wkr * w1i + wki * w1r;
        wkr = 0.5 - wkr;
        i = i0;
    }
    xr = a[2] - a[n - 2];
    xi = a[3] + a[n - 1];
    yr = wdr * xr + wdi * xi;
    yi = wdr * xi - wdi * xr;
    a[2] -= yr;
    a[3] = yi - a[3];
    a[n - 2] += yr;
    a[n - 1] = yi - a[n - 1];
    a[1] = -a[1];
}


void dctsub(int n, double *a)
{
    int i, i0, j, k, m;
    double ec, w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi;
    
    ec = M_PI_2 / n;
    wkr = 0.5;
    wki = 0.5;
    w1r = cos(ec);
    w1i = sin(ec);
    wdr = 0.5 * (w1r - w1i);
    wdi = 0.5 * (w1r + w1i);
    ss = 2 * w1i;
    m = n >> 1;
    i = 0;
    for (;;) {
        i0 = i + 2 * DCST_LOOP_DIV;
        if (i0 > m - 2) {
            i0 = m - 2;
        }
        for (j = i + 2; j <= i0; j += 2) {
            k = n - j;
            xr = wdi * a[j - 1] - wdr * a[k + 1];
            xi = wdr * a[j - 1] + wdi * a[k + 1];
            wkr -= ss * wdi;
            wki += ss * wdr;
            yr = wki * a[j] - wkr * a[k];
            yi = wkr * a[j] + wki * a[k];
            wdr -= ss * wki;
            wdi += ss * wkr;
            a[k + 1] = xr;
            a[k] = yr;
            a[j - 1] = xi;
            a[j] = yi;
        }
        if (i0 == m - 2) {
            break;
        }
        wdr = cos(ec * i0);
        wdi = sin(ec * i0);
        wkr = 0.5 * (wdr - wdi);
        wki = 0.5 * (wdr + wdi);
        wdr = wkr * w1r - wki * w1i;
        wdi = wkr * w1i + wki * w1r;
        i = i0;
    }
    xr = wdi * a[m - 1] - wdr * a[m + 1];
    a[m - 1] = wdr * a[m - 1] + wdi * a[m + 1];
    a[m + 1] = xr;
    a[m] *= wki + ss * wdr;
}


void dstsub(int n, double *a)
{
    int i, i0, j, k, m;
    double ec, w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi;
    
    ec = M_PI_2 / n;
    wkr = 0.5;
    wki = 0.5;
    w1r = cos(ec);
    w1i = sin(ec);
    wdr = 0.5 * (w1r - w1i);
    wdi = 0.5 * (w1r + w1i);
    ss = 2 * w1i;
    m = n >> 1;
    i = 0;
    for (;;) {
        i0 = i + 2 * DCST_LOOP_DIV;
        if (i0 > m - 2) {
            i0 = m - 2;
        }
        for (j = i + 2; j <= i0; j += 2) {
            k = n - j;
            xr = wdi * a[k + 1] - wdr * a[j - 1];
            xi = wdr * a[k + 1] + wdi * a[j - 1];
            wkr -= ss * wdi;
            wki += ss * wdr;
            yr = wki * a[k] - wkr * a[j];
            yi = wkr * a[k] + wki * a[j];
            wdr -= ss * wki;
            wdi += ss * wkr;
            a[j - 1] = xr;
            a[j] = yr;
            a[k + 1] = xi;
            a[k] = yi;
        }
        if (i0 == m - 2) {
            break;
        }
        wdr = cos(ec * i0);
        wdi = sin(ec * i0);
        wkr = 0.5 * (wdr - wdi);
        wki = 0.5 * (wdr + wdi);
        wdr = wkr * w1r - wki * w1i;
        wdi = wkr * w1i + wki * w1r;
        i = i0;
    }
    xr = wdi * a[m + 1] - wdr * a[m - 1];
    a[m + 1] = wdr * a[m + 1] + wdi * a[m - 1];
    a[m - 1] = xr;
    a[m] *= wki + ss * wdr;
}


void dctsub4(int n, double *a)
{
    int m;
    double wki, wdr, wdi, xr;
    
    wki = WR5000;
    m = n >> 1;
    if (m == 2) {
        wdr = wki * WI2500;
        wdi = wki * WR2500;
        xr = wdi * a[1] - wdr * a[3];
        a[1] = wdr * a[1] + wdi * a[3];
        a[3] = xr;
    }
    a[m] *= wki;
}


void dstsub4(int n, double *a)
{
    int m;
    double wki, wdr, wdi, xr;
    
    wki = WR5000;
    m = n >> 1;
    if (m == 2) {
        wdr = wki * WI2500;
        wdi = wki * WR2500;
        xr = wdi * a[3] - wdr * a[1];
        a[3] = wdr * a[3] + wdi * a[1];
        a[1] = xr;
    }
    a[m] *= wki;
}