[Scipy-svn] r4283 - branches/refactor_fft/scipy/fftpack/src/djbfft
scipy-svn at scipy.org
scipy-svn at scipy.org
Mon May 12 09:08:16 EDT 2008
Author: cdavid
Date: 2008-05-12 08:08:07 -0500 (Mon, 12 May 2008)
New Revision: 4283
Modified:
branches/refactor_fft/scipy/fftpack/src/djbfft/convolve.cxx
Log:
More cleanup for djbfft backend for convolve.
Modified: branches/refactor_fft/scipy/fftpack/src/djbfft/convolve.cxx
===================================================================
--- branches/refactor_fft/scipy/fftpack/src/djbfft/convolve.cxx 2008-05-12 13:07:26 UTC (rev 4282)
+++ branches/refactor_fft/scipy/fftpack/src/djbfft/convolve.cxx 2008-05-12 13:08:07 UTC (rev 4283)
@@ -23,12 +23,72 @@
}
/**************** convolve **********************/
-extern "C"
- void convolve(int n, double *inout, double *omega, int swap_real_imag)
+static void convolve_djbfft(int n, double *inout, double *omega, int swap_real_imag)
{
int i;
double *ptr = NULL;
+
+ i = get_cache_id_ddjbfft(n);
+ ptr = caches_ddjbfft[i].ptr;
+ COPYSTD2DJB(inout, ptr, n);
switch (n) {
+#define TMPCASE(N) case N: fftr8_##N(ptr); break
+ TMPCASE(2);
+ TMPCASE(4);
+ TMPCASE(8);
+ TMPCASE(16);
+ TMPCASE(32);
+ TMPCASE(64);
+ TMPCASE(128);
+ TMPCASE(256);
+ TMPCASE(512);
+ TMPCASE(1024);
+ TMPCASE(2048);
+ TMPCASE(4096);
+ TMPCASE(8192);
+#undef TMPCASE
+ }
+ if (swap_real_imag) {
+ int n1 = n - 1;
+ double c;
+ ptr[0] *= omega[0];
+ ptr[1] *= omega[1];
+ for (i = 2; i < n1; i += 2) {
+ c = ptr[i] * omega[i];
+ ptr[i] = ptr[i + 1] * omega[i + 1];
+ ptr[i + 1] = c;
+ }
+ } else {
+ for (i = 0; i < n; ++i) {
+ ptr[i] *= omega[i];
+ }
+ }
+ switch (n) {
+#define TMPCASE(N)case N:fftr8_un##N(ptr);break
+ TMPCASE(2);
+ TMPCASE(4);
+ TMPCASE(8);
+ TMPCASE(16);
+ TMPCASE(32);
+ TMPCASE(64);
+ TMPCASE(128);
+ TMPCASE(256);
+ TMPCASE(512);
+ TMPCASE(1024);
+ TMPCASE(2048);
+ TMPCASE(4096);
+ TMPCASE(8192);
+#undef TMPCASE
+ }
+ COPYINVDJB2STD2(ptr, inout, n);
+}
+
+extern "C"
+void convolve(int n, double *inout, double *omega, int swap_real_imag)
+{
+ bool use_def = true;
+
+ switch (n) {
case 2:;
case 4:;
case 8:;
@@ -42,71 +102,83 @@
case 2048:;
case 4096:;
case 8192:
- i = get_cache_id_ddjbfft(n);
- ptr = caches_ddjbfft[i].ptr;
- COPYSTD2DJB(inout, ptr, n);
- switch (n) {
+ use_def = false;
+ }
+
+ if (!use_def) {
+ convolve_djbfft(n, inout, omega, swap_real_imag);
+ } else {
+ convolve_def(n, inout, omega, swap_real_imag);
+ }
+}
+
+/**************** convolve **********************/
+static void convolve_z_djbfft(int n, double *inout, double *omega_real,
+ double *omega_imag)
+{
+ int i;
+ double *ptr = NULL;
+
+ i = get_cache_id_ddjbfft(n);
+ ptr = caches_ddjbfft[i].ptr;
+ COPYSTD2DJB(inout, ptr, n);
+ switch (n) {
#define TMPCASE(N) case N: fftr8_##N(ptr); break
- TMPCASE(2);
- TMPCASE(4);
- TMPCASE(8);
- TMPCASE(16);
- TMPCASE(32);
- TMPCASE(64);
- TMPCASE(128);
- TMPCASE(256);
- TMPCASE(512);
- TMPCASE(1024);
- TMPCASE(2048);
- TMPCASE(4096);
- TMPCASE(8192);
+ TMPCASE(2);
+ TMPCASE(4);
+ TMPCASE(8);
+ TMPCASE(16);
+ TMPCASE(32);
+ TMPCASE(64);
+ TMPCASE(128);
+ TMPCASE(256);
+ TMPCASE(512);
+ TMPCASE(1024);
+ TMPCASE(2048);
+ TMPCASE(4096);
+ TMPCASE(8192);
#undef TMPCASE
+ }
+ {
+ int n1 = n - 1;
+ double c;
+ ptr[0] *= (omega_real[0] + omega_imag[0]);
+ ptr[1] *= (omega_real[1] + omega_imag[1]);
+ for (i = 2; i < n1; i += 2) {
+ c = ptr[i] * omega_imag[i];
+ ptr[i] *= omega_real[i];
+ ptr[i] += ptr[i + 1] * omega_imag[i + 1];
+ ptr[i + 1] *= omega_real[i + 1];
+ ptr[i + 1] += c;
}
- if (swap_real_imag) {
- int n1 = n - 1;
- double c;
- ptr[0] *= omega[0];
- ptr[1] *= omega[1];
- for (i = 2; i < n1; i += 2) {
- c = ptr[i] * omega[i];
- ptr[i] = ptr[i + 1] * omega[i + 1];
- ptr[i + 1] = c;
- }
- } else
- for (i = 0; i < n; ++i)
- ptr[i] *= omega[i];
- switch (n) {
+ }
+ switch (n) {
#define TMPCASE(N)case N:fftr8_un##N(ptr);break
- TMPCASE(2);
- TMPCASE(4);
- TMPCASE(8);
- TMPCASE(16);
- TMPCASE(32);
- TMPCASE(64);
- TMPCASE(128);
- TMPCASE(256);
- TMPCASE(512);
- TMPCASE(1024);
- TMPCASE(2048);
- TMPCASE(4096);
- TMPCASE(8192);
+ TMPCASE(2);
+ TMPCASE(4);
+ TMPCASE(8);
+ TMPCASE(16);
+ TMPCASE(32);
+ TMPCASE(64);
+ TMPCASE(128);
+ TMPCASE(256);
+ TMPCASE(512);
+ TMPCASE(1024);
+ TMPCASE(2048);
+ TMPCASE(4096);
+ TMPCASE(8192);
#undef TMPCASE
- }
- COPYINVDJB2STD2(ptr, inout, n);
- return;
}
- {
- convolve_def(n, inout, omega, swap_real_imag);
- }
+ COPYINVDJB2STD2(ptr, inout, n);
+ return;
}
-/**************** convolve **********************/
extern "C"
void convolve_z(int n, double *inout, double *omega_real,
double *omega_imag)
{
- int i;
- double *ptr = NULL;
+ bool use_def = true;
+
switch (n) {
case 2:;
case 4:;
@@ -121,69 +193,79 @@
case 2048:;
case 4096:;
case 8192:
- i = get_cache_id_ddjbfft(n);
- ptr = caches_ddjbfft[i].ptr;
- COPYSTD2DJB(inout, ptr, n);
- switch (n) {
-#define TMPCASE(N) case N: fftr8_##N(ptr); break
- TMPCASE(2);
- TMPCASE(4);
- TMPCASE(8);
- TMPCASE(16);
- TMPCASE(32);
- TMPCASE(64);
- TMPCASE(128);
- TMPCASE(256);
- TMPCASE(512);
- TMPCASE(1024);
- TMPCASE(2048);
- TMPCASE(4096);
- TMPCASE(8192);
-#undef TMPCASE
- }
- {
- int n1 = n - 1;
- double c;
- ptr[0] *= (omega_real[0] + omega_imag[0]);
- ptr[1] *= (omega_real[1] + omega_imag[1]);
- for (i = 2; i < n1; i += 2) {
- c = ptr[i] * omega_imag[i];
- ptr[i] *= omega_real[i];
- ptr[i] += ptr[i + 1] * omega_imag[i + 1];
- ptr[i + 1] *= omega_real[i + 1];
- ptr[i + 1] += c;
- }
- }
- switch (n) {
-#define TMPCASE(N)case N:fftr8_un##N(ptr);break
- TMPCASE(2);
- TMPCASE(4);
- TMPCASE(8);
- TMPCASE(16);
- TMPCASE(32);
- TMPCASE(64);
- TMPCASE(128);
- TMPCASE(256);
- TMPCASE(512);
- TMPCASE(1024);
- TMPCASE(2048);
- TMPCASE(4096);
- TMPCASE(8192);
-#undef TMPCASE
- }
- COPYINVDJB2STD2(ptr, inout, n);
- return;
+ use_def = false;
}
- {
+
+ if (use_def) {
convolve_z_def(n, inout, omega_real, omega_imag);
+ } else {
+ convolve_z_djbfft(n, inout, omega_real, omega_imag);
}
}
+static void init_convolution_kernel_djbfft(int n, double *omega, int d,
+ double (*kernel_func) (int),
+ int zero_nyquist)
+{
+ int k, n2 = n / 2;
+ unsigned int *f =
+ (unsigned int *) malloc(sizeof(int) * (n));
+ fftfreq_rtable(f, n);
+ for (k = 1; k < n; ++k)
+ if (f[k] > n2)
+ f[k] -= n;
+ omega[0] = (*kernel_func) (0) / n;
+ switch (d % 4) {
+ case 0:
+ for (k = 2; k < n - 1; k += 2) {
+ omega[k] =
+ (*kernel_func) (f[k]) / n2;
+ omega[k + 1] = -omega[k];
+ }
+ omega[1] =
+ (zero_nyquist ? 0.0
+ : (*kernel_func) (n2) / n);
+ break;
+ case 1:;
+ case -3:
+ for (k = 2; k < n - 1; k += 2)
+ omega[k] = omega[k + 1] =
+ -(*kernel_func) (f[k]) / n2;
+ omega[1] =
+ (zero_nyquist ? 0.0
+ : (*kernel_func) (n2) / n);
+ break;
+ case 2:;
+ case -2:
+ for (k = 2; k < n - 1; k += 2) {
+ omega[k] =
+ -(*kernel_func) (f[k]) / n2;
+ omega[k + 1] = -omega[k];
+ }
+ omega[1] =
+ (zero_nyquist ? 0.0 :
+ -(*kernel_func) (n2) / n);
+ break;
+ case 3:;
+ case -1:
+ for (k = 2; k < n - 1; k += 2)
+ omega[k] = omega[k + 1] =
+ (*kernel_func) (f[k]) / n2;
+ omega[1] =
+ (zero_nyquist ? 0.0 :
+ -(*kernel_func) (n2) / n);
+ break;
+ }
+ free(f);
+}
+
extern "C"
void init_convolution_kernel(int n, double *omega, int d,
double (*kernel_func) (int),
int zero_nyquist)
{
+ bool use_def = true;
+
/*
omega[k] = pow(sqrt(-1),d) * kernel_func(k)
omega[0] = kernel_func(0)
@@ -203,60 +285,14 @@
case 2048:;
case 4096:;
case 8192:
- {
- int k, n2 = n / 2;
- unsigned int *f =
- (unsigned int *) malloc(sizeof(int) * (n));
- fftfreq_rtable(f, n);
- for (k = 1; k < n; ++k)
- if (f[k] > n2)
- f[k] -= n;
- omega[0] = (*kernel_func) (0) / n;
- switch (d % 4) {
- case 0:
- for (k = 2; k < n - 1; k += 2) {
- omega[k] =
- (*kernel_func) (f[k]) / n2;
- omega[k + 1] = -omega[k];
- }
- omega[1] =
- (zero_nyquist ? 0.0
- : (*kernel_func) (n2) / n);
- break;
- case 1:;
- case -3:
- for (k = 2; k < n - 1; k += 2)
- omega[k] = omega[k + 1] =
- -(*kernel_func) (f[k]) / n2;
- omega[1] =
- (zero_nyquist ? 0.0
- : (*kernel_func) (n2) / n);
- break;
- case 2:;
- case -2:
- for (k = 2; k < n - 1; k += 2) {
- omega[k] =
- -(*kernel_func) (f[k]) / n2;
- omega[k + 1] = -omega[k];
- }
- omega[1] =
- (zero_nyquist ? 0.0 :
- -(*kernel_func) (n2) / n);
- break;
- case 3:;
- case -1:
- for (k = 2; k < n - 1; k += 2)
- omega[k] = omega[k + 1] =
- (*kernel_func) (f[k]) / n2;
- omega[1] =
- (zero_nyquist ? 0.0 :
- -(*kernel_func) (n2) / n);
- break;
- }
- free(f);
- }
- return;
+ use_def = false;
}
- init_convolution_kernel_def(n, omega, d, kernel_func,
- zero_nyquist);
+
+ if (use_def) {
+ init_convolution_kernel_def(n, omega, d, kernel_func,
+ zero_nyquist);
+ } else {
+ init_convolution_kernel_djbfft(n, omega, d, kernel_func,
+ zero_nyquist);
+ }
}
More information about the Scipy-svn
mailing list