?? fftsg.c
字號:
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) { cftfsub(m, a, ip, nw, w); rftfsub(m, a, nc, w + nw); } else if (m == 4) { cftfsub(m, a, ip, nw, 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) { cftfsub(m, t, ip, nw, w); rftfsub(m, t, nc, w + nw); } else if (m == 4) { cftfsub(m, t, ip, nw, 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; a[k - l] = t[j] - t[j + 1]; a[k + l] = t[j] + t[j + 1]; } l <<= 1; mh = m >> 1; for (j = 0; j < mh; j++) { k = m - j; t[j] = t[m + k] - t[m + j]; t[k] = t[m + k] + t[m + j]; } t[mh] = t[m + mh]; m = mh; } a[l] = t[0]; a[n] = t[2] - t[1]; a[0] = t[2] + t[1]; } else { a[1] = a[0]; a[2] = t[0]; a[0] = t[1]; }}void dfst(int n, 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 cftfsub(int n, double *a, int *ip, int nw, double *w); void rftfsub(int n, double *a, int nc, double *c); void dstsub(int n, double *a, int nc, double *c); int 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); } if (n > 2) { m = n >> 1; 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[0] = a[mh] - a[n - mh]; a[mh] += a[n - mh]; a[0] = a[m]; dstsub(m, a, nc, w + nw); if (m > 4) { cftfsub(m, a, ip, nw, w); rftfsub(m, a, nc, w + nw); } else if (m == 4) { cftfsub(m, a, ip, nw, w); } a[n - 1] = a[1] - a[0]; 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) { dstsub(m, t, nc, w + nw); if (m > 4) { cftfsub(m, t, ip, nw, w); rftfsub(m, t, nc, w + nw); } else if (m == 4) { cftfsub(m, t, ip, nw, w); } a[n - l] = t[1] - t[0]; a[l] = t[0] + t[1]; k = 0; for (j = 2; j < m; j += 2) { k += l << 2; a[k - l] = -t[j] - t[j + 1]; a[k + l] = t[j] - t[j + 1]; } l <<= 1; mh = m >> 1; for (j = 1; j < mh; j++) { k = m - j; t[j] = t[m + k] + t[m + j]; t[k] = t[m + k] - t[m + j]; } t[0] = t[m + mh]; m = mh; } a[l] = t[0]; } a[0] = 0;}/* -------- initializing routines -------- */#include <math.h>void makewt(int nw, int *ip, double *w){ void makeipt(int nw, int *ip); int j, nwh, nw0, nw1; double delta, wn4r, wk1r, wk1i, wk3r, wk3i; ip[0] = nw; ip[1] = 1; if (nw > 2) { nwh = nw >> 1; delta = atan(1.0) / nwh; wn4r = cos(delta * nwh); w[0] = 1; w[1] = wn4r; if (nwh == 4) { w[2] = cos(delta * 2); w[3] = sin(delta * 2); } else if (nwh > 4) { makeipt(nw, ip); w[2] = 0.5 / cos(delta * 2); w[3] = 0.5 / cos(delta * 6); for (j = 4; j < nwh; j += 4) { w[j] = cos(delta * j); w[j + 1] = sin(delta * j); w[j + 2] = cos(3 * delta * j); w[j + 3] = -sin(3 * delta * j); } } nw0 = 0; while (nwh > 2) { nw1 = nw0 + nwh; nwh >>= 1; w[nw1] = 1; w[nw1 + 1] = wn4r; if (nwh == 4) { wk1r = w[nw0 + 4]; wk1i = w[nw0 + 5]; w[nw1 + 2] = wk1r; w[nw1 + 3] = wk1i; } else if (nwh > 4) { wk1r = w[nw0 + 4]; wk3r = w[nw0 + 6]; w[nw1 + 2] = 0.5 / wk1r; w[nw1 + 3] = 0.5 / wk3r; for (j = 4; j < nwh; j += 4) { wk1r = w[nw0 + 2 * j]; wk1i = w[nw0 + 2 * j + 1]; wk3r = w[nw0 + 2 * j + 2]; wk3i = w[nw0 + 2 * j + 3]; w[nw1 + j] = wk1r; w[nw1 + j + 1] = wk1i; w[nw1 + j + 2] = wk3r; w[nw1 + j + 3] = wk3i; } } nw0 = nw1; } }}void makeipt(int nw, int *ip){ int j, l, m, m2, p, q; ip[2] = 0; ip[3] = 16; m = 2; for (l = nw; l > 32; l >>= 2) { m2 = m << 1; q = m2 << 3; for (j = m; j < m2; j++) { p = ip[j] << 2; ip[m + j] = p; ip[m2 + j] = p + q; } m = m2; }}void makect(int nc, int *ip, double *c){ int j, nch; double delta; ip[1] = nc; if (nc > 1) { nch = nc >> 1; delta = atan(1.0) / nch; c[0] = cos(delta * nch); c[nch] = 0.5 * c[0]; for (j = 1; j < nch; j++) { c[j] = 0.5 * cos(delta * j); c[nc - j] = 0.5 * sin(delta * j); } }}/* -------- child routines -------- */#ifdef USE_CDFT_PTHREADS#define USE_CDFT_THREADS#ifndef CDFT_THREADS_BEGIN_N#define CDFT_THREADS_BEGIN_N 8192#endif#ifndef CDFT_4THREADS_BEGIN_N#define CDFT_4THREADS_BEGIN_N 65536#endif#include <pthread.h>#include <stdio.h>#include <stdlib.h>#define cdft_thread_t pthread_t#define cdft_thread_create(thp,func,argp) { \ if (pthread_create(thp, NULL, func, (void *) argp) != 0) { \ fprintf(stderr, "cdft thread error\n"); \ exit(1); \ } \}#define cdft_thread_wait(th) { \ if (pthread_join(th, NULL) != 0) { \ fprintf(stderr, "cdft thread error\n"); \ exit(1); \ } \}#endif /* USE_CDFT_PTHREADS */#ifdef USE_CDFT_WINTHREADS#define USE_CDFT_THREADS#ifndef CDFT_THREADS_BEGIN_N#define CDFT_THREADS_BEGIN_N 32768#endif#ifndef CDFT_4THREADS_BEGIN_N#define CDFT_4THREADS_BEGIN_N 524288#endif#include <windows.h>#include <stdio.h>#include <stdlib.h>#define cdft_thread_t HANDLE#define cdft_thread_create(thp,func,argp) { \ DWORD thid; \ *(thp) = CreateThread(NULL, 0, (LPTHREAD_START_ROUTINE) func, (LPVOID) argp, 0, &thid); \ if (*(thp) == 0) { \ fprintf(stderr, "cdft thread error\n"); \ exit(1); \ } \}#define cdft_thread_wait(th) { \ WaitForSingleObject(th, INFINITE); \ CloseHandle(th); \}#endif /* USE_CDFT_WINTHREADS */void cftfsub(int n, double *a, int *ip, int nw, double *w){ void bitrv2(int n, int *ip, double *a); void bitrv216(double *a); void bitrv208(double *a); void cftf1st(int n, double *a, double *w); void cftrec4(int n, double *a, int nw, double *w); void cftleaf(int n, int isplt, double *a, int nw, double *w); void cftfx41(int n, double *a, int nw, double *w); void cftf161(double *a, double *w); void cftf081(double *a, double *w); void cftf040(double *a); void cftx020(double *a);#ifdef USE_CDFT_THREADS void cftrec4_th(int n, double *a, int nw, double *w);#endif /* USE_CDFT_THREADS */ if (n > 8) { if (n > 32) { cftf1st(n, a, &w[nw - (n >> 2)]);#ifdef USE_CDFT_THREADS if (n > CDFT_THREADS_BEGIN_N) { cftrec4_th(n, a, nw, w); } else #endif /* USE_CDFT_THREADS */ if (n > 512) { cftrec4(n, a, nw, w); } else if (n > 128) { cftleaf(n, 1, a, nw, w); } else { cftfx41(n, a, nw, w); } bitrv2(n, ip, a); } else if (n == 32) { cftf161(a, &w[nw - 8]); bitrv216(a); } else { cftf081(a, w); bitrv208(a); } } else if (n == 8) { cftf040(a); } else if (n == 4) { cftx020(a); }}void cftbsub(int n, double *a, int *ip, int nw, double *w){ void bitrv2conj(int n, int *ip, double *a); void bitrv216neg(double *a); void bitrv208neg(double *a); void cftb1st(int n, double *a, double *w); void cftrec4(int n, double *a, int nw, double *w); void cftleaf(int n, int isplt, double *a, int nw, double *w); void cftfx41(int n, double *a, int nw, double *w); void cftf161(double *a, double *w); void cftf081(double *a, double *w); void cftb040(double *a); void cftx020(double *a);#ifdef USE_CDFT_THREADS void cftrec4_th(int n, double *a, int nw, double *w);#endif /* USE_CDFT_THREADS */ if (n > 8) { if (n > 32) { cftb1st(n, a, &w[nw - (n >> 2)]);#ifdef USE_CDFT_THREADS if (n > CDFT_THREADS_BEGIN_N) { cftrec4_th(n, a, nw, w); } else #endif /* USE_CDFT_THREADS */ if (n > 512) { cftrec4(n, a, nw, w); } else if (n > 128) { cftleaf(n, 1, a, nw, w); } else { cftfx41(n, a, nw, w); } bitrv2conj(n, ip, a); } else if (n == 32) { cftf161(a, &w[nw - 8]); bitrv216neg(a); } else { cftf081(a, w); bitrv208neg(a); } } else if (n == 8) { cftb040(a); } else if (n == 4) { cftx020(a); }}void bitrv2(int n, int *ip, double *a){ int j, j1, k, k1, l, m, nh, nm; double xr, xi, yr, yi; m = 1; for (l = n >> 2; l > 8; l >>= 2) { m <<= 1; } nh = n >> 1; nm = 4 * m; if (l == 8) { for (k = 0; k < m; k++) { for (j = 0; j < k; j++) { j1 = 4 * j + 2 * ip[m + k]; k1 = 4 * k + 2 * ip[m + j]; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 += 2 * nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 -= nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr; a[j1 + 1] = yi; a[k1] = xr; a[k1 + 1] = xi; j1 += nm; k1 += 2 * nm; xr = a[j1]; xi = a[j1 + 1]; yr = a[k1]; yi = a[k1 + 1]; a[j1] = yr;?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -