/*Fast Fourier/Cosine/Sine Transform    dimension   :one    data length :power of 2    decimation  :frequency    radix       :8, 4, 2    data        :inplace    table       :not usefunctions    cdft: Complex Discrete Fourier Transform    rdft: Real Discrete Fourier Transform    ddct: Discrete Cosine Transform    ddst: Discrete Sine Transform    dfct: Cosine Transform of RDFT (Real Symmetric DFT)    dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)function prototypes    void cdft(int, int, double *);    void rdft(int, int, double *);    void ddct(int, int, double *);    void ddst(int, int, double *);    void dfct(int, double *);    void dfst(int, double *);-------- Complex DFT (Discrete Fourier Transform) --------    [definition]        <case1>            X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n        <case2>            X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n        (notes: sum_j=0^n-1 is a summation from j=0 to n-1)    [usage]        <case1>            cdft(2*n, 1, a);        <case2>            cdft(2*n, -1, a);    [parameters]        2*n            :data length (int)                        n >= 1, n = power of 2        a[0...2*n-1]   :input/output data (double *)                        input data                            a[2*j] = Re(x[j]),                             a[2*j+1] = Im(x[j]), 0<=j<n                        output data                            a[2*k] = Re(X[k]),                             a[2*k+1] = Im(X[k]), 0<=k<n    [remark]        Inverse of             cdft(2*n, -1, a);        is             cdft(2*n, 1, a);            for (j = 0; j <= 2 * n - 1; j++) {                a[j] *= 1.0 / n;            }        .-------- Real DFT / Inverse of Real DFT --------    [definition]        <case1> RDFT            R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2            I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2        <case2> IRDFT (excluding scale)            a[k] = (R[0] + R[n/2]*cos(pi*k))/2 +                    sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) +                    sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n    [usage]        <case1>            rdft(n, 1, a);        <case2>            rdft(n, -1, a);    [parameters]        n              :data length (int)                        n >= 2, n = power of 2        a[0...n-1]     :input/output data (double *)                        <case1>                            output data                                a[2*k] = R[k], 0<=k<n/2                                a[2*k+1] = I[k], 0<k<n/2                                a[1] = R[n/2]                        <case2>                            input data                                a[2*j] = R[j], 0<=j<n/2                                a[2*j+1] = I[j], 0<j<n/2                                a[1] = R[n/2]    [remark]        Inverse of             rdft(n, 1, a);        is             rdft(n, -1, a);            for (j = 0; j <= n - 1; j++) {                a[j] *= 2.0 / n;            }        .-------- DCT (Discrete Cosine Transform) / Inverse of DCT --------    [definition]        <case1> IDCT (excluding scale)            C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n        <case2> DCT            C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n    [usage]        <case1>            ddct(n, 1, a);        <case2>            ddct(n, -1, a);    [parameters]        n              :data length (int)                        n >= 2, n = power of 2        a[0...n-1]     :input/output data (double *)                        output data                            a[k] = C[k], 0<=k<n    [remark]        Inverse of             ddct(n, -1, a);        is             a[0] *= 0.5;            ddct(n, 1, a);            for (j = 0; j <= n - 1; j++) {                a[j] *= 2.0 / n;            }        .-------- DST (Discrete Sine Transform) / Inverse of DST --------    [definition]        <case1> IDST (excluding scale)            S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n        <case2> DST            S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n    [usage]        <case1>            ddst(n, 1, a);        <case2>            ddst(n, -1, a);    [parameters]        n              :data length (int)                        n >= 2, n = power of 2        a[0...n-1]     :input/output data (double *)                        <case1>                            input data                                a[j] = A[j], 0<j<n                                a[0] = A[n]                            output data                                a[k] = S[k], 0<=k<n                        <case2>                            output data                                a[k] = S[k], 0<k<n                                a[0] = S[n]    [remark]        Inverse of             ddst(n, -1, a);        is             a[0] *= 0.5;            ddst(n, 1, a);            for (j = 0; j <= n - 1; j++) {                a[j] *= 2.0 / n;            }        .-------- Cosine Transform of RDFT (Real Symmetric DFT) --------    [definition]        C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n    [usage]        dfct(n, a);    [parameters]        n              :data length - 1 (int)                        n >= 2, n = power of 2        a[0...n]       :input/output data (double *)                        output data                            a[k] = C[k], 0<=k<=n    [remark]        Inverse of             a[0] *= 0.5;            a[n] *= 0.5;            dfct(n, a);        is             a[0] *= 0.5;            a[n] *= 0.5;            dfct(n, a);            for (j = 0; j <= n; j++) {                a[j] *= 2.0 / n;            }        .-------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------    [definition]        S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n    [usage]        dfst(n, a);    [parameters]        n              :data length + 1 (int)                        n >= 2, n = power of 2        a[0...n-1]     :input/output data (double *)                        output data                            a[k] = S[k], 0<k<n                        (a[0] is used for work area)    [remark]        Inverse of             dfst(n, a);        is             dfst(n, a);            for (j = 1; j <= n - 1; j++) {                a[j] *= 2.0 / n;            }        .*/void cdft(int n, int isgn, double *a){    void bitrv2(int n, double *a);    void bitrv2conj(int n, double *a);    void cftfsub(int n, double *a);    void cftbsub(int n, double *a);        if (n > 4) {        if (isgn >= 0) {            bitrv2(n, a);            cftfsub(n, a);        } else {            bitrv2conj(n, a);            cftbsub(n, a);        }    } else if (n == 4) {        cftfsub(n, a);    }}void rdft(int n, int isgn, double *a){    void bitrv2(int n, double *a);    void cftfsub(int n, double *a);    void cftbsub(int n, double *a);    void rftfsub(int n, double *a);    void rftbsub(int n, double *a);    double xi;        if (isgn >= 0) {        if (n > 4) {            bitrv2(n, a);            cftfsub(n, a);            rftfsub(n, a);        } else if (n == 4) {            cftfsub(n, a);        }        xi = a[0] - a[1];        a[0] += a[1];        a[1] = xi;    } else {        a[1] = 0.5 * (a[0] - a[1]);        a[0] -= a[1];        if (n > 4) {            rftbsub(n, a);            bitrv2(n, a);            cftbsub(n, a);        } else if (n == 4) {            cftfsub(n, a);        }    }}void ddct(int n, int isgn, double *a){    void bitrv2(int n, double *a);    void cftfsub(int n, double *a);    void cftbsub(int n, double *a);    void rftfsub(int n, double *a);    void rftbsub(int n, double *a);    void dctsub(int n, double *a);    void dctsub4(int n, double *a);    int j;    double xr;        if (isgn < 0) {        xr = a[n - 1];        for (j = n - 2; j >= 2; j -= 2) {            a[j + 1] = a[j] - a[j - 1];            a[j] += a[j - 1];        }        a[1] = a[0] - xr;        a[0] += xr;        if (n > 4) {            rftbsub(n, a);            bitrv2(n, a);            cftbsub(n, a);        } else if (n == 4) {            cftfsub(n, a);        }    }    if (n > 4) {        dctsub(n, a);    } else {        dctsub4(n, a);    }    if (isgn >= 0) {        if (n > 4) {            bitrv2(n, a);            cftfsub(n, a);            rftfsub(n, a);        } else if (n == 4) {            cftfsub(n, a);        }        xr = a[0] - a[1];        a[0] += a[1];        for (j = 2; j < n; j += 2) {            a[j - 1] = a[j] - a[j + 1];            a[j] += a[j + 1];        }        a[n - 1] = xr;    }}void ddst(int n, int isgn, double *a){    void bitrv2(int n, double *a);    void cftfsub(int n, double *a);    void cftbsub(int n, double *a);    void rftfsub(int n, double *a);    void rftbsub(int n, double *a);    void dstsub(int n, double *a);    void dstsub4(int n, double *a);    int j;    double xr;        if (isgn < 0) {        xr = a[n - 1];        for (j = n - 2; j >= 2; j -= 2) {            a[j + 1] = -a[j] - a[j - 1];            a[j] -= a[j - 1];        }        a[1] = a[0] + xr;        a[0] -= xr;        if (n > 4) {            rftbsub(n, a);            bitrv2(n, a);            cftbsub(n, a);        } else if (n == 4) {            cftfsub(n, a);        }    }    if (n > 4) {        dstsub(n, a);    } else {        dstsub4(n, a);    }    if (isgn >= 0) {        if (n > 4) {            bitrv2(n, a);            cftfsub(n, a);            rftfsub(n, a);        } else if (n == 4) {            cftfsub(n, a);        }        xr = a[0] - a[1];        a[0] += a[1];        for (j = 2; j < n; j += 2) {            a[j - 1] = -a[j] - a[j + 1];            a[j] -= a[j + 1];        }        a[n - 1] = -xr;    }}void dfct(int n, double *a){    void ddct(int n, int isgn, double *a);    void bitrv1(int n, double *a);    int j, k, m, mh;    double xr, xi, yr, yi, an;        m = n >> 1;    for (j = 0; j < m; j++) {        k = n - j;        xr = a[j] + a[k];        a[j] -= a[k];        a[k] = xr;    }    an = a[n];    while (m >= 2) {        ddct(m, 1, a);        bitrv1(m, a);        mh = m >> 1;        xi = a[m];        a[m] = a[0];        a[0] = an - xi;        an += xi;        for (j = 1; j < mh; j++) {            k = m - j;            xr = a[m + k];            xi = a[m + j];            yr = a[j];            yi = a[k];            a[m + j] = yr;            a[m + k] = yi;            a[j] = xr - xi;            a[k] = xr + xi;        }        xr = a[mh];        a[mh] = a[m + mh];        a[m + mh] = xr;        m = mh;    }    xi = a[1];    a[1] = a[0];    a[0] = an + xi;    a[n] = an - xi;    bitrv1(n, a);}void dfst(int n, double *a){    void ddst(int n, int isgn, double *a);    void bitrv1(int n, double *a);    int j, k, m, mh;    double xr, xi, yr, yi;        m = n >> 1;    for (j = 1; j < m; j++) {        k = n - j;        xr = a[j] - a[k];        a[j] += a[k];        a[k] = xr;    }    a[0] = a[m];    while (m >= 2) {        ddst(m, 1, a);        bitrv1(m, a);        mh = m >> 1;        for (j = 1; j < mh; j++) {            k = m - j;            xr = a[m + k];            xi = a[m + j];            yr = a[j];            yi = a[k];            a[m + j] = yr;            a[m + k] = yi;            a[j] = xr + xi;            a[k] = xr - xi;        }        a[m] = a[0];        a[0] = a[m + mh];        a[m + mh] = a[mh];        m = mh;    }    a[1] = a[0];    a[0] = 0;    bitrv1(n, a);}/* -------- child routines -------- */#include <math.h>#ifndef M_PI_2#define M_PI_2      1.570796326794896619231321691639751442098584699687#endif#ifndef WR5000  /* cos(M_PI_2*0.5000) */#define WR5000      0.707106781186547524400844362104849039284835937688#endif#ifndef WR2500  /* cos(M_PI_2*0.2500) */#define WR2500      0.923879532511286756128183189396788286822416625863#endif#ifndef WI2500  /* sin(M_PI_2*0.2500) */#define WI2500      0.382683432365089771728459984030398866761344562485#endif#ifndef RDFT_LOOP_DIV  /* control of the RDFT's speed & tolerance */#define RDFT_LOOP_DIV 64#endif#ifndef DCST_LOOP_DIV  /* control of the DCT,DST's speed & tolerance */#define DCST_LOOP_DIV 64#endifvoid bitrv2(int n, double *a){    int j0, k0, j1, k1, l, m, i, j, k;    double xr, xi, yr, yi;        l = n >> 2;    m = 2;    while (m < l) {        l >>= 1;        m <<= 1;    }    if (m == l) {        j0 = 0;        for (k0 = 0; k0 < m; k0 += 2) {            k = k0;            for (j = j0; j < j0 + k0; j += 2) {                xr = a[j];                xi = a[j + 1];                yr = a[k];                yi = a[k + 1];                a[j] = yr;                a[j + 1] = yi;                a[k] = xr;                a[k + 1] = xi;                j1 = j + m;                k1 = k + 2 * m;                xr = a[j1];                xi = a[j1 + 1];                yr = a[k1];                yi = a[k1 + 1];                a[j1] = yr;                a[j1 + 1] = yi;                a[k1] = xr;                a[k1 + 1] = xi;                j1 += m;                k1 -= m;                xr = a[j1];                xi = a[j1 + 1];                yr = a[k1];                yi = a[k1 + 1];                a[j1] = yr;                a[j1 + 1] = yi;                a[k1] = xr;                a[k1 + 1] = xi;                j1 += m;                k1 += 2 * m;                xr = a[j1];                xi = a[j1 + 1];                yr = a[k1];                yi = a[k1 + 1];                a[j1] = yr;                a[j1 + 1] = yi;                a[k1] = xr;                a[k1 + 1] = xi;                for (i = n >> 1; i > (k ^= i); i >>= 1);            }            j1 = j0 + k0 + m;            k1 = j1 + m;            xr = a[j1];            xi = a[j1 + 1];            yr = a[k1];            yi = a[k1 + 1];            a[j1] = yr;            a[j1 + 1] = yi;            a[k1] = xr;            a[k1 + 1] = xi;            for (i = n >> 1; i > (j0 ^= i); i >>= 1);        }