package ooura:bindings

Source

-------- Complex DFT (Discrete Fourier Transform) -------- [definition] <case1> X[k] = sum_j=0^n-1 x[j]exp(2piij*k/n), 0<=k<n <case2> X[k] = sum_j=0^n-1 x[j]exp(-2piij*k/n), 0<=k<n (notes: sum_j=0^n-1 is a summation from j=0 to n-1) [usage] <case1> ip[0] = 0; // first time only cdft(2*n, 1, a, ip, w); <case2> ip[0] = 0; // first time only cdft(2*n, -1, a, ip, w); [parameters] 2*n :data length (int) n >= 1, n = power of 2 a[0...2n-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 ip[0...] :work area for bit reversal (int ) length of ip >= 2+sqrt(n) strictly, length of ip >= 2+(1<<(int)(log(n+0.5)/log(2))/2). ip[0],ip[1] are pointers of the cos/sin table. w[0...n/2-1] :cos/sin table (double *) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of cdft(2*n, -1, a, ip, w); is cdft(2*n, 1, a, ip, w); for (j = 0; j <= 2 * n - 1; j++) { a[j] *= 1.0 / n; }

-------- DCT (Discrete Cosine Transform) / Inverse of DCT -------- [definition] <case1> IDCT (excluding scale) C[k] = sum_j=0^n-1 a[j]cos(pij*(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> ip[0] = 0; // first time only ddct(n, 1, a, ip, w); <case2> ip[0] = 0; // first time only ddct(n, -1, a, ip, w); [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 ip[0...] :work area for bit reversal (int ) length of ip >= 2+sqrt(n/2) strictly, length of ip >= 2+(1<<(int)(log(n/2+0.5)/log(2))/2). ip[0],ip[1] are pointers of the cos/sin table. w[0...n5/4-1] :cos/sin table (double ) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of ddct(n, -1, a, ip, w); is a[0] *= 0.5; ddct(n, 1, a, ip, w); 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(pij*(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> ip[0] = 0; // first time only ddst(n, 1, a, ip, w); <case2> ip[0] = 0; // first time only ddst(n, -1, a, ip, w); [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] ip[0...] :work area for bit reversal (int ) length of ip >= 2+sqrt(n/2) strictly, length of ip >= 2+(1<<(int)(log(n/2+0.5)/log(2))/2). ip[0],ip[1] are pointers of the cos/sin table. w[0...n5/4-1] :cos/sin table (double ) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of ddst(n, -1, a, ip, w); is a[0] *= 0.5; ddst(n, 1, a, ip, w); 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(pij*k/n), 0<=k<=n [usage] ip[0] = 0; // first time only dfct(n, a, t, ip, w); [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 t[0...n/2] :work area (double *) ip[0...] :work area for bit reversal (int ) length of ip >= 2+sqrt(n/4) strictly, length of ip >= 2+(1<<(int)(log(n/4+0.5)/log(2))/2). ip[0],ip[1] are pointers of the cos/sin table. w[0...n5/8-1] :cos/sin table (double ) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of a[0] *= 0.5; a[n] *= 0.5; dfct(n, a, t, ip, w); is a[0] *= 0.5; a[n] *= 0.5; dfct(n, a, t, ip, w); 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(pij*k/n), 0<k<n [usage] ip[0] = 0; // first time only dfst(n, a, t, ip, w); [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) t[0...n/2-1] :work area (double *) ip[0...] :work area for bit reversal (int ) length of ip >= 2+sqrt(n/4) strictly, length of ip >= 2+(1<<(int)(log(n/4+0.5)/log(2))/2). ip[0],ip[1] are pointers of the cos/sin table. w[0...n5/8-1] :cos/sin table (double ) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of dfst(n, a, t, ip, w); is dfst(n, a, t, ip, w); for (j = 1; j <= n - 1; j++) { a[j] *= 2.0 / n; }

-------- Real DFT / Inverse of Real DFT --------

[definition] <case1> RDFT R[k] = sum_j=0^n-1 a[j]cos(2pijk/n), 0<=k<=n/2 I[k] = sum_j=0^n-1 a[j]sin(2pijk/n), 0<k<n/2 <case2> IRDFT (excluding scale) a[k] = (R[0] + R[n/2]cos(pik))/2 + sum_j=1^n/2-1 R[j]cos(2pijk/n) + sum_j=1^n/2-1 I[j]sin(2pijk/n), 0<=k<n

[usage] <case1> ip[0] = 0; // first time only rdft(n, 1, a, ip, w); <case2> ip[0] = 0; // first time only rdft(n, -1, a, ip, w);

[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] ip[0...] :work area for bit reversal (int ) length of ip >= 2+sqrt(n/2) strictly, length of ip >= 2+(1<<(int)(log(n/2+0.5)/log(2))/2). ip[0],ip[1] are pointers of the cos/sin table. w[0...n/2-1] :cos/sin table (double *) w[],ip[] are initialized if ip[0] == 0.

[remark] Inverse of rdft(n, 1, a, ip, w); is rdft(n, -1, a, ip, w); for (j = 0; j <= n - 1; j++) { a[j] *= 2.0 / n; } .

-------- Complex DFT (Discrete Fourier Transform) -------- [definition] <case1> X[k] = sum_j=0^n-1 x[j]exp(2piij*k/n), 0<=k<n <case2> X[k] = sum_j=0^n-1 x[j]exp(-2piij*k/n), 0<=k<n (notes: sum_j=0^n-1 is a summation from j=0 to n-1) [usage] <case1> ip[0] = 0; // first time only cdft(2*n, 1, a, ip, w); <case2> ip[0] = 0; // first time only cdft(2*n, -1, a, ip, w); [parameters] 2*n :data length (int) n >= 1, n = power of 2 a[0...2n-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 ip[0...] :work area for bit reversal (int ) length of ip >= 2+sqrt(n) strictly, length of ip >= 2+(1<<(int)(log(n+0.5)/log(2))/2). ip[0],ip[1] are pointers of the cos/sin table. w[0...n/2-1] :cos/sin table (double *) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of cdft(2*n, -1, a, ip, w); is cdft(2*n, 1, a, ip, w); for (j = 0; j <= 2 * n - 1; j++) { a[j] *= 1.0 / n; }

-------- DCT (Discrete Cosine Transform) / Inverse of DCT -------- [definition] <case1> IDCT (excluding scale) C[k] = sum_j=0^n-1 a[j]cos(pij*(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> ip[0] = 0; // first time only ddct(n, 1, a, ip, w); <case2> ip[0] = 0; // first time only ddct(n, -1, a, ip, w); [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 ip[0...] :work area for bit reversal (int ) length of ip >= 2+sqrt(n/2) strictly, length of ip >= 2+(1<<(int)(log(n/2+0.5)/log(2))/2). ip[0],ip[1] are pointers of the cos/sin table. w[0...n5/4-1] :cos/sin table (double ) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of ddct(n, -1, a, ip, w); is a[0] *= 0.5; ddct(n, 1, a, ip, w); 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(pij*(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> ip[0] = 0; // first time only ddst(n, 1, a, ip, w); <case2> ip[0] = 0; // first time only ddst(n, -1, a, ip, w); [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] ip[0...] :work area for bit reversal (int ) length of ip >= 2+sqrt(n/2) strictly, length of ip >= 2+(1<<(int)(log(n/2+0.5)/log(2))/2). ip[0],ip[1] are pointers of the cos/sin table. w[0...n5/4-1] :cos/sin table (double ) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of ddst(n, -1, a, ip, w); is a[0] *= 0.5; ddst(n, 1, a, ip, w); 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(pij*k/n), 0<=k<=n [usage] ip[0] = 0; // first time only dfct(n, a, t, ip, w); [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 t[0...n/2] :work area (double *) ip[0...] :work area for bit reversal (int ) length of ip >= 2+sqrt(n/4) strictly, length of ip >= 2+(1<<(int)(log(n/4+0.5)/log(2))/2). ip[0],ip[1] are pointers of the cos/sin table. w[0...n5/8-1] :cos/sin table (double ) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of a[0] *= 0.5; a[n] *= 0.5; dfct(n, a, t, ip, w); is a[0] *= 0.5; a[n] *= 0.5; dfct(n, a, t, ip, w); 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(pij*k/n), 0<k<n [usage] ip[0] = 0; // first time only dfst(n, a, t, ip, w); [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) t[0...n/2-1] :work area (double *) ip[0...] :work area for bit reversal (int ) length of ip >= 2+sqrt(n/4) strictly, length of ip >= 2+(1<<(int)(log(n/4+0.5)/log(2))/2). ip[0],ip[1] are pointers of the cos/sin table. w[0...n5/8-1] :cos/sin table (double ) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of dfst(n, a, t, ip, w); is dfst(n, a, t, ip, w); for (j = 1; j <= n - 1; j++) { a[j] *= 2.0 / n; }

-------- Real DFT / Inverse of Real DFT -------- [definition] <case1> RDFT R[k] = sum_j=0^n-1 a[j]cos(2pijk/n), 0<=k<=n/2 I[k] = sum_j=0^n-1 a[j]sin(2pijk/n), 0<k<n/2 <case2> IRDFT (excluding scale) a[k] = (R[0] + R[n/2]cos(pik))/2 + sum_j=1^n/2-1 R[j]cos(2pijk/n) + sum_j=1^n/2-1 I[j]sin(2pijk/n), 0<=k<n [usage] <case1> ip[0] = 0; // first time only rdft(n, 1, a, ip, w); <case2> ip[0] = 0; // first time only rdft(n, -1, a, ip, w); [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] ip[0...] :work area for bit reversal (int ) length of ip >= 2+sqrt(n/2) strictly, length of ip >= 2+(1<<(int)(log(n/2+0.5)/log(2))/2). ip[0],ip[1] are pointers of the cos/sin table. w[0...n/2-1] :cos/sin table (double *) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of rdft(n, 1, a, ip, w); is rdft(n, -1, a, ip, w); for (j = 0; j <= n - 1; j++) { a[j] *= 2.0 / n; } .

-------- Complex DFT (Discrete Fourier Transform) -------- [definition] <case1> X[k] = sum_j=0^n-1 x[j]exp(2piij*k/n), 0<=k<n <case2> X[k] = sum_j=0^n-1 x[j]exp(-2piij*k/n), 0<=k<n (notes: sum_j=0^n-1 is a summation from j=0 to n-1) [usage] <case1> ip[0] = 0; // first time only cdft(2*n, 1, a, ip, w); <case2> ip[0] = 0; // first time only cdft(2*n, -1, a, ip, w); [parameters] 2*n :data length (int) n >= 1, n = power of 2 a[0...2n-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 ip[0...] :work area for bit reversal (int ) length of ip >= 2+sqrt(n) strictly, length of ip >= 2+(1<<(int)(log(n+0.5)/log(2))/2). ip[0],ip[1] are pointers of the cos/sin table. w[0...n/2-1] :cos/sin table (double *) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of cdft(2*n, -1, a, ip, w); is cdft(2*n, 1, a, ip, w); for (j = 0; j <= 2 * n - 1; j++) { a[j] *= 1.0 / n; }

-------- Complex DFT (Discrete Fourier Transform) -------- [definition] <case1> X[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 x[j1][j2] * exp(2piij1k1/n1) * exp(2piij2k2/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(-2piij1k1/n1) * exp(-2piij2k2/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, t, ip, w); <case2> ip[0] = 0; // first time only cdft2d(n1, 2*n2, -1, a, t, 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 t[0...*] :work area (double *) length of t >= 8*n1, if single thread, length of t >= 8n1FFT2D_MAX_THREADS, if multi threads, t is dynamically allocated, if t == NULL. 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, t, ip, w); is cdft2d(n1, 2*n2, 1, a, t, ip, w); for (j1 = 0; j1 <= n1 - 1; j1++) { for (j2 = 0; j2 <= 2 * n2 - 1; j2++) { a[j1][j2] *= 1.0 / n1 / n2; } }

-------- Complex DFT (Discrete Fourier Transform) -------- [definition] <case1> X[k1][k2][k3] = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 x[j1][j2][j3] * exp(2piij1k1/n1) * exp(2piij2k2/n2) * exp(2piij3k3/n3), 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 <case2> X[k1][k2][k3] = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 x[j1][j2][j3] * exp(-2piij1k1/n1) * exp(-2piij2k2/n2) * exp(-2piij3k3/n3), 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 (notes: sum_j=0^n-1 is a summation from j=0 to n-1) [usage] <case1> ip[0] = 0; // first time only cdft3d(n1, n2, 2*n3, 1, a, t, ip, w); <case2> ip[0] = 0; // first time only cdft3d(n1, n2, 2*n3, -1, a, t, ip, w); [parameters] n1 :data length (int) n1 >= 1, n1 = power of 2 n2 :data length (int) n2 >= 1, n2 = power of 2 2*n3 :data length (int) n3 >= 1, n3 = power of 2 a[0...n1-1][0...n2-1][0...2*n3-1] :input/output data (double ***) input data a[j1][j2][2*j3] = Re(x[j1][j2][j3]), a[j1][j2][2*j3+1] = Im(x[j1][j2][j3]), 0<=j1<n1, 0<=j2<n2, 0<=j3<n3 output data a[k1][k2][2*k3] = Re(X[k1][k2][k3]), a[k1][k2][2*k3+1] = Im(X[k1][k2][k3]), 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 t[0...*] :work area (double *) length of t >= max(8n1, 8n2), if single thread, length of t >= max(8n1, 8n2)*FFT3D_MAX_THREADS, if multi threads, t is dynamically allocated, if t == NULL. ip[0...*] :work area for bit reversal (int *) length of ip >= 2+sqrt(n) (n = max(n1, n2, n3)) 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, n3/2) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of cdft3d(n1, n2, 2*n3, -1, a, t, ip, w); is cdft3d(n1, n2, 2*n3, 1, a, t, ip, w); for (j1 = 0; j1 <= n1 - 1; j1++) { for (j2 = 0; j2 <= n2 - 1; j2++) { for (j3 = 0; j3 <= 2 * n3 - 1; j3++) { a[j1][j2][j3] *= 1.0 / n1 / n2 / n3; } } }

-------- DCT (Discrete Cosine Transform) / Inverse of DCT -------- [definition] <case1> IDCT (excluding scale) C[k] = sum_j=0^n-1 a[j]cos(pij*(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> ip[0] = 0; // first time only ddct(n, 1, a, ip, w); <case2> ip[0] = 0; // first time only ddct(n, -1, a, ip, w); [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 ip[0...] :work area for bit reversal (int ) length of ip >= 2+sqrt(n/2) strictly, length of ip >= 2+(1<<(int)(log(n/2+0.5)/log(2))/2). ip[0],ip[1] are pointers of the cos/sin table. w[0...n5/4-1] :cos/sin table (double ) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of ddct(n, -1, a, ip, w); is a[0] *= 0.5; ddct(n, 1, a, ip, w); 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[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 a[j1][j2] * cos(pij1(k1+1/2)/n1) * cos(pij2(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...*] :work area (double *) length of t >= 4*n1, if single thread, length of t >= 4n1FFT2D_MAX_THREADS, if multi threads, t is dynamically allocated, if t == NULL. ip[0...*] :work area for bit reversal (int *) length of ip >= 2+sqrt(n) (n = max(n1/2, n2/2)) ip[0],ip[1] are pointers of the cos/sin table. w[0...*] :cos/sin table (double *) length of w >= max(n13/2, n23/2) 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; } }

-------- DCT (Discrete Cosine Transform) / Inverse of DCT -------- [definition] <case1> IDCT (excluding scale) C[k1][k2][k3] = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 a[j1][j2][j3] * cos(pij1(k1+1/2)/n1) * cos(pij2(k2+1/2)/n2) * cos(pij3(k3+1/2)/n3), 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 <case2> DCT C[k1][k2][k3] = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 a[j1][j2][j3] * cos(pi(j1+1/2)k1/n1) * cos(pi(j2+1/2)k2/n2) * cos(pi(j3+1/2)k3/n3), 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 [usage] <case1> ip[0] = 0; // first time only ddct3d(n1, n2, n3, 1, a, t, ip, w); <case2> ip[0] = 0; // first time only ddct3d(n1, n2, n3, -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 n3 :data length (int) n3 >= 2, n3 = power of 2 a[0...n1-1][0...n2-1][0...n3-1] :input/output data (double ***) output data a[k1][k2][k3] = C[k1][k2][k3], 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 t[0...*] :work area (double *) length of t >= max(4n1, 4n2), if single thread, length of t >= max(4n1, 4n2)*FFT3D_MAX_THREADS, if multi threads, t is dynamically allocated, if t == NULL. ip[0...*] :work area for bit reversal (int *) length of ip >= 2+sqrt(n) (n = max(n1/2, n2/2, n3/2)) ip[0],ip[1] are pointers of the cos/sin table. w[0...*] :cos/sin table (double *) length of w >= max(n13/2, n23/2, n3*3/2) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of ddct3d(n1, n2, n3, -1, a, t, ip, w); is for (j1 = 0; j1 <= n1 - 1; j1++) { for (j2 = 0; j2 <= n2 - 1; j2++) { a[j1][j2][0] *= 0.5; } for (j3 = 0; j3 <= n3 - 1; j3++) { a[j1][0][j3] *= 0.5; } } for (j2 = 0; j2 <= n2 - 1; j2++) { for (j3 = 0; j3 <= n3 - 1; j3++) { a[0][j2][j3] *= 0.5; } } ddct3d(n1, n2, n3, 1, a, t, ip, w); for (j1 = 0; j1 <= n1 - 1; j1++) { for (j2 = 0; j2 <= n2 - 1; j2++) { for (j3 = 0; j3 <= n3 - 1; j3++) { a[j1][j2][j3] *= 8.0 / n1 / n2 / n3; } } }

-------- DST (Discrete Sine Transform) / Inverse of DST -------- [definition] <case1> IDST (excluding scale) S[k] = sum_j=1^n A[j]sin(pij*(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> ip[0] = 0; // first time only ddst(n, 1, a, ip, w); <case2> ip[0] = 0; // first time only ddst(n, -1, a, ip, w); [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] ip[0...] :work area for bit reversal (int ) length of ip >= 2+sqrt(n/2) strictly, length of ip >= 2+(1<<(int)(log(n/2+0.5)/log(2))/2). ip[0],ip[1] are pointers of the cos/sin table. w[0...n5/4-1] :cos/sin table (double ) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of ddst(n, -1, a, ip, w); is a[0] *= 0.5; ddst(n, 1, a, ip, w); 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[k1][k2] = sum_j1=1^n1 sum_j2=1^n2 A[j1][j2] * sin(pij1(k1+1/2)/n1) * sin(pij2(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...*] :work area (double *) length of t >= 4*n1, if single thread, length of t >= 4n1FFT2D_MAX_THREADS, if multi threads, t is dynamically allocated, if t == NULL. ip[0...*] :work area for bit reversal (int *) length of ip >= 2+sqrt(n) (n = max(n1/2, n2/2)) ip[0],ip[1] are pointers of the cos/sin table. w[0...*] :cos/sin table (double *) length of w >= max(n13/2, n23/2) 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; } }

-------- DST (Discrete Sine Transform) / Inverse of DST -------- [definition] <case1> IDST (excluding scale) S[k1][k2][k3] = sum_j1=1^n1 sum_j2=1^n2 sum_j3=1^n3 A[j1][j2][j3] * sin(pij1(k1+1/2)/n1) * sin(pij2(k2+1/2)/n2) * sin(pij3(k3+1/2)/n3), 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 <case2> DST S[k1][k2][k3] = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 a[j1][j2][j3] * sin(pi(j1+1/2)k1/n1) * sin(pi(j2+1/2)k2/n2) * sin(pi(j3+1/2)k3/n3), 0<k1<=n1, 0<k2<=n2, 0<k3<=n3 [usage] <case1> ip[0] = 0; // first time only ddst3d(n1, n2, n3, 1, a, t, ip, w); <case2> ip[0] = 0; // first time only ddst3d(n1, n2, n3, -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 n3 :data length (int) n3 >= 2, n3 = power of 2 a[0...n1-1][0...n2-1][0...n3-1] :input/output data (double ***) <case1> input data a[j1%n1][j2%n2][j3%n3] = A[j1][j2][j3], 0<j1<=n1, 0<j2<=n2, 0<j3<=n3, (n%m : n mod m) output data a[k1][k2][k3] = S[k1][k2][k3], 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 <case2> output data a[k1%n1][k2%n2][k3%n3] = S[k1][k2][k3], 0<k1<=n1, 0<k2<=n2, 0<k3<=n3 t[0...*] :work area (double *) length of t >= max(4n1, 4n2), if single thread, length of t >= max(4n1, 4n2)*FFT3D_MAX_THREADS, if multi threads, t is dynamically allocated, if t == NULL. ip[0...*] :work area for bit reversal (int *) length of ip >= 2+sqrt(n) (n = max(n1/2, n2/2, n3/2)) ip[0],ip[1] are pointers of the cos/sin table. w[0...*] :cos/sin table (double *) length of w >= max(n13/2, n23/2, n3*3/2) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of ddst3d(n1, n2, n3, -1, a, t, ip, w); is for (j1 = 0; j1 <= n1 - 1; j1++) { for (j2 = 0; j2 <= n2 - 1; j2++) { a[j1][j2][0] *= 0.5; } for (j3 = 0; j3 <= n3 - 1; j3++) { a[j1][0][j3] *= 0.5; } } for (j2 = 0; j2 <= n2 - 1; j2++) { for (j3 = 0; j3 <= n3 - 1; j3++) { a[0][j2][j3] *= 0.5; } } ddst3d(n1, n2, n3, 1, a, t, ip, w); for (j1 = 0; j1 <= n1 - 1; j1++) { for (j2 = 0; j2 <= n2 - 1; j2++) { for (j3 = 0; j3 <= n3 - 1; j3++) { a[j1][j2][j3] *= 8.0 / n1 / n2 / n3; } } }

-------- Cosine Transform of RDFT (Real Symmetric DFT) -------- [definition] C[k] = sum_j=0^n a[j]cos(pij*k/n), 0<=k<=n [usage] ip[0] = 0; // first time only dfct(n, a, t, ip, w); [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 t[0...n/2] :work area (double *) ip[0...] :work area for bit reversal (int ) length of ip >= 2+sqrt(n/4) strictly, length of ip >= 2+(1<<(int)(log(n/4+0.5)/log(2))/2). ip[0],ip[1] are pointers of the cos/sin table. w[0...n5/8-1] :cos/sin table (double ) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of a[0] *= 0.5; a[n] *= 0.5; dfct(n, a, t, ip, w); is a[0] *= 0.5; a[n] *= 0.5; dfct(n, a, t, ip, w); 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(pij*k/n), 0<k<n [usage] ip[0] = 0; // first time only dfst(n, a, t, ip, w); [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) t[0...n/2-1] :work area (double *) ip[0...] :work area for bit reversal (int ) length of ip >= 2+sqrt(n/4) strictly, length of ip >= 2+(1<<(int)(log(n/4+0.5)/log(2))/2). ip[0],ip[1] are pointers of the cos/sin table. w[0...n5/8-1] :cos/sin table (double ) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of dfst(n, a, t, ip, w); is dfst(n, a, t, ip, w); for (j = 1; j <= n - 1; j++) { a[j] *= 2.0 / n; }

-------- Real DFT / Inverse of Real DFT -------- [definition] <case1> RDFT R[k] = sum_j=0^n-1 a[j]cos(2pijk/n), 0<=k<=n/2 I[k] = sum_j=0^n-1 a[j]sin(2pijk/n), 0<k<n/2 <case2> IRDFT (excluding scale) a[k] = (R[0] + R[n/2]cos(pik))/2 + sum_j=1^n/2-1 R[j]cos(2pijk/n) + sum_j=1^n/2-1 I[j]sin(2pijk/n), 0<=k<n [usage] <case1> ip[0] = 0; // first time only rdft(n, 1, a, ip, w); <case2> ip[0] = 0; // first time only rdft(n, -1, a, ip, w); [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] ip[0...] :work area for bit reversal (int ) length of ip >= 2+sqrt(n/2) strictly, length of ip >= 2+(1<<(int)(log(n/2+0.5)/log(2))/2). ip[0],ip[1] are pointers of the cos/sin table. w[0...n/2-1] :cos/sin table (double *) w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of rdft(n, 1, a, ip, w); is rdft(n, -1, a, ip, w); for (j = 0; j <= n - 1; j++) { a[j] *= 2.0 / n; } .

-------- 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(2pij1k1/n1 + 2pij2k2/n2), 0<=k1<n1, 0<=k2<n2 I[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 a[j1][j2] * sin(2pij1k1/n1 + 2pij2k2/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(2pij1k1/n1 + 2pij2k2/n2) + I[j1][j2] * sin(2pij1k1/n1 + 2pij2k2/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, t, ip, w); <case2> ip[0] = 0; // first time only rdft2d(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> 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] ---- output ordering ---- rdft2d(n1, n2, 1, a, t, ip, w); rdft2dsort(n1, n2, 1, a); // stored data is a[0...n1-1][0...n2+1]: // a[k1][2*k2] = R[k1][k2], // a[k1][2*k2+1] = I[k1][k2], // 0<=k1<n1, 0<=k2<=n2/2. // the stored data is larger than the input data! ---- input ordering ---- rdft2dsort(n1, n2, -1, a); rdft2d(n1, n2, -1, a, t, ip, w); t[0...*] :work area (double *) length of t >= 8*n1, if single thread, length of t >= 8n1FFT2D_MAX_THREADS, if multi threads, t is dynamically allocated, if t == NULL. 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, t, ip, w); is rdft2d(n1, n2, -1, a, t, ip, w); for (j1 = 0; j1 <= n1 - 1; j1++) { for (j2 = 0; j2 <= n2 - 1; j2++) { a[j1][j2] *= 2.0 / n1 / n2; } }

-------- Real DFT / Inverse of Real DFT -------- [definition] <case1> RDFT R[k1][k2][k3] = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 a[j1][j2][j3] * cos(2pij1k1/n1 + 2pij2k2/n2 + 2pij3*k3/n3), 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 I[k1][k2][k3] = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 a[j1][j2][j3] * sin(2pij1k1/n1 + 2pij2k2/n2 + 2pij3*k3/n3), 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 <case2> IRDFT (excluding scale) a[k1][k2][k3] = (1/2) * sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 (R[j1][j2][j3] * cos(2pij1k1/n1 + 2pij2k2/n2 + 2pij3*k3/n3) + I[j1][j2][j3] * sin(2pij1k1/n1 + 2pij2k2/n2 + 2pij3*k3/n3)), 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 (notes: R[(n1-k1)%n1][(n2-k2)%n2][(n3-k3)%n3] = R[k1][k2][k3], I[(n1-k1)%n1][(n2-k2)%n2][(n3-k3)%n3] = -I[k1][k2][k3], 0<=k1<n1, 0<=k2<n2, 0<=k3<n3) [usage] <case1> ip[0] = 0; // first time only rdft3d(n1, n2, n3, 1, a, t, ip, w); <case2> ip[0] = 0; // first time only rdft3d(n1, n2, n3, -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 n3 :data length (int) n3 >= 2, n3 = power of 2 a[0...n1-1][0...n2-1][0...n3-1] :input/output data (double ***) <case1> output data a[k1][k2][2*k3] = R[k1][k2][k3] = R[(n1-k1)%n1][(n2-k2)%n2][n3-k3], a[k1][k2][2*k3+1] = I[k1][k2][k3] = -I[(n1-k1)%n1][(n2-k2)%n2][n3-k3], 0<=k1<n1, 0<=k2<n2, 0<k3<n3/2, (n%m : n mod m), a[k1][k2][0] = R[k1][k2][0] = R[(n1-k1)%n1][n2-k2][0], a[k1][k2][1] = I[k1][k2][0] = -I[(n1-k1)%n1][n2-k2][0], a[k1][n2-k2][1] = R[(n1-k1)%n1][k2][n3/2] = R[k1][n2-k2][n3/2], a[k1][n2-k2][0] = -I[(n1-k1)%n1][k2][n3/2] = I[k1][n2-k2][n3/2], 0<=k1<n1, 0<k2<n2/2, a[k1][0][0] = R[k1][0][0] = R[n1-k1][0][0], a[k1][0][1] = I[k1][0][0] = -I[n1-k1][0][0], a[k1][n2/2][0] = R[k1][n2/2][0] = R[n1-k1][n2/2][0], a[k1][n2/2][1] = I[k1][n2/2][0] = -I[n1-k1][n2/2][0], a[n1-k1][0][1] = R[k1][0][n3/2] = R[n1-k1][0][n3/2], a[n1-k1][0][0] = -I[k1][0][n3/2] = I[n1-k1][0][n3/2], a[n1-k1][n2/2][1] = R[k1][n2/2][n3/2] = R[n1-k1][n2/2][n3/2], a[n1-k1][n2/2][0] = -I[k1][n2/2][n3/2] = I[n1-k1][n2/2][n3/2], 0<k1<n1/2, a[0][0][0] = R[0][0][0], a[0][0][1] = R[0][0][n3/2], a[0][n2/2][0] = R[0][n2/2][0], a[0][n2/2][1] = R[0][n2/2][n3/2], a[n1/2][0][0] = R[n1/2][0][0], a[n1/2][0][1] = R[n1/2][0][n3/2], a[n1/2][n2/2][0] = R[n1/2][n2/2][0], a[n1/2][n2/2][1] = R[n1/2][n2/2][n3/2] <case2> input data a[j1][j2][2*j3] = R[j1][j2][j3] = R[(n1-j1)%n1][(n2-j2)%n2][n3-j3], a[j1][j2][2*j3+1] = I[j1][j2][j3] = -I[(n1-j1)%n1][(n2-j2)%n2][n3-j3], 0<=j1<n1, 0<=j2<n2, 0<j3<n3/2, a[j1][j2][0] = R[j1][j2][0] = R[(n1-j1)%n1][n2-j2][0], a[j1][j2][1] = I[j1][j2][0] = -I[(n1-j1)%n1][n2-j2][0], a[j1][n2-j2][1] = R[(n1-j1)%n1][j2][n3/2] = R[j1][n2-j2][n3/2], a[j1][n2-j2][0] = -I[(n1-j1)%n1][j2][n3/2] = I[j1][n2-j2][n3/2], 0<=j1<n1, 0<j2<n2/2, a[j1][0][0] = R[j1][0][0] = R[n1-j1][0][0], a[j1][0][1] = I[j1][0][0] = -I[n1-j1][0][0], a[j1][n2/2][0] = R[j1][n2/2][0] = R[n1-j1][n2/2][0], a[j1][n2/2][1] = I[j1][n2/2][0] = -I[n1-j1][n2/2][0], a[n1-j1][0][1] = R[j1][0][n3/2] = R[n1-j1][0][n3/2], a[n1-j1][0][0] = -I[j1][0][n3/2] = I[n1-j1][0][n3/2], a[n1-j1][n2/2][1] = R[j1][n2/2][n3/2] = R[n1-j1][n2/2][n3/2], a[n1-j1][n2/2][0] = -I[j1][n2/2][n3/2] = I[n1-j1][n2/2][n3/2], 0<j1<n1/2, a[0][0][0] = R[0][0][0], a[0][0][1] = R[0][0][n3/2], a[0][n2/2][0] = R[0][n2/2][0], a[0][n2/2][1] = R[0][n2/2][n3/2], a[n1/2][0][0] = R[n1/2][0][0], a[n1/2][0][1] = R[n1/2][0][n3/2], a[n1/2][n2/2][0] = R[n1/2][n2/2][0], a[n1/2][n2/2][1] = R[n1/2][n2/2][n3/2] ---- output ordering ---- rdft3d(n1, n2, n3, 1, a, t, ip, w); rdft3dsort(n1, n2, n3, 1, a); // stored data is a[0...n1-1][0...n2-1][0...n3+1]: // a[k1][k2][2*k3] = R[k1][k2][k3], // a[k1][k2][2*k3+1] = I[k1][k2][k3], // 0<=k1<n1, 0<=k2<n2, 0<=k3<=n3/2. // the stored data is larger than the input data! ---- input ordering ---- rdft3dsort(n1, n2, n3, -1, a); rdft3d(n1, n2, n3, -1, a, t, ip, w); t[0...*] :work area (double *) length of t >= max(8n1, 8n2), if single thread, length of t >= max(8n1, 8n2)*FFT3D_MAX_THREADS, if multi threads, t is dynamically allocated, if t == NULL. ip[0...*] :work area for bit reversal (int *) length of ip >= 2+sqrt(n) (n = max(n1, n2, n3/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/2, n3/4) + n3/4 w[],ip[] are initialized if ip[0] == 0. [remark] Inverse of rdft3d(n1, n2, n3, 1, a, t, ip, w); is rdft3d(n1, n2, n3, -1, a, t, ip, w); for (j1 = 0; j1 <= n1 - 1; j1++) { for (j2 = 0; j2 <= n2 - 1; j2++) { for (j3 = 0; j3 <= n3 - 1; j3++) { a[j1][j2][j3] *= 2.0 / n1 / n2 / n3; } } }