| /* |
| Fast Fourier/Cosine/Sine Transform |
| dimension :two |
| data length :power of 2 |
| decimation :frequency |
| radix :4, 2, row-column |
| data :inplace |
| table :use |
| functions |
| cdft2d: Complex Discrete Fourier Transform |
| rdft2d: Real Discrete Fourier Transform |
| ddct2d: Discrete Cosine Transform |
| ddst2d: Discrete Sine Transform |
| function prototypes |
| void cdft2d(int, int, int, double **, int *, double *); |
| void rdft2d(int, int, int, double **, int *, double *); |
| void ddct2d(int, int, int, double **, double **, int *, double *); |
| void ddst2d(int, int, int, double **, double **, int *, double *); |
| |
| |
| -------- Complex DFT (Discrete Fourier Transform) -------- |
| [definition] |
| <case1> |
| X[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 x[j1][j2] * |
| exp(2*pi*i*j1*k1/n1) * |
| exp(2*pi*i*j2*k2/n2), 0<=k1<n1, 0<=k2<n2 |
| <case2> |
| X[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 x[j1][j2] * |
| exp(-2*pi*i*j1*k1/n1) * |
| exp(-2*pi*i*j2*k2/n2), 0<=k1<n1, 0<=k2<n2 |
| (notes: sum_j=0^n-1 is a summation from j=0 to n-1) |
| [usage] |
| <case1> |
| ip[0] = 0; // first time only |
| cdft2d(n1, 2*n2, 1, a, ip, w); |
| <case2> |
| ip[0] = 0; // first time only |
| cdft2d(n1, 2*n2, -1, a, ip, w); |
| [parameters] |
| n1 :data length (int) |
| n1 >= 1, n1 = power of 2 |
| 2*n2 :data length (int) |
| n2 >= 1, n2 = power of 2 |
| a[0...n1-1][0...2*n2-1] |
| :input/output data (double **) |
| input data |
| a[j1][2*j2] = Re(x[j1][j2]), |
| a[j1][2*j2+1] = Im(x[j1][j2]), |
| 0<=j1<n1, 0<=j2<n2 |
| output data |
| a[k1][2*k2] = Re(X[k1][k2]), |
| a[k1][2*k2+1] = Im(X[k1][k2]), |
| 0<=k1<n1, 0<=k2<n2 |
| ip[0...*] |
| :work area for bit reversal (int *) |
| length of ip >= 2+sqrt(n) |
| (n = max(n1, n2)) |
| ip[0],ip[1] are pointers of the cos/sin table. |
| w[0...*] |
| :cos/sin table (double *) |
| length of w >= max(n1/2, n2/2) |
| w[],ip[] are initialized if ip[0] == 0. |
| [remark] |
| Inverse of |
| cdft2d(n1, 2*n2, -1, a, ip, w); |
| is |
| cdft2d(n1, 2*n2, 1, a, ip, w); |
| for (j1 = 0; j1 <= n1 - 1; j1++) { |
| for (j2 = 0; j2 <= 2 * n2 - 1; j2++) { |
| a[j1][j2] *= 1.0 / (n1 * n2); |
| } |
| } |
| . |
| |
| |
| -------- Real DFT / Inverse of Real DFT -------- |
| [definition] |
| <case1> RDFT |
| R[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 a[j1][j2] * |
| cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2), |
| 0<=k1<n1, 0<=k2<n2 |
| I[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 a[j1][j2] * |
| sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2), |
| 0<=k1<n1, 0<=k2<n2 |
| <case2> IRDFT (excluding scale) |
| a[k1][k2] = (1/2) * sum_j1=0^n1-1 sum_j2=0^n2-1 |
| (R[j1][j2] * |
| cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2) + |
| I[j1][j2] * |
| sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2)), |
| 0<=k1<n1, 0<=k2<n2 |
| (notes: R[n1-k1][n2-k2] = R[k1][k2], |
| I[n1-k1][n2-k2] = -I[k1][k2], |
| R[n1-k1][0] = R[k1][0], |
| I[n1-k1][0] = -I[k1][0], |
| R[0][n2-k2] = R[0][k2], |
| I[0][n2-k2] = -I[0][k2], |
| 0<k1<n1, 0<k2<n2) |
| [usage] |
| <case1> |
| ip[0] = 0; // first time only |
| rdft2d(n1, n2, 1, a, ip, w); |
| <case2> |
| ip[0] = 0; // first time only |
| rdft2d(n1, n2, -1, a, ip, w); |
| [parameters] |
| n1 :data length (int) |
| n1 >= 2, n1 = power of 2 |
| n2 :data length (int) |
| n2 >= 2, n2 = power of 2 |
| a[0...n1-1][0...n2-1] |
| :input/output data (double **) |
| <case1> |
| output data |
| a[k1][2*k2] = R[k1][k2] = R[n1-k1][n2-k2], |
| a[k1][2*k2+1] = I[k1][k2] = -I[n1-k1][n2-k2], |
| 0<k1<n1, 0<k2<n2/2, |
| a[0][2*k2] = R[0][k2] = R[0][n2-k2], |
| a[0][2*k2+1] = I[0][k2] = -I[0][n2-k2], |
| 0<k2<n2/2, |
| a[k1][0] = R[k1][0] = R[n1-k1][0], |
| a[k1][1] = I[k1][0] = -I[n1-k1][0], |
| a[n1-k1][1] = R[k1][n2/2] = R[n1-k1][n2/2], |
| a[n1-k1][0] = -I[k1][n2/2] = I[n1-k1][n2/2], |
| 0<k1<n1/2, |
| a[0][0] = R[0][0], |
| a[0][1] = R[0][n2/2], |
| a[n1/2][0] = R[n1/2][0], |
| a[n1/2][1] = R[n1/2][n2/2] |
| <case2> |
| input data |
| a[j1][2*j2] = R[j1][j2] = R[n1-j1][n2-j2], |
| a[j1][2*j2+1] = I[j1][j2] = -I[n1-j1][n2-j2], |
| 0<j1<n1, 0<j2<n2/2, |
| a[0][2*j2] = R[0][j2] = R[0][n2-j2], |
| a[0][2*j2+1] = I[0][j2] = -I[0][n2-j2], |
| 0<j2<n2/2, |
| a[j1][0] = R[j1][0] = R[n1-j1][0], |
| a[j1][1] = I[j1][0] = -I[n1-j1][0], |
| a[n1-j1][1] = R[j1][n2/2] = R[n1-j1][n2/2], |
| a[n1-j1][0] = -I[j1][n2/2] = I[n1-j1][n2/2], |
| 0<j1<n1/2, |
| a[0][0] = R[0][0], |
| a[0][1] = R[0][n2/2], |
| a[n1/2][0] = R[n1/2][0], |
| a[n1/2][1] = R[n1/2][n2/2] |
| ip[0...*] |
| :work area for bit reversal (int *) |
| length of ip >= 2+sqrt(n) |
| (n = max(n1, n2/2)) |
| ip[0],ip[1] are pointers of the cos/sin table. |
| w[0...*] |
| :cos/sin table (double *) |
| length of w >= max(n1/2, n2/4) + n2/4 |
| w[],ip[] are initialized if ip[0] == 0. |
| [remark] |
| Inverse of |
| rdft2d(n1, n2, 1, a, ip, w); |
| is |
| rdft2d(n1, n2, -1, a, ip, w); |
| for (j1 = 0; j1 <= n1 - 1; j1++) { |
| for (j2 = 0; j2 <= n2 - 1; j2++) { |
| a[j1][j2] *= 2.0 / (n1 * n2); |
| } |
| } |
| . |
| |
| |
| -------- DCT (Discrete Cosine Transform) / Inverse of DCT -------- |
| [definition] |
| <case1> IDCT (excluding scale) |
| C[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 a[j1][j2] * |
| cos(pi*j1*(k1+1/2)/n1) * |
| cos(pi*j2*(k2+1/2)/n2), |
| 0<=k1<n1, 0<=k2<n2 |
| <case2> DCT |
| C[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 a[j1][j2] * |
| cos(pi*(j1+1/2)*k1/n1) * |
| cos(pi*(j2+1/2)*k2/n2), |
| 0<=k1<n1, 0<=k2<n2 |
| [usage] |
| <case1> |
| ip[0] = 0; // first time only |
| ddct2d(n1, n2, 1, a, t, ip, w); |
| <case2> |
| ip[0] = 0; // first time only |
| ddct2d(n1, n2, -1, a, t, ip, w); |
| [parameters] |
| n1 :data length (int) |
| n1 >= 2, n1 = power of 2 |
| n2 :data length (int) |
| n2 >= 2, n2 = power of 2 |
| a[0...n1-1][0...n2-1] |
| :input/output data (double **) |
| output data |
| a[k1][k2] = C[k1][k2], 0<=k1<n1, 0<=k2<n2 |
| t[0...n1-1][0...n2-1] |
| :work area (double **) |
| ip[0...*] |
| :work area for bit reversal (int *) |
| length of ip >= 2+sqrt(n) |
| (n = max(n1, n2/2)) |
| ip[0],ip[1] are pointers of the cos/sin table. |
| w[0...*] |
| :cos/sin table (double *) |
| length of w >= max(n1/2, n2/4) + max(n1, n2) |
| w[],ip[] are initialized if ip[0] == 0. |
| [remark] |
| Inverse of |
| ddct2d(n1, n2, -1, a, t, ip, w); |
| is |
| for (j1 = 0; j1 <= n1 - 1; j1++) { |
| a[j1][0] *= 0.5; |
| } |
| for (j2 = 0; j2 <= n2 - 1; j2++) { |
| a[0][j2] *= 0.5; |
| } |
| ddct2d(n1, n2, 1, a, t, ip, w); |
| for (j1 = 0; j1 <= n1 - 1; j1++) { |
| for (j2 = 0; j2 <= n2 - 1; j2++) { |
| a[j1][j2] *= 4.0 / (n1 * n2); |
| } |
| } |
| . |
| |
| |
| -------- DST (Discrete Sine Transform) / Inverse of DST -------- |
| [definition] |
| <case1> IDST (excluding scale) |
| S[k1][k2] = sum_j1=1^n1 sum_j2=1^n2 A[j1][j2] * |
| sin(pi*j1*(k1+1/2)/n1) * |
| sin(pi*j2*(k2+1/2)/n2), |
| 0<=k1<n1, 0<=k2<n2 |
| <case2> DST |
| S[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 a[j1][j2] * |
| sin(pi*(j1+1/2)*k1/n1) * |
| sin(pi*(j2+1/2)*k2/n2), |
| 0<k1<=n1, 0<k2<=n2 |
| [usage] |
| <case1> |
| ip[0] = 0; // first time only |
| ddst2d(n1, n2, 1, a, t, ip, w); |
| <case2> |
| ip[0] = 0; // first time only |
| ddst2d(n1, n2, -1, a, t, ip, w); |
| [parameters] |
| n1 :data length (int) |
| n1 >= 2, n1 = power of 2 |
| n2 :data length (int) |
| n2 >= 2, n2 = power of 2 |
| a[0...n1-1][0...n2-1] |
| :input/output data (double **) |
| <case1> |
| input data |
| a[j1][j2] = A[j1][j2], 0<j1<n1, 0<j2<n2, |
| a[j1][0] = A[j1][n2], 0<j1<n1, |
| a[0][j2] = A[n1][j2], 0<j2<n2, |
| a[0][0] = A[n1][n2] |
| (i.e. A[j1][j2] = a[j1 % n1][j2 % n2]) |
| output data |
| a[k1][k2] = S[k1][k2], 0<=k1<n1, 0<=k2<n2 |
| <case2> |
| output data |
| a[k1][k2] = S[k1][k2], 0<k1<n1, 0<k2<n2, |
| a[k1][0] = S[k1][n2], 0<k1<n1, |
| a[0][k2] = S[n1][k2], 0<k2<n2, |
| a[0][0] = S[n1][n2] |
| (i.e. S[k1][k2] = a[k1 % n1][k2 % n2]) |
| t[0...n1-1][0...n2-1] |
| :work area (double **) |
| ip[0...*] |
| :work area for bit reversal (int *) |
| length of ip >= 2+sqrt(n) |
| (n = max(n1, n2/2)) |
| ip[0],ip[1] are pointers of the cos/sin table. |
| w[0...*] |
| :cos/sin table (double *) |
| length of w >= max(n1/2, n2/4) + max(n1, n2) |
| w[],ip[] are initialized if ip[0] == 0. |
| [remark] |
| Inverse of |
| ddst2d(n1, n2, -1, a, t, ip, w); |
| is |
| for (j1 = 0; j1 <= n1 - 1; j1++) { |
| a[j1][0] *= 0.5; |
| } |
| for (j2 = 0; j2 <= n2 - 1; j2++) { |
| a[0][j2] *= 0.5; |
| } |
| ddst2d(n1, n2, 1, a, t, ip, w); |
| for (j1 = 0; j1 <= n1 - 1; j1++) { |
| for (j2 = 0; j2 <= n2 - 1; j2++) { |
| a[j1][j2] *= 4.0 / (n1 * n2); |
| } |
| } |
| . |
| */ |
| |
| |
| void cdft2d(int n1, int n2, int isgn, double **a, int *ip, double *w) |
| { |
| void makewt(int nw, int *ip, double *w); |
| void bitrv2col(int n1, int n, int *ip, double **a); |
| void bitrv2row(int n, int n2, int *ip, double **a); |
| void cftbcol(int n1, int n, double **a, double *w); |
| void cftbrow(int n, int n2, double **a, double *w); |
| void cftfcol(int n1, int n, double **a, double *w); |
| void cftfrow(int n, int n2, double **a, double *w); |
| int n; |
| |
| n = n1 << 1; |
| if (n < n2) { |
| n = n2; |
| } |
| if (n > (ip[0] << 2)) { |
| makewt(n >> 2, ip, w); |
| } |
| if (n2 > 4) { |
| bitrv2col(n1, n2, ip + 2, a); |
| } |
| if (n1 > 2) { |
| bitrv2row(n1, n2, ip + 2, a); |
| } |
| if (isgn < 0) { |
| cftfcol(n1, n2, a, w); |
| cftfrow(n1, n2, a, w); |
| } else { |
| cftbcol(n1, n2, a, w); |
| cftbrow(n1, n2, a, w); |
| } |
| } |
| |
| |
| void rdft2d(int n1, int n2, int isgn, double **a, int *ip, double *w) |
| { |
| void makewt(int nw, int *ip, double *w); |
| void makect(int nc, int *ip, double *c); |
| void bitrv2col(int n1, int n, int *ip, double **a); |
| void bitrv2row(int n, int n2, int *ip, double **a); |
| void cftbcol(int n1, int n, double **a, double *w); |
| void cftbrow(int n, int n2, double **a, double *w); |
| void cftfcol(int n1, int n, double **a, double *w); |
| void cftfrow(int n, int n2, double **a, double *w); |
| void rftbcol(int n1, int n, double **a, int nc, double *c); |
| void rftfcol(int n1, int n, double **a, int nc, double *c); |
| int n, nw, nc, n1h, i, j; |
| double xi; |
| |
| n = n1 << 1; |
| if (n < n2) { |
| n = n2; |
| } |
| nw = ip[0]; |
| if (n > (nw << 2)) { |
| nw = n >> 2; |
| makewt(nw, ip, w); |
| } |
| nc = ip[1]; |
| if (n2 > (nc << 2)) { |
| nc = n2 >> 2; |
| makect(nc, ip, w + nw); |
| } |
| n1h = n1 >> 1; |
| if (isgn < 0) { |
| for (i = 1; i <= n1h - 1; i++) { |
| j = n1 - i; |
| xi = a[i][0] - a[j][0]; |
| a[i][0] += a[j][0]; |
| a[j][0] = xi; |
| xi = a[j][1] - a[i][1]; |
| a[i][1] += a[j][1]; |
| a[j][1] = xi; |
| } |
| if (n1 > 2) { |
| bitrv2row(n1, n2, ip + 2, a); |
| } |
| cftfrow(n1, n2, a, w); |
| for (i = 0; i <= n1 - 1; i++) { |
| a[i][1] = 0.5 * (a[i][0] - a[i][1]); |
| a[i][0] -= a[i][1]; |
| } |
| if (n2 > 4) { |
| rftfcol(n1, n2, a, nc, w + nw); |
| bitrv2col(n1, n2, ip + 2, a); |
| } |
| cftfcol(n1, n2, a, w); |
| } else { |
| if (n2 > 4) { |
| bitrv2col(n1, n2, ip + 2, a); |
| } |
| cftbcol(n1, n2, a, w); |
| if (n2 > 4) { |
| rftbcol(n1, n2, a, nc, w + nw); |
| } |
| for (i = 0; i <= n1 - 1; i++) { |
| xi = a[i][0] - a[i][1]; |
| a[i][0] += a[i][1]; |
| a[i][1] = xi; |
| } |
| if (n1 > 2) { |
| bitrv2row(n1, n2, ip + 2, a); |
| } |
| cftbrow(n1, n2, a, w); |
| for (i = 1; i <= n1h - 1; i++) { |
| j = n1 - i; |
| a[j][0] = 0.5 * (a[i][0] - a[j][0]); |
| a[i][0] -= a[j][0]; |
| a[j][1] = 0.5 * (a[i][1] + a[j][1]); |
| a[i][1] -= a[j][1]; |
| } |
| } |
| } |
| |
| |
| void ddct2d(int n1, int n2, int isgn, double **a, double **t, |
| int *ip, double *w) |
| { |
| void makewt(int nw, int *ip, double *w); |
| void makect(int nc, int *ip, double *c); |
| void bitrv2col(int n1, int n, int *ip, double **a); |
| void bitrv2row(int n, int n2, int *ip, double **a); |
| void cftbcol(int n1, int n, double **a, double *w); |
| void cftbrow(int n, int n2, double **a, double *w); |
| void cftfcol(int n1, int n, double **a, double *w); |
| void cftfrow(int n, int n2, double **a, double *w); |
| void rftbcol(int n1, int n, double **a, int nc, double *c); |
| void rftfcol(int n1, int n, double **a, int nc, double *c); |
| void dctbsub(int n1, int n2, double **a, int nc, double *c); |
| void dctfsub(int n1, int n2, double **a, int nc, double *c); |
| int n, nw, nc, n1h, n2h, i, ix, ic, j, jx, jc; |
| double xi; |
| |
| n = n1 << 1; |
| if (n < n2) { |
| n = n2; |
| } |
| nw = ip[0]; |
| if (n > (nw << 2)) { |
| nw = n >> 2; |
| makewt(nw, ip, w); |
| } |
| nc = ip[1]; |
| if (n1 > nc || n2 > nc) { |
| if (n1 > n2) { |
| nc = n1; |
| } else { |
| nc = n2; |
| } |
| makect(nc, ip, w + nw); |
| } |
| n1h = n1 >> 1; |
| n2h = n2 >> 1; |
| if (isgn >= 0) { |
| for (i = 0; i <= n1 - 1; i++) { |
| for (j = 1; j <= n2h - 1; j++) { |
| jx = j << 1; |
| t[i][jx] = a[i][j]; |
| t[i][jx + 1] = a[i][n2 - j]; |
| } |
| } |
| t[0][0] = a[0][0]; |
| t[0][1] = a[0][n2h]; |
| t[n1h][0] = a[n1h][0]; |
| t[n1h][1] = a[n1h][n2h]; |
| for (i = 1; i <= n1h - 1; i++) { |
| ic = n1 - i; |
| t[i][0] = a[i][0]; |
| t[ic][1] = a[i][n2h]; |
| t[i][1] = a[ic][0]; |
| t[ic][0] = a[ic][n2h]; |
| } |
| dctfsub(n1, n2, t, nc, w + nw); |
| if (n1 > 2) { |
| bitrv2row(n1, n2, ip + 2, t); |
| } |
| cftfrow(n1, n2, t, w); |
| for (i = 0; i <= n1 - 1; i++) { |
| t[i][1] = 0.5 * (t[i][0] - t[i][1]); |
| t[i][0] -= t[i][1]; |
| } |
| if (n2 > 4) { |
| rftfcol(n1, n2, t, nc, w + nw); |
| bitrv2col(n1, n2, ip + 2, t); |
| } |
| cftfcol(n1, n2, t, w); |
| for (i = 0; i <= n1h - 1; i++) { |
| ix = i << 1; |
| ic = n1 - 1 - i; |
| for (j = 0; j <= n2h - 1; j++) { |
| jx = j << 1; |
| jc = n2 - 1 - j; |
| a[ix][jx] = t[i][j]; |
| a[ix][jx + 1] = t[i][jc]; |
| a[ix + 1][jx] = t[ic][j]; |
| a[ix + 1][jx + 1] = t[ic][jc]; |
| } |
| } |
| } else { |
| for (i = 0; i <= n1h - 1; i++) { |
| ix = i << 1; |
| ic = n1 - 1 - i; |
| for (j = 0; j <= n2h - 1; j++) { |
| jx = j << 1; |
| jc = n2 - 1 - j; |
| t[i][j] = a[ix][jx]; |
| t[i][jc] = a[ix][jx + 1]; |
| t[ic][j] = a[ix + 1][jx]; |
| t[ic][jc] = a[ix + 1][jx + 1]; |
| } |
| } |
| if (n2 > 4) { |
| bitrv2col(n1, n2, ip + 2, t); |
| } |
| cftbcol(n1, n2, t, w); |
| if (n2 > 4) { |
| rftbcol(n1, n2, t, nc, w + nw); |
| } |
| for (i = 0; i <= n1 - 1; i++) { |
| xi = t[i][0] - t[i][1]; |
| t[i][0] += t[i][1]; |
| t[i][1] = xi; |
| } |
| if (n1 > 2) { |
| bitrv2row(n1, n2, ip + 2, t); |
| } |
| cftbrow(n1, n2, t, w); |
| dctbsub(n1, n2, t, nc, w + nw); |
| for (i = 0; i <= n1 - 1; i++) { |
| for (j = 1; j <= n2h - 1; j++) { |
| jx = j << 1; |
| a[i][j] = t[i][jx]; |
| a[i][n2 - j] = t[i][jx + 1]; |
| } |
| } |
| a[0][0] = t[0][0]; |
| a[0][n2h] = t[0][1]; |
| a[n1h][0] = t[n1h][0]; |
| a[n1h][n2h] = t[n1h][1]; |
| for (i = 1; i <= n1h - 1; i++) { |
| ic = n1 - i; |
| a[i][0] = t[i][0]; |
| a[i][n2h] = t[ic][1]; |
| a[ic][0] = t[i][1]; |
| a[ic][n2h] = t[ic][0]; |
| } |
| } |
| } |
| |
| |
| void ddst2d(int n1, int n2, int isgn, double **a, double **t, |
| int *ip, double *w) |
| { |
| void makewt(int nw, int *ip, double *w); |
| void makect(int nc, int *ip, double *c); |
| void bitrv2col(int n1, int n, int *ip, double **a); |
| void bitrv2row(int n, int n2, int *ip, double **a); |
| void cftbcol(int n1, int n, double **a, double *w); |
| void cftbrow(int n, int n2, double **a, double *w); |
| void cftfcol(int n1, int n, double **a, double *w); |
| void cftfrow(int n, int n2, double **a, double *w); |
| void rftbcol(int n1, int n, double **a, int nc, double *c); |
| void rftfcol(int n1, int n, double **a, int nc, double *c); |
| void dstbsub(int n1, int n2, double **a, int nc, double *c); |
| void dstfsub(int n1, int n2, double **a, int nc, double *c); |
| int n, nw, nc, n1h, n2h, i, ix, ic, j, jx, jc; |
| double xi; |
| |
| n = n1 << 1; |
| if (n < n2) { |
| n = n2; |
| } |
| nw = ip[0]; |
| if (n > (nw << 2)) { |
| nw = n >> 2; |
| makewt(nw, ip, w); |
| } |
| nc = ip[1]; |
| if (n1 > nc || n2 > nc) { |
| if (n1 > n2) { |
| nc = n1; |
| } else { |
| nc = n2; |
| } |
| makect(nc, ip, w + nw); |
| } |
| n1h = n1 >> 1; |
| n2h = n2 >> 1; |
| if (isgn >= 0) { |
| for (i = 0; i <= n1 - 1; i++) { |
| for (j = 1; j <= n2h - 1; j++) { |
| jx = j << 1; |
| t[i][jx] = a[i][j]; |
| t[i][jx + 1] = a[i][n2 - j]; |
| } |
| } |
| t[0][0] = a[0][0]; |
| t[0][1] = a[0][n2h]; |
| t[n1h][0] = a[n1h][0]; |
| t[n1h][1] = a[n1h][n2h]; |
| for (i = 1; i <= n1h - 1; i++) { |
| ic = n1 - i; |
| t[i][0] = a[i][0]; |
| t[ic][1] = a[i][n2h]; |
| t[i][1] = a[ic][0]; |
| t[ic][0] = a[ic][n2h]; |
| } |
| dstfsub(n1, n2, t, nc, w + nw); |
| if (n1 > 2) { |
| bitrv2row(n1, n2, ip + 2, t); |
| } |
| cftfrow(n1, n2, t, w); |
| for (i = 0; i <= n1 - 1; i++) { |
| t[i][1] = 0.5 * (t[i][0] - t[i][1]); |
| t[i][0] -= t[i][1]; |
| } |
| if (n2 > 4) { |
| rftfcol(n1, n2, t, nc, w + nw); |
| bitrv2col(n1, n2, ip + 2, t); |
| } |
| cftfcol(n1, n2, t, w); |
| for (i = 0; i <= n1h - 1; i++) { |
| ix = i << 1; |
| ic = n1 - 1 - i; |
| for (j = 0; j <= n2h - 1; j++) { |
| jx = j << 1; |
| jc = n2 - 1 - j; |
| a[ix][jx] = t[i][j]; |
| a[ix][jx + 1] = -t[i][jc]; |
| a[ix + 1][jx] = -t[ic][j]; |
| a[ix + 1][jx + 1] = t[ic][jc]; |
| } |
| } |
| } else { |
| for (i = 0; i <= n1h - 1; i++) { |
| ix = i << 1; |
| ic = n1 - 1 - i; |
| for (j = 0; j <= n2h - 1; j++) { |
| jx = j << 1; |
| jc = n2 - 1 - j; |
| t[i][j] = a[ix][jx]; |
| t[i][jc] = -a[ix][jx + 1]; |
| t[ic][j] = -a[ix + 1][jx]; |
| t[ic][jc] = a[ix + 1][jx + 1]; |
| } |
| } |
| if (n2 > 4) { |
| bitrv2col(n1, n2, ip + 2, t); |
| } |
| cftbcol(n1, n2, t, w); |
| if (n2 > 4) { |
| rftbcol(n1, n2, t, nc, w + nw); |
| } |
| for (i = 0; i <= n1 - 1; i++) { |
| xi = t[i][0] - t[i][1]; |
| t[i][0] += t[i][1]; |
| t[i][1] = xi; |
| } |
| if (n1 > 2) { |
| bitrv2row(n1, n2, ip + 2, t); |
| } |
| cftbrow(n1, n2, t, w); |
| dstbsub(n1, n2, t, nc, w + nw); |
| for (i = 0; i <= n1 - 1; i++) { |
| for (j = 1; j <= n2h - 1; j++) { |
| jx = j << 1; |
| a[i][j] = t[i][jx]; |
| a[i][n2 - j] = t[i][jx + 1]; |
| } |
| } |
| a[0][0] = t[0][0]; |
| a[0][n2h] = t[0][1]; |
| a[n1h][0] = t[n1h][0]; |
| a[n1h][n2h] = t[n1h][1]; |
| for (i = 1; i <= n1h - 1; i++) { |
| ic = n1 - i; |
| a[i][0] = t[i][0]; |
| a[i][n2h] = t[ic][1]; |
| a[ic][0] = t[i][1]; |
| a[ic][n2h] = t[ic][0]; |
| } |
| } |
| } |
| |
| |
| /* -------- initializing routines -------- */ |
| |
| |
| #include <math.h> |
| |
| void makewt(int nw, int *ip, double *w) |
| { |
| void bitrv2(int n, int *ip, double *a); |
| int nwh, j; |
| double delta, x, y; |
| |
| ip[0] = nw; |
| ip[1] = 1; |
| if (nw > 2) { |
| nwh = nw >> 1; |
| delta = atan(1.0) / nwh; |
| w[0] = 1; |
| w[1] = 0; |
| w[nwh] = cos(delta * nwh); |
| w[nwh + 1] = w[nwh]; |
| for (j = 2; j <= nwh - 2; j += 2) { |
| x = cos(delta * j); |
| y = sin(delta * j); |
| w[j] = x; |
| w[j + 1] = y; |
| w[nw - j] = y; |
| w[nw - j + 1] = x; |
| } |
| bitrv2(nw, ip + 2, w); |
| } |
| } |
| |
| |
| void makect(int nc, int *ip, double *c) |
| { |
| int nch, j; |
| double delta; |
| |
| ip[1] = nc; |
| if (nc > 1) { |
| nch = nc >> 1; |
| delta = atan(1.0) / nch; |
| c[0] = 0.5; |
| c[nch] = 0.5 * cos(delta * nch); |
| for (j = 1; j <= nch - 1; j++) { |
| c[j] = 0.5 * cos(delta * j); |
| c[nc - j] = 0.5 * sin(delta * j); |
| } |
| } |
| } |
| |
| |
| /* -------- child routines -------- */ |
| |
| |
| void bitrv2(int n, int *ip, double *a) |
| { |
| int j, j1, k, k1, l, m, m2; |
| double xr, xi; |
| |
| ip[0] = 0; |
| l = n; |
| m = 1; |
| while ((m << 2) < l) { |
| l >>= 1; |
| for (j = 0; j <= m - 1; j++) { |
| ip[m + j] = ip[j] + l; |
| } |
| m <<= 1; |
| } |
| if ((m << 2) > l) { |
| for (k = 1; k <= m - 1; k++) { |
| for (j = 0; j <= k - 1; j++) { |
| j1 = (j << 1) + ip[k]; |
| k1 = (k << 1) + ip[j]; |
| xr = a[j1]; |
| xi = a[j1 + 1]; |
| a[j1] = a[k1]; |
| a[j1 + 1] = a[k1 + 1]; |
| a[k1] = xr; |
| a[k1 + 1] = xi; |
| } |
| } |
| } else { |
| m2 = m << 1; |
| for (k = 1; k <= m - 1; k++) { |
| for (j = 0; j <= k - 1; j++) { |
| j1 = (j << 1) + ip[k]; |
| k1 = (k << 1) + ip[j]; |
| xr = a[j1]; |
| xi = a[j1 + 1]; |
| a[j1] = a[k1]; |
| a[j1 + 1] = a[k1 + 1]; |
| a[k1] = xr; |
| a[k1 + 1] = xi; |
| j1 += m2; |
| k1 += m2; |
| xr = a[j1]; |
| xi = a[j1 + 1]; |
| a[j1] = a[k1]; |
| a[j1 + 1] = a[k1 + 1]; |
| a[k1] = xr; |
| a[k1 + 1] = xi; |
| } |
| } |
| } |
| } |
| |
| |
| void bitrv2col(int n1, int n, int *ip, double **a) |
| { |
| int i, j, j1, k, k1, l, m, m2; |
| double xr, xi; |
| |
| ip[0] = 0; |
| l = n; |
| m = 1; |
| while ((m << 2) < l) { |
| l >>= 1; |
| for (j = 0; j <= m - 1; j++) { |
| ip[m + j] = ip[j] + l; |
| } |
| m <<= 1; |
| } |
| if ((m << 2) > l) { |
| for (i = 0; i <= n1 - 1; i++) { |
| for (k = 1; k <= m - 1; k++) { |
| for (j = 0; j <= k - 1; j++) { |
| j1 = (j << 1) + ip[k]; |
| k1 = (k << 1) + ip[j]; |
| xr = a[i][j1]; |
| xi = a[i][j1 + 1]; |
| a[i][j1] = a[i][k1]; |
| a[i][j1 + 1] = a[i][k1 + 1]; |
| a[i][k1] = xr; |
| a[i][k1 + 1] = xi; |
| } |
| } |
| } |
| } else { |
| m2 = m << 1; |
| for (i = 0; i <= n1 - 1; i++) { |
| for (k = 1; k <= m - 1; k++) { |
| for (j = 0; j <= k - 1; j++) { |
| j1 = (j << 1) + ip[k]; |
| k1 = (k << 1) + ip[j]; |
| xr = a[i][j1]; |
| xi = a[i][j1 + 1]; |
| a[i][j1] = a[i][k1]; |
| a[i][j1 + 1] = a[i][k1 + 1]; |
| a[i][k1] = xr; |
| a[i][k1 + 1] = xi; |
| j1 += m2; |
| k1 += m2; |
| xr = a[i][j1]; |
| xi = a[i][j1 + 1]; |
| a[i][j1] = a[i][k1]; |
| a[i][j1 + 1] = a[i][k1 + 1]; |
| a[i][k1] = xr; |
| a[i][k1 + 1] = xi; |
| } |
| } |
| } |
| } |
| } |
| |
| |
| void bitrv2row(int n, int n2, int *ip, double **a) |
| { |
| int i, j, j1, k, k1, l, m; |
| double xr, xi; |
| |
| ip[0] = 0; |
| l = n; |
| m = 1; |
| while ((m << 1) < l) { |
| l >>= 1; |
| for (j = 0; j <= m - 1; j++) { |
| ip[m + j] = ip[j] + l; |
| } |
| m <<= 1; |
| } |
| if ((m << 1) > l) { |
| for (k = 1; k <= m - 1; k++) { |
| for (j = 0; j <= k - 1; j++) { |
| j1 = j + ip[k]; |
| k1 = k + ip[j]; |
| for (i = 0; i <= n2 - 2; i += 2) { |
| xr = a[j1][i]; |
| xi = a[j1][i + 1]; |
| a[j1][i] = a[k1][i]; |
| a[j1][i + 1] = a[k1][i + 1]; |
| a[k1][i] = xr; |
| a[k1][i + 1] = xi; |
| } |
| } |
| } |
| } else { |
| for (k = 1; k <= m - 1; k++) { |
| for (j = 0; j <= k - 1; j++) { |
| j1 = j + ip[k]; |
| k1 = k + ip[j]; |
| for (i = 0; i <= n2 - 2; i += 2) { |
| xr = a[j1][i]; |
| xi = a[j1][i + 1]; |
| a[j1][i] = a[k1][i]; |
| a[j1][i + 1] = a[k1][i + 1]; |
| a[k1][i] = xr; |
| a[k1][i + 1] = xi; |
| } |
| j1 += m; |
| k1 += m; |
| for (i = 0; i <= n2 - 2; i += 2) { |
| xr = a[j1][i]; |
| xi = a[j1][i + 1]; |
| a[j1][i] = a[k1][i]; |
| a[j1][i + 1] = a[k1][i + 1]; |
| a[k1][i] = xr; |
| a[k1][i + 1] = xi; |
| } |
| } |
| } |
| } |
| } |
| |
| |
| void cftbcol(int n1, int n, double **a, double *w) |
| { |
| int i, j, j1, j2, j3, k, k1, ks, l, m; |
| double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; |
| double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; |
| |
| for (i = 0; i <= n1 - 1; i++) { |
| l = 2; |
| while ((l << 1) < n) { |
| m = l << 2; |
| for (j = 0; j <= l - 2; j += 2) { |
| j1 = j + l; |
| j2 = j1 + l; |
| j3 = j2 + l; |
| x0r = a[i][j] + a[i][j1]; |
| x0i = a[i][j + 1] + a[i][j1 + 1]; |
| x1r = a[i][j] - a[i][j1]; |
| x1i = a[i][j + 1] - a[i][j1 + 1]; |
| x2r = a[i][j2] + a[i][j3]; |
| x2i = a[i][j2 + 1] + a[i][j3 + 1]; |
| x3r = a[i][j2] - a[i][j3]; |
| x3i = a[i][j2 + 1] - a[i][j3 + 1]; |
| a[i][j] = x0r + x2r; |
| a[i][j + 1] = x0i + x2i; |
| a[i][j2] = x0r - x2r; |
| a[i][j2 + 1] = x0i - x2i; |
| a[i][j1] = x1r - x3i; |
| a[i][j1 + 1] = x1i + x3r; |
| a[i][j3] = x1r + x3i; |
| a[i][j3 + 1] = x1i - x3r; |
| } |
| if (m < n) { |
| wk1r = w[2]; |
| for (j = m; j <= l + m - 2; j += 2) { |
| j1 = j + l; |
| j2 = j1 + l; |
| j3 = j2 + l; |
| x0r = a[i][j] + a[i][j1]; |
| x0i = a[i][j + 1] + a[i][j1 + 1]; |
| x1r = a[i][j] - a[i][j1]; |
| x1i = a[i][j + 1] - a[i][j1 + 1]; |
| x2r = a[i][j2] + a[i][j3]; |
| x2i = a[i][j2 + 1] + a[i][j3 + 1]; |
| x3r = a[i][j2] - a[i][j3]; |
| x3i = a[i][j2 + 1] - a[i][j3 + 1]; |
| a[i][j] = x0r + x2r; |
| a[i][j + 1] = x0i + x2i; |
| a[i][j2] = x2i - x0i; |
| a[i][j2 + 1] = x0r - x2r; |
| x0r = x1r - x3i; |
| x0i = x1i + x3r; |
| a[i][j1] = wk1r * (x0r - x0i); |
| a[i][j1 + 1] = wk1r * (x0r + x0i); |
| x0r = x3i + x1r; |
| x0i = x3r - x1i; |
| a[i][j3] = wk1r * (x0i - x0r); |
| a[i][j3 + 1] = wk1r * (x0i + x0r); |
| } |
| k1 = 1; |
| ks = -1; |
| for (k = (m << 1); k <= n - m; k += m) { |
| k1++; |
| ks = -ks; |
| wk1r = w[k1 << 1]; |
| wk1i = w[(k1 << 1) + 1]; |
| wk2r = ks * w[k1]; |
| wk2i = w[k1 + ks]; |
| wk3r = wk1r - 2 * wk2i * wk1i; |
| wk3i = 2 * wk2i * wk1r - wk1i; |
| for (j = k; j <= l + k - 2; j += 2) { |
| j1 = j + l; |
| j2 = j1 + l; |
| j3 = j2 + l; |
| x0r = a[i][j] + a[i][j1]; |
| x0i = a[i][j + 1] + a[i][j1 + 1]; |
| x1r = a[i][j] - a[i][j1]; |
| x1i = a[i][j + 1] - a[i][j1 + 1]; |
| x2r = a[i][j2] + a[i][j3]; |
| x2i = a[i][j2 + 1] + a[i][j3 + 1]; |
| x3r = a[i][j2] - a[i][j3]; |
| x3i = a[i][j2 + 1] - a[i][j3 + 1]; |
| a[i][j] = x0r + x2r; |
| a[i][j + 1] = x0i + x2i; |
| x0r -= x2r; |
| x0i -= x2i; |
| a[i][j2] = wk2r * x0r - wk2i * x0i; |
| a[i][j2 + 1] = wk2r * x0i + wk2i * x0r; |
| x0r = x1r - x3i; |
| x0i = x1i + x3r; |
| a[i][j1] = wk1r * x0r - wk1i * x0i; |
| a[i][j1 + 1] = wk1r * x0i + wk1i * x0r; |
| x0r = x1r + x3i; |
| x0i = x1i - x3r; |
| a[i][j3] = wk3r * x0r - wk3i * x0i; |
| a[i][j3 + 1] = wk3r * x0i + wk3i * x0r; |
| } |
| } |
| } |
| l = m; |
| } |
| if (l < n) { |
| for (j = 0; j <= l - 2; j += 2) { |
| j1 = j + l; |
| x0r = a[i][j] - a[i][j1]; |
| x0i = a[i][j + 1] - a[i][j1 + 1]; |
| a[i][j] += a[i][j1]; |
| a[i][j + 1] += a[i][j1 + 1]; |
| a[i][j1] = x0r; |
| a[i][j1 + 1] = x0i; |
| } |
| } |
| } |
| } |
| |
| |
| void cftbrow(int n, int n2, double **a, double *w) |
| { |
| int i, j, j1, j2, j3, k, k1, ks, l, m; |
| double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; |
| double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; |
| |
| l = 1; |
| while ((l << 1) < n) { |
| m = l << 2; |
| for (j = 0; j <= l - 1; j++) { |
| j1 = j + l; |
| j2 = j1 + l; |
| j3 = j2 + l; |
| for (i = 0; i <= n2 - 2; i += 2) { |
| x0r = a[j][i] + a[j1][i]; |
| x0i = a[j][i + 1] + a[j1][i + 1]; |
| x1r = a[j][i] - a[j1][i]; |
| x1i = a[j][i + 1] - a[j1][i + 1]; |
| x2r = a[j2][i] + a[j3][i]; |
| x2i = a[j2][i + 1] + a[j3][i + 1]; |
| x3r = a[j2][i] - a[j3][i]; |
| x3i = a[j2][i + 1] - a[j3][i + 1]; |
| a[j][i] = x0r + x2r; |
| a[j][i + 1] = x0i + x2i; |
| a[j2][i] = x0r - x2r; |
| a[j2][i + 1] = x0i - x2i; |
| a[j1][i] = x1r - x3i; |
| a[j1][i + 1] = x1i + x3r; |
| a[j3][i] = x1r + x3i; |
| a[j3][i + 1] = x1i - x3r; |
| } |
| } |
| if (m < n) { |
| wk1r = w[2]; |
| for (j = m; j <= l + m - 1; j++) { |
| j1 = j + l; |
| j2 = j1 + l; |
| j3 = j2 + l; |
| for (i = 0; i <= n2 - 2; i += 2) { |
| x0r = a[j][i] + a[j1][i]; |
| x0i = a[j][i + 1] + a[j1][i + 1]; |
| x1r = a[j][i] - a[j1][i]; |
| x1i = a[j][i + 1] - a[j1][i + 1]; |
| x2r = a[j2][i] + a[j3][i]; |
| x2i = a[j2][i + 1] + a[j3][i + 1]; |
| x3r = a[j2][i] - a[j3][i]; |
| x3i = a[j2][i + 1] - a[j3][i + 1]; |
| a[j][i] = x0r + x2r; |
| a[j][i + 1] = x0i + x2i; |
| a[j2][i] = x2i - x0i; |
| a[j2][i + 1] = x0r - x2r; |
| x0r = x1r - x3i; |
| x0i = x1i + x3r; |
| a[j1][i] = wk1r * (x0r - x0i); |
| a[j1][i + 1] = wk1r * (x0r + x0i); |
| x0r = x3i + x1r; |
| x0i = x3r - x1i; |
| a[j3][i] = wk1r * (x0i - x0r); |
| a[j3][i + 1] = wk1r * (x0i + x0r); |
| } |
| } |
| k1 = 1; |
| ks = -1; |
| for (k = (m << 1); k <= n - m; k += m) { |
| k1++; |
| ks = -ks; |
| wk1r = w[k1 << 1]; |
| wk1i = w[(k1 << 1) + 1]; |
| wk2r = ks * w[k1]; |
| wk2i = w[k1 + ks]; |
| wk3r = wk1r - 2 * wk2i * wk1i; |
| wk3i = 2 * wk2i * wk1r - wk1i; |
| for (j = k; j <= l + k - 1; j++) { |
| j1 = j + l; |
| j2 = j1 + l; |
| j3 = j2 + l; |
| for (i = 0; i <= n2 - 2; i += 2) { |
| x0r = a[j][i] + a[j1][i]; |
| x0i = a[j][i + 1] + a[j1][i + 1]; |
| x1r = a[j][i] - a[j1][i]; |
| x1i = a[j][i + 1] - a[j1][i + 1]; |
| x2r = a[j2][i] + a[j3][i]; |
| x2i = a[j2][i + 1] + a[j3][i + 1]; |
| x3r = a[j2][i] - a[j3][i]; |
| x3i = a[j2][i + 1] - a[j3][i + 1]; |
| a[j][i] = x0r + x2r; |
| a[j][i + 1] = x0i + x2i; |
| x0r -= x2r; |
| x0i -= x2i; |
| a[j2][i] = wk2r * x0r - wk2i * x0i; |
| a[j2][i + 1] = wk2r * x0i + wk2i * x0r; |
| x0r = x1r - x3i; |
| x0i = x1i + x3r; |
| a[j1][i] = wk1r * x0r - wk1i * x0i; |
| a[j1][i + 1] = wk1r * x0i + wk1i * x0r; |
| x0r = x1r + x3i; |
| x0i = x1i - x3r; |
| a[j3][i] = wk3r * x0r - wk3i * x0i; |
| a[j3][i + 1] = wk3r * x0i + wk3i * x0r; |
| } |
| } |
| } |
| } |
| l = m; |
| } |
| if (l < n) { |
| for (j = 0; j <= l - 1; j++) { |
| j1 = j + l; |
| for (i = 0; i <= n2 - 2; i += 2) { |
| x0r = a[j][i] - a[j1][i]; |
| x0i = a[j][i + 1] - a[j1][i + 1]; |
| a[j][i] += a[j1][i]; |
| a[j][i + 1] += a[j1][i + 1]; |
| a[j1][i] = x0r; |
| a[j1][i + 1] = x0i; |
| } |
| } |
| } |
| } |
| |
| |
| void cftfcol(int n1, int n, double **a, double *w) |
| { |
| int i, j, j1, j2, j3, k, k1, ks, l, m; |
| double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; |
| double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; |
| |
| for (i = 0; i <= n1 - 1; i++) { |
| l = 2; |
| while ((l << 1) < n) { |
| m = l << 2; |
| for (j = 0; j <= l - 2; j += 2) { |
| j1 = j + l; |
| j2 = j1 + l; |
| j3 = j2 + l; |
| x0r = a[i][j] + a[i][j1]; |
| x0i = a[i][j + 1] + a[i][j1 + 1]; |
| x1r = a[i][j] - a[i][j1]; |
| x1i = a[i][j + 1] - a[i][j1 + 1]; |
| x2r = a[i][j2] + a[i][j3]; |
| x2i = a[i][j2 + 1] + a[i][j3 + 1]; |
| x3r = a[i][j2] - a[i][j3]; |
| x3i = a[i][j2 + 1] - a[i][j3 + 1]; |
| a[i][j] = x0r + x2r; |
| a[i][j + 1] = x0i + x2i; |
| a[i][j2] = x0r - x2r; |
| a[i][j2 + 1] = x0i - x2i; |
| a[i][j1] = x1r + x3i; |
| a[i][j1 + 1] = x1i - x3r; |
| a[i][j3] = x1r - x3i; |
| a[i][j3 + 1] = x1i + x3r; |
| } |
| if (m < n) { |
| wk1r = w[2]; |
| for (j = m; j <= l + m - 2; j += 2) { |
| j1 = j + l; |
| j2 = j1 + l; |
| j3 = j2 + l; |
| x0r = a[i][j] + a[i][j1]; |
| x0i = a[i][j + 1] + a[i][j1 + 1]; |
| x1r = a[i][j] - a[i][j1]; |
| x1i = a[i][j + 1] - a[i][j1 + 1]; |
| x2r = a[i][j2] + a[i][j3]; |
| x2i = a[i][j2 + 1] + a[i][j3 + 1]; |
| x3r = a[i][j2] - a[i][j3]; |
| x3i = a[i][j2 + 1] - a[i][j3 + 1]; |
| a[i][j] = x0r + x2r; |
| a[i][j + 1] = x0i + x2i; |
| a[i][j2] = x0i - x2i; |
| a[i][j2 + 1] = x2r - x0r; |
| x0r = x1r + x3i; |
| x0i = x1i - x3r; |
| a[i][j1] = wk1r * (x0i + x0r); |
| a[i][j1 + 1] = wk1r * (x0i - x0r); |
| x0r = x3i - x1r; |
| x0i = x3r + x1i; |
| a[i][j3] = wk1r * (x0r + x0i); |
| a[i][j3 + 1] = wk1r * (x0r - x0i); |
| } |
| k1 = 1; |
| ks = -1; |
| for (k = (m << 1); k <= n - m; k += m) { |
| k1++; |
| ks = -ks; |
| wk1r = w[k1 << 1]; |
| wk1i = w[(k1 << 1) + 1]; |
| wk2r = ks * w[k1]; |
| wk2i = w[k1 + ks]; |
| wk3r = wk1r - 2 * wk2i * wk1i; |
| wk3i = 2 * wk2i * wk1r - wk1i; |
| for (j = k; j <= l + k - 2; j += 2) { |
| j1 = j + l; |
| j2 = j1 + l; |
| j3 = j2 + l; |
| x0r = a[i][j] + a[i][j1]; |
| x0i = a[i][j + 1] + a[i][j1 + 1]; |
| x1r = a[i][j] - a[i][j1]; |
| x1i = a[i][j + 1] - a[i][j1 + 1]; |
| x2r = a[i][j2] + a[i][j3]; |
| x2i = a[i][j2 + 1] + a[i][j3 + 1]; |
| x3r = a[i][j2] - a[i][j3]; |
| x3i = a[i][j2 + 1] - a[i][j3 + 1]; |
| a[i][j] = x0r + x2r; |
| a[i][j + 1] = x0i + x2i; |
| x0r -= x2r; |
| x0i -= x2i; |
| a[i][j2] = wk2r * x0r + wk2i * x0i; |
| a[i][j2 + 1] = wk2r * x0i - wk2i * x0r; |
| x0r = x1r + x3i; |
| x0i = x1i - x3r; |
| a[i][j1] = wk1r * x0r + wk1i * x0i; |
| a[i][j1 + 1] = wk1r * x0i - wk1i * x0r; |
| x0r = x1r - x3i; |
| x0i = x1i + x3r; |
| a[i][j3] = wk3r * x0r + wk3i * x0i; |
| a[i][j3 + 1] = wk3r * x0i - wk3i * x0r; |
| } |
| } |
| } |
| l = m; |
| } |
| if (l < n) { |
| for (j = 0; j <= l - 2; j += 2) { |
| j1 = j + l; |
| x0r = a[i][j] - a[i][j1]; |
| x0i = a[i][j + 1] - a[i][j1 + 1]; |
| a[i][j] += a[i][j1]; |
| a[i][j + 1] += a[i][j1 + 1]; |
| a[i][j1] = x0r; |
| a[i][j1 + 1] = x0i; |
| } |
| } |
| } |
| } |
| |
| |
| void cftfrow(int n, int n2, double **a, double *w) |
| { |
| int i, j, j1, j2, j3, k, k1, ks, l, m; |
| double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; |
| double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; |
| |
| l = 1; |
| while ((l << 1) < n) { |
| m = l << 2; |
| for (j = 0; j <= l - 1; j++) { |
| j1 = j + l; |
| j2 = j1 + l; |
| j3 = j2 + l; |
| for (i = 0; i <= n2 - 2; i += 2) { |
| x0r = a[j][i] + a[j1][i]; |
| x0i = a[j][i + 1] + a[j1][i + 1]; |
| x1r = a[j][i] - a[j1][i]; |
| x1i = a[j][i + 1] - a[j1][i + 1]; |
| x2r = a[j2][i] + a[j3][i]; |
| x2i = a[j2][i + 1] + a[j3][i + 1]; |
| x3r = a[j2][i] - a[j3][i]; |
| x3i = a[j2][i + 1] - a[j3][i + 1]; |
| a[j][i] = x0r + x2r; |
| a[j][i + 1] = x0i + x2i; |
| a[j2][i] = x0r - x2r; |
| a[j2][i + 1] = x0i - x2i; |
| a[j1][i] = x1r + x3i; |
| a[j1][i + 1] = x1i - x3r; |
| a[j3][i] = x1r - x3i; |
| a[j3][i + 1] = x1i + x3r; |
| } |
| } |
| if (m < n) { |
| wk1r = w[2]; |
| for (j = m; j <= l + m - 1; j++) { |
| j1 = j + l; |
| j2 = j1 + l; |
| j3 = j2 + l; |
| for (i = 0; i <= n2 - 2; i += 2) { |
| x0r = a[j][i] + a[j1][i]; |
| x0i = a[j][i + 1] + a[j1][i + 1]; |
| x1r = a[j][i] - a[j1][i]; |
| x1i = a[j][i + 1] - a[j1][i + 1]; |
| x2r = a[j2][i] + a[j3][i]; |
| x2i = a[j2][i + 1] + a[j3][i + 1]; |
| x3r = a[j2][i] - a[j3][i]; |
| x3i = a[j2][i + 1] - a[j3][i + 1]; |
| a[j][i] = x0r + x2r; |
| a[j][i + 1] = x0i + x2i; |
| a[j2][i] = x0i - x2i; |
| a[j2][i + 1] = x2r - x0r; |
| x0r = x1r + x3i; |
| x0i = x1i - x3r; |
| a[j1][i] = wk1r * (x0i + x0r); |
| a[j1][i + 1] = wk1r * (x0i - x0r); |
| x0r = x3i - x1r; |
| x0i = x3r + x1i; |
| a[j3][i] = wk1r * (x0r + x0i); |
| a[j3][i + 1] = wk1r * (x0r - x0i); |
| } |
| } |
| k1 = 1; |
| ks = -1; |
| for (k = (m << 1); k <= n - m; k += m) { |
| k1++; |
| ks = -ks; |
| wk1r = w[k1 << 1]; |
| wk1i = w[(k1 << 1) + 1]; |
| wk2r = ks * w[k1]; |
| wk2i = w[k1 + ks]; |
| wk3r = wk1r - 2 * wk2i * wk1i; |
| wk3i = 2 * wk2i * wk1r - wk1i; |
| for (j = k; j <= l + k - 1; j++) { |
| j1 = j + l; |
| j2 = j1 + l; |
| j3 = j2 + l; |
| for (i = 0; i <= n2 - 2; i += 2) { |
| x0r = a[j][i] + a[j1][i]; |
| x0i = a[j][i + 1] + a[j1][i + 1]; |
| x1r = a[j][i] - a[j1][i]; |
| x1i = a[j][i + 1] - a[j1][i + 1]; |
| x2r = a[j2][i] + a[j3][i]; |
| x2i = a[j2][i + 1] + a[j3][i + 1]; |
| x3r = a[j2][i] - a[j3][i]; |
| x3i = a[j2][i + 1] - a[j3][i + 1]; |
| a[j][i] = x0r + x2r; |
| a[j][i + 1] = x0i + x2i; |
| x0r -= x2r; |
| x0i -= x2i; |
| a[j2][i] = wk2r * x0r + wk2i * x0i; |
| a[j2][i + 1] = wk2r * x0i - wk2i * x0r; |
| x0r = x1r + x3i; |
| x0i = x1i - x3r; |
| a[j1][i] = wk1r * x0r + wk1i * x0i; |
| a[j1][i + 1] = wk1r * x0i - wk1i * x0r; |
| x0r = x1r - x3i; |
| x0i = x1i + x3r; |
| a[j3][i] = wk3r * x0r + wk3i * x0i; |
| a[j3][i + 1] = wk3r * x0i - wk3i * x0r; |
| } |
| } |
| } |
| } |
| l = m; |
| } |
| if (l < n) { |
| for (j = 0; j <= l - 1; j++) { |
| j1 = j + l; |
| for (i = 0; i <= n2 - 2; i += 2) { |
| x0r = a[j][i] - a[j1][i]; |
| x0i = a[j][i + 1] - a[j1][i + 1]; |
| a[j][i] += a[j1][i]; |
| a[j][i + 1] += a[j1][i + 1]; |
| a[j1][i] = x0r; |
| a[j1][i + 1] = x0i; |
| } |
| } |
| } |
| } |
| |
| |
| void rftbcol(int n1, int n, double **a, int nc, double *c) |
| { |
| int i, j, k, kk, ks; |
| double wkr, wki, xr, xi, yr, yi; |
| |
| ks = (nc << 2) / n; |
| for (i = 0; i <= n1 - 1; i++) { |
| kk = 0; |
| for (k = (n >> 1) - 2; k >= 2; k -= 2) { |
| j = n - k; |
| kk += ks; |
| wkr = 0.5 - c[kk]; |
| wki = c[nc - kk]; |
| xr = a[i][k] - a[i][j]; |
| xi = a[i][k + 1] + a[i][j + 1]; |
| yr = wkr * xr - wki * xi; |
| yi = wkr * xi + wki * xr; |
| a[i][k] -= yr; |
| a[i][k + 1] -= yi; |
| a[i][j] += yr; |
| a[i][j + 1] -= yi; |
| } |
| } |
| } |
| |
| |
| void rftfcol(int n1, int n, double **a, int nc, double *c) |
| { |
| int i, j, k, kk, ks; |
| double wkr, wki, xr, xi, yr, yi; |
| |
| ks = (nc << 2) / n; |
| for (i = 0; i <= n1 - 1; i++) { |
| kk = 0; |
| for (k = (n >> 1) - 2; k >= 2; k -= 2) { |
| j = n - k; |
| kk += ks; |
| wkr = 0.5 - c[kk]; |
| wki = c[nc - kk]; |
| xr = a[i][k] - a[i][j]; |
| xi = a[i][k + 1] + a[i][j + 1]; |
| yr = wkr * xr + wki * xi; |
| yi = wkr * xi - wki * xr; |
| a[i][k] -= yr; |
| a[i][k + 1] -= yi; |
| a[i][j] += yr; |
| a[i][j + 1] -= yi; |
| } |
| } |
| } |
| |
| |
| void dctbsub(int n1, int n2, double **a, int nc, double *c) |
| { |
| int kk1, kk2, ks1, ks2, n1h, j1, k1, k2; |
| double w1r, w1i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i; |
| |
| ks1 = nc / n1; |
| ks2 = nc / n2; |
| n1h = n1 >> 1; |
| kk1 = ks1; |
| for (k1 = 1; k1 <= n1h - 1; k1++) { |
| j1 = n1 - k1; |
| w1r = 2 * c[kk1]; |
| w1i = 2 * c[nc - kk1]; |
| kk1 += ks1; |
| kk2 = ks2; |
| for (k2 = 2; k2 <= n2 - 2; k2 += 2) { |
| x0r = w1r * c[kk2]; |
| x0i = w1i * c[kk2]; |
| x1r = w1r * c[nc - kk2]; |
| x1i = w1i * c[nc - kk2]; |
| wkr = x0r - x1i; |
| wki = x0i + x1r; |
| wji = x0r + x1i; |
| wjr = x0i - x1r; |
| kk2 += ks2; |
| x0r = wkr * a[k1][k2] - wki * a[k1][k2 + 1]; |
| x0i = wkr * a[k1][k2 + 1] + wki * a[k1][k2]; |
| x1r = wjr * a[j1][k2] - wji * a[j1][k2 + 1]; |
| x1i = wjr * a[j1][k2 + 1] + wji * a[j1][k2]; |
| a[k1][k2] = x0r + x1i; |
| a[k1][k2 + 1] = x0i - x1r; |
| a[j1][k2] = x1r + x0i; |
| a[j1][k2 + 1] = x1i - x0r; |
| } |
| wkr = w1r * 0.5; |
| wki = w1i * 0.5; |
| wjr = w1r * c[kk2]; |
| wji = w1i * c[kk2]; |
| x0r = a[k1][0] + a[j1][0]; |
| x0i = a[k1][1] - a[j1][1]; |
| x1r = a[k1][0] - a[j1][0]; |
| x1i = a[k1][1] + a[j1][1]; |
| a[k1][0] = wkr * x0r - wki * x0i; |
| a[k1][1] = wkr * x0i + wki * x0r; |
| a[j1][0] = -wjr * x1r + wji * x1i; |
| a[j1][1] = wjr * x1i + wji * x1r; |
| } |
| w1r = 2 * c[kk1]; |
| kk2 = ks2; |
| for (k2 = 2; k2 <= n2 - 2; k2 += 2) { |
| wkr = 2 * c[kk2]; |
| wki = 2 * c[nc - kk2]; |
| wjr = w1r * wkr; |
| wji = w1r * wki; |
| kk2 += ks2; |
| x0i = wkr * a[0][k2 + 1] + wki * a[0][k2]; |
| a[0][k2] = wkr * a[0][k2] - wki * a[0][k2 + 1]; |
| a[0][k2 + 1] = x0i; |
| x0i = wjr * a[n1h][k2 + 1] + wji * a[n1h][k2]; |
| a[n1h][k2] = wjr * a[n1h][k2] - wji * a[n1h][k2 + 1]; |
| a[n1h][k2 + 1] = x0i; |
| } |
| a[0][1] *= w1r; |
| a[n1h][0] *= w1r; |
| a[n1h][1] *= 0.5; |
| } |
| |
| |
| void dctfsub(int n1, int n2, double **a, int nc, double *c) |
| { |
| int kk1, kk2, ks1, ks2, n1h, j1, k1, k2; |
| double w1r, w1i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i; |
| |
| ks1 = nc / n1; |
| ks2 = nc / n2; |
| n1h = n1 >> 1; |
| kk1 = ks1; |
| for (k1 = 1; k1 <= n1h - 1; k1++) { |
| j1 = n1 - k1; |
| w1r = 2 * c[kk1]; |
| w1i = 2 * c[nc - kk1]; |
| kk1 += ks1; |
| kk2 = ks2; |
| for (k2 = 2; k2 <= n2 - 2; k2 += 2) { |
| x0r = w1r * c[kk2]; |
| x0i = w1i * c[kk2]; |
| x1r = w1r * c[nc - kk2]; |
| x1i = w1i * c[nc - kk2]; |
| wkr = x0r - x1i; |
| wki = x0i + x1r; |
| wji = x0r + x1i; |
| wjr = x0i - x1r; |
| kk2 += ks2; |
| x0r = a[k1][k2] - a[j1][k2 + 1]; |
| x0i = a[j1][k2] + a[k1][k2 + 1]; |
| x1r = a[j1][k2] - a[k1][k2 + 1]; |
| x1i = a[k1][k2] + a[j1][k2 + 1]; |
| a[k1][k2] = wkr * x0r + wki * x0i; |
| a[k1][k2 + 1] = wkr * x0i - wki * x0r; |
| a[j1][k2] = wjr * x1r + wji * x1i; |
| a[j1][k2 + 1] = wjr * x1i - wji * x1r; |
| } |
| x0r = 2 * c[kk2]; |
| wjr = x0r * w1r; |
| wji = x0r * w1i; |
| x0r = w1r * a[k1][0] + w1i * a[k1][1]; |
| x0i = w1r * a[k1][1] - w1i * a[k1][0]; |
| x1r = -wjr * a[j1][0] + wji * a[j1][1]; |
| x1i = wjr * a[j1][1] + wji * a[j1][0]; |
| a[k1][0] = x0r + x1r; |
| a[k1][1] = x1i + x0i; |
| a[j1][0] = x0r - x1r; |
| a[j1][1] = x1i - x0i; |
| } |
| w1r = 2 * c[kk1]; |
| kk2 = ks2; |
| for (k2 = 2; k2 <= n2 - 2; k2 += 2) { |
| wkr = 2 * c[kk2]; |
| wki = 2 * c[nc - kk2]; |
| wjr = w1r * wkr; |
| wji = w1r * wki; |
| kk2 += ks2; |
| x0i = wkr * a[0][k2 + 1] - wki * a[0][k2]; |
| a[0][k2] = wkr * a[0][k2] + wki * a[0][k2 + 1]; |
| a[0][k2 + 1] = x0i; |
| x0i = wjr * a[n1h][k2 + 1] - wji * a[n1h][k2]; |
| a[n1h][k2] = wjr * a[n1h][k2] + wji * a[n1h][k2 + 1]; |
| a[n1h][k2 + 1] = x0i; |
| } |
| w1r *= 2; |
| a[0][0] *= 2; |
| a[0][1] *= w1r; |
| a[n1h][0] *= w1r; |
| } |
| |
| |
| void dstbsub(int n1, int n2, double **a, int nc, double *c) |
| { |
| int kk1, kk2, ks1, ks2, n1h, j1, k1, k2; |
| double w1r, w1i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i; |
| |
| ks1 = nc / n1; |
| ks2 = nc / n2; |
| n1h = n1 >> 1; |
| kk1 = ks1; |
| for (k1 = 1; k1 <= n1h - 1; k1++) { |
| j1 = n1 - k1; |
| w1r = 2 * c[kk1]; |
| w1i = 2 * c[nc - kk1]; |
| kk1 += ks1; |
| kk2 = ks2; |
| for (k2 = 2; k2 <= n2 - 2; k2 += 2) { |
| x0r = w1r * c[kk2]; |
| x0i = w1i * c[kk2]; |
| x1r = w1r * c[nc - kk2]; |
| x1i = w1i * c[nc - kk2]; |
| wkr = x0r - x1i; |
| wki = x0i + x1r; |
| wji = x0r + x1i; |
| wjr = x0i - x1r; |
| kk2 += ks2; |
| x0r = wkr * a[k1][k2] - wki * a[k1][k2 + 1]; |
| x0i = wkr * a[k1][k2 + 1] + wki * a[k1][k2]; |
| x1r = wjr * a[j1][k2] - wji * a[j1][k2 + 1]; |
| x1i = wjr * a[j1][k2 + 1] + wji * a[j1][k2]; |
| a[k1][k2] = x1i - x0r; |
| a[k1][k2 + 1] = x1r + x0i; |
| a[j1][k2] = x0i - x1r; |
| a[j1][k2 + 1] = x0r + x1i; |
| } |
| wkr = w1r * 0.5; |
| wki = w1i * 0.5; |
| wjr = w1r * c[kk2]; |
| wji = w1i * c[kk2]; |
| x0r = a[k1][0] + a[j1][0]; |
| x0i = a[k1][1] - a[j1][1]; |
| x1r = a[k1][0] - a[j1][0]; |
| x1i = a[k1][1] + a[j1][1]; |
| a[k1][1] = wkr * x0r - wki * x0i; |
| a[k1][0] = wkr * x0i + wki * x0r; |
| a[j1][1] = -wjr * x1r + wji * x1i; |
| a[j1][0] = wjr * x1i + wji * x1r; |
| } |
| w1r = 2 * c[kk1]; |
| kk2 = ks2; |
| for (k2 = 2; k2 <= n2 - 2; k2 += 2) { |
| wkr = 2 * c[kk2]; |
| wki = 2 * c[nc - kk2]; |
| wjr = w1r * wkr; |
| wji = w1r * wki; |
| kk2 += ks2; |
| x0i = wkr * a[0][k2 + 1] + wki * a[0][k2]; |
| a[0][k2 + 1] = wkr * a[0][k2] - wki * a[0][k2 + 1]; |
| a[0][k2] = x0i; |
| x0i = wjr * a[n1h][k2 + 1] + wji * a[n1h][k2]; |
| a[n1h][k2 + 1] = wjr * a[n1h][k2] - wji * a[n1h][k2 + 1]; |
| a[n1h][k2] = x0i; |
| } |
| a[0][1] *= w1r; |
| a[n1h][0] *= w1r; |
| a[n1h][1] *= 0.5; |
| } |
| |
| |
| void dstfsub(int n1, int n2, double **a, int nc, double *c) |
| { |
| int kk1, kk2, ks1, ks2, n1h, j1, k1, k2; |
| double w1r, w1i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i; |
| |
| ks1 = nc / n1; |
| ks2 = nc / n2; |
| n1h = n1 >> 1; |
| kk1 = ks1; |
| for (k1 = 1; k1 <= n1h - 1; k1++) { |
| j1 = n1 - k1; |
| w1r = 2 * c[kk1]; |
| w1i = 2 * c[nc - kk1]; |
| kk1 += ks1; |
| kk2 = ks2; |
| for (k2 = 2; k2 <= n2 - 2; k2 += 2) { |
| x0r = w1r * c[kk2]; |
| x0i = w1i * c[kk2]; |
| x1r = w1r * c[nc - kk2]; |
| x1i = w1i * c[nc - kk2]; |
| wkr = x0r - x1i; |
| wki = x0i + x1r; |
| wji = x0r + x1i; |
| wjr = x0i - x1r; |
| kk2 += ks2; |
| x0r = a[j1][k2 + 1] - a[k1][k2]; |
| x0i = a[k1][k2 + 1] + a[j1][k2]; |
| x1r = a[k1][k2 + 1] - a[j1][k2]; |
| x1i = a[j1][k2 + 1] + a[k1][k2]; |
| a[k1][k2] = wkr * x0r + wki * x0i; |
| a[k1][k2 + 1] = wkr * x0i - wki * x0r; |
| a[j1][k2] = wjr * x1r + wji * x1i; |
| a[j1][k2 + 1] = wjr * x1i - wji * x1r; |
| } |
| x0r = 2 * c[kk2]; |
| wjr = x0r * w1r; |
| wji = x0r * w1i; |
| x0r = w1r * a[k1][1] + w1i * a[k1][0]; |
| x0i = w1r * a[k1][0] - w1i * a[k1][1]; |
| x1r = -wjr * a[j1][1] + wji * a[j1][0]; |
| x1i = wjr * a[j1][0] + wji * a[j1][1]; |
| a[k1][0] = x0r + x1r; |
| a[k1][1] = x1i + x0i; |
| a[j1][0] = x0r - x1r; |
| a[j1][1] = x1i - x0i; |
| } |
| w1r = 2 * c[kk1]; |
| kk2 = ks2; |
| for (k2 = 2; k2 <= n2 - 2; k2 += 2) { |
| wkr = 2 * c[kk2]; |
| wki = 2 * c[nc - kk2]; |
| wjr = w1r * wkr; |
| wji = w1r * wki; |
| kk2 += ks2; |
| x0i = wkr * a[0][k2] - wki * a[0][k2 + 1]; |
| a[0][k2] = wkr * a[0][k2 + 1] + wki * a[0][k2]; |
| a[0][k2 + 1] = x0i; |
| x0i = wjr * a[n1h][k2] - wji * a[n1h][k2 + 1]; |
| a[n1h][k2] = wjr * a[n1h][k2 + 1] + wji * a[n1h][k2]; |
| a[n1h][k2 + 1] = x0i; |
| } |
| w1r *= 2; |
| a[0][0] *= 2; |
| a[0][1] *= w1r; |
| a[n1h][0] *= w1r; |
| } |
| |
