?? fftsg_h.c
字號:
/*Fast Fourier/Cosine/Sine Transform dimension :one data length :power of 2 decimation :frequency radix :split-radix 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 *);macro definitions USE_CDFT_PTHREADS : default=not defined CDFT_THREADS_BEGIN_N : must be >= 512, default=8192 CDFT_4THREADS_BEGIN_N : must be >= 512, default=65536 USE_CDFT_WINTHREADS : default=not defined CDFT_THREADS_BEGIN_N : must be >= 512, default=32768 CDFT_4THREADS_BEGIN_N : must be >= 512, default=524288-------- 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 cftfsub(int n, double *a); void cftbsub(int n, double *a); if (isgn >= 0) { cftfsub(n, a); } else { cftbsub(n, a); }}void rdft(int n, int isgn, 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)
{ 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); cftbsub(n, a); }
else if (n == 4)
{ cftbsub(n, a); } }}void ddct(int n, int isgn, 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); cftbsub(n, a); } else if (n == 4) { cftbsub(n, a); } } if (n > 4) { dctsub(n, a); } else { dctsub4(n, a); } if (isgn >= 0) { if (n > 4) { 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 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); cftbsub(n, a); } else if (n == 4) { cftbsub(n, a); } } if (n > 4) { dstsub(n, a); } else { dstsub4(n, a); } if (isgn >= 0) { if (n > 4) { 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); if (m > 2) { 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; if (n > 2) { 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); if (m > 2) { 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; if (n > 2) { 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 WR1250 /* cos(M_PI_2*0.1250) */#define WR1250 0.980785280403230449126182236134239036973933730893#endif#ifndef WI1250 /* sin(M_PI_2*0.1250) */#define WI1250 0.195090322016128267848284868477022240927691617751#endif#ifndef WR3750 /* cos(M_PI_2*0.3750) */#define WR3750 0.831469612302545237078788377617905756738560811987#endif#ifndef WI3750 /* sin(M_PI_2*0.3750) */#define WI3750 0.555570233019602224742830813948532874374937190754#endif?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -