transform8.cc

各种工程计算的库函数

            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);
        }
        k1 = 4;
        for (k = 2 * m; k < n; k += m) {
            k1 += 4;
            wk1r = w[k1];
            wk1i = w[k1 + 1];
            wk2r = w[k1 + 2];
            wk2i = w[k1 + 3];
            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;
            }
        }
    }
}


T8_INLINE void clTransform8::rftfsub(long n, float *a, long nc, float *c)
{
    long j, k, kk, ks, m;
    float wkr, wki, xr, xi, yr, yi;
    
    m = n >> 1;
    ks = 2 * nc / m;
    kk = 0;
    for (j = 2; j < m; j += 2) {
        k = n - j;
        kk += ks;
        wkr = 0.5f - c[nc - kk];
        wki = c[kk];
        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;
    }
}


T8_INLINE void clTransform8::rftbsub(long n, float *a, long nc, float *c)
{
    long j, k, kk, ks, m;
    float wkr, wki, xr, xi, yr, yi;
    
    a[1] = -a[1];
    m = n >> 1;
    ks = 2 * nc / m;
    kk = 0;
    for (j = 2; j < m; j += 2) {
        k = n - j;
        kk += ks;
        wkr = 0.5f - c[nc - kk];
        wki = c[kk];
        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];
    }
    a[m + 1] = -a[m + 1];
}


T8_INLINE void clTransform8::dctsub(long n, float *a, long nc, float *c)
{
    long j, k, kk, ks, m;
    float wkr, wki, xr;
    
    m = n >> 1;
    ks = nc / n;
    kk = 0;
    for (j = 1; j < m; j++) {
        k = n - j;
        kk += ks;
        wkr = c[kk] - c[nc - kk];
        wki = c[kk] + c[nc - kk];
        xr = wki * a[j] - wkr * a[k];
        a[j] = wkr * a[j] + wki * a[k];
        a[k] = xr;
    }
    a[m] *= c[0];
}


T8_INLINE void clTransform8::dstsub(long n, float *a, long nc, float *c)
{
    long j, k, kk, ks, m;
    float wkr, wki, xr;
    
    m = n >> 1;
    ks = nc / n;
    kk = 0;
    for (j = 1; j < m; j++) {
        k = n - j;
        kk += ks;
        wkr = c[kk] - c[nc - kk];
        wki = c[kk] + c[nc - kk];
        xr = wki * a[k] - wkr * a[j];
        a[k] = wkr * a[k] + wki * a[j];
        a[j] = xr;
    }
    a[m] *= c[0];
}


// ----- double precision version starts here


void clTransform8::cdft(long n, long isgn, double *a, long *ip, double *w)
{
    if (n > (ip[0] << 2)) {
        makewt(n >> 2, ip, w);
    }
    if (n > 4) {
        if (isgn >= 0) {
            bitrv2(n, ip + 2, a);
            cftfsub(n, a, w);
        } else {
            bitrv2conj(n, ip + 2, a);
            cftbsub(n, a, w);
        }
    } else if (n == 4) {
        cftfsub(n, a, w);
    }
}


void clTransform8::rdft(long n, long isgn, double *a, long *ip, double *w)
{
    long nw, nc;
    double xi;
    
    nw = ip[0];
    if (n > (nw << 2)) {
        nw = n >> 2;
        makewt(nw, ip, w);
    }
    nc = ip[1];
    if (n > (nc << 2)) {
        nc = n >> 2;
        makect(nc, ip, w + nw);
    }
    if (isgn >= 0) {
        if (n > 4) {
            bitrv2(n, ip + 2, a);
            cftfsub(n, a, w);
            rftfsub(n, a, nc, w + nw);
        } else if (n == 4) {
            cftfsub(n, a, w);
        }
        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, nc, w + nw);
            bitrv2(n, ip + 2, a);
            cftbsub(n, a, w);
        } else if (n == 4) {
            cftfsub(n, a, w);
        }
    }
}


void clTransform8::ddct(long n, long isgn, double *a, long *ip, double *w)
{
    long j, nw, nc;
    double xr;
    
    nw = ip[0];
    if (n > (nw << 2)) {
        nw = n >> 2;
        makewt(nw, ip, w);
    }
    nc = ip[1];
    if (n > nc) {
        nc = n;
        makect(nc, ip, w + nw);
    }
    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, nc, w + nw);
            bitrv2(n, ip + 2, a);
            cftbsub(n, a, w);
        } else if (n == 4) {
            cftfsub(n, a, w);
        }
    }
    dctsub(n, a, nc, w + nw);
    if (isgn >= 0) {
        if (n > 4) {
            bitrv2(n, ip + 2, a);
            cftfsub(n, a, w);
            rftfsub(n, a, nc, w + nw);
        } else if (n == 4) {
            cftfsub(n, a, w);
        }
        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 clTransform8::ddst(long n, long isgn, double *a, long *ip, double *w)
{
    long j, nw, nc;
    double xr;
    
    nw = ip[0];
    if (n > (nw << 2)) {
        nw = n >> 2;
        makewt(nw, ip, w);
    }
    nc = ip[1];
    if (n > nc) {
        nc = n;
        makect(nc, ip, w + nw);
    }
    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, nc, w + nw);
            bitrv2(n, ip + 2, a);
            cftbsub(n, a, w);
        } else if (n == 4) {
            cftfsub(n, a, w);
        }
    }
    dstsub(n, a, nc, w + nw);
    if (isgn >= 0) {
        if (n > 4) {
            bitrv2(n, ip + 2, a);
            cftfsub(n, a, w);
            rftfsub(n, a, nc, w + nw);
        } else if (n == 4) {
            cftfsub(n, a, w);
        }
        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 clTransform8::dfct(long n, double *a, double *t, long *ip, double *w)
{
    long j, k, l, m, mh, nw, nc;
    double xr, xi, yr, yi;
    
    nw = ip[0];
    if (n > (nw << 3)) {
        nw = n >> 3;
        makewt(nw, ip, w);
    }
    nc = ip[1];
    if (n > (nc << 1)) {
        nc = n >> 1;
        makect(nc, ip, w + nw);
    }
    m = n >> 1;
    yi = a[m];
    xi = a[0] + a[n];
    a[0] -= a[n];
    t[0] = xi - yi;
    t[m] = xi + yi;
    if (n > 2) {
        mh = m >> 1;
        for (j = 1; j < mh; j++) {
            k = m - j;
            xr = a[j] - a[n - j];
            xi = a[j] + a[n - j];
            yr = a[k] - a[n - k];
            yi = a[k] + a[n - k];
            a[j] = xr;
            a[k] = yr;
            t[j] = xi - yi;
            t[k] = xi + yi;
        }
        t[mh] = a[mh] + a[n - mh];
        a[mh] -= a[n - mh];
        dctsub(m, a, nc, w + nw);
        if (m > 4) {
            bitrv2(m, ip + 2, a);
            cftfsub(m, a, w);
            rftfsub(m, a, nc, w + nw);
        } else if (m == 4) {
            cftfsub(m, a, w);
        }
        a[n - 1] = a[0] - a[1];
        a[1] = a[0] + a[1];
        for (j = m - 2; j >= 2; j -= 2) {
            a[2 * j + 1] = a[j] + a[j + 1];
            a[2 * j - 1] = a[j] - a[j + 1];
        }
        l = 2;
        m = mh;
        while (m >= 2) {
            dctsub(m, t, nc, w + nw);
            if (m > 4) {
                bitrv2(m, ip + 2, t);
                cftfsub(m, t, w);
                rftfsub(m, t, nc, w + nw);
            } else if (m == 4) {
                cftfsub(m, t, w);
            }
            a[n - l] = t[0] - t[1];
            a[l] = t[0] + t[1];
            k = 0;
            for (j = 2; j < m; j += 2) {
                k += l << 2;