summaryrefslogtreecommitdiff
path: root/numpy/fft
diff options
context:
space:
mode:
Diffstat (limited to 'numpy/fft')
-rw-r--r--numpy/fft/SConscript8
-rw-r--r--numpy/fft/SConstruct2
-rw-r--r--numpy/fft/__init__.py43
-rw-r--r--numpy/fft/fftpack.c1501
-rw-r--r--numpy/fft/fftpack.h28
-rw-r--r--numpy/fft/fftpack.py550
-rw-r--r--numpy/fft/fftpack_litemodule.c330
-rw-r--r--numpy/fft/helper.py116
-rw-r--r--numpy/fft/info.py29
-rw-r--r--numpy/fft/setup.py19
-rw-r--r--numpy/fft/setupscons.py15
-rw-r--r--numpy/fft/tests/test_fftpack.py23
-rw-r--r--numpy/fft/tests/test_helper.py42
13 files changed, 2706 insertions, 0 deletions
diff --git a/numpy/fft/SConscript b/numpy/fft/SConscript
new file mode 100644
index 000000000..42764229a
--- /dev/null
+++ b/numpy/fft/SConscript
@@ -0,0 +1,8 @@
+# Last Change: Thu Jun 12 06:00 PM 2008 J
+# vim:syntax=python
+from numscons import GetNumpyEnvironment
+
+env = GetNumpyEnvironment(ARGUMENTS)
+
+env.NumpyPythonExtension('fftpack_lite',
+ source = ['fftpack_litemodule.c', 'fftpack.c'])
diff --git a/numpy/fft/SConstruct b/numpy/fft/SConstruct
new file mode 100644
index 000000000..a377d8391
--- /dev/null
+++ b/numpy/fft/SConstruct
@@ -0,0 +1,2 @@
+from numscons import GetInitEnvironment
+GetInitEnvironment(ARGUMENTS).DistutilsSConscript('SConscript')
diff --git a/numpy/fft/__init__.py b/numpy/fft/__init__.py
new file mode 100644
index 000000000..1171a6026
--- /dev/null
+++ b/numpy/fft/__init__.py
@@ -0,0 +1,43 @@
+"""
+Discrete Fast Fourier Transform (FFT)
+=====================================
+
+========= =========================================================
+Standard FFTs
+===================================================================
+fft Discrete Fourier transform.
+ifft Inverse discrete Fourier transform.
+fft2 Discrete Fourier transform in two dimensions.
+ifft2 Inverse discrete Fourier transform in two dimensions.
+fftn Discrete Fourier transform in N-dimensions.
+ifftn Inverse discrete Fourier transform in N dimensions.
+========= =========================================================
+
+========= ==========================================================
+Real FFTs
+====================================================================
+rfft Real discrete Fourier transform.
+irfft Inverse real discrete Fourier transform.
+rfft2 Real discrete Fourier transform in two dimensions.
+irfft2 Inverse real discrete Fourier transform in two dimensions.
+rfftn Real discrete Fourier transform in N dimensions.
+irfftn Inverse real discrete Fourier transform in N dimensions.
+========= ==========================================================
+
+========= =========================================================
+Hermite FFTs
+===================================================================
+hfft Hermite discrete Fourier transform.
+ihfft Inverse hermite discrete Fourier transform.
+========= =========================================================
+
+"""
+# To get sub-modules
+from info import __doc__
+
+from fftpack import *
+from helper import *
+
+from numpy.testing import Tester
+test = Tester().test
+bench = Tester().bench
diff --git a/numpy/fft/fftpack.c b/numpy/fft/fftpack.c
new file mode 100644
index 000000000..9c8fd118a
--- /dev/null
+++ b/numpy/fft/fftpack.c
@@ -0,0 +1,1501 @@
+/*
+fftpack.c : A set of FFT routines in C.
+Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber (Version 4, 1985).
+
+*/
+
+/* isign is +1 for backward and -1 for forward transforms */
+
+#include <math.h>
+#include <stdio.h>
+#define DOUBLE
+
+#ifdef DOUBLE
+#define Treal double
+#else
+#define Treal float
+#endif
+
+
+#define ref(u,a) u[a]
+
+#define MAXFAC 13 /* maximum number of factors in factorization of n */
+#define NSPECIAL 4 /* number of factors for which we have special-case routines */
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+/* ----------------------------------------------------------------------
+ passf2, passf3, passf4, passf5, passf. Complex FFT passes fwd and bwd.
+---------------------------------------------------------------------- */
+
+static void passf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[], int isign)
+ /* isign==+1 for backward transform */
+ {
+ int i, k, ah, ac;
+ Treal ti2, tr2;
+ if (ido <= 2) {
+ for (k=0; k<l1; k++) {
+ ah = k*ido;
+ ac = 2*k*ido;
+ ch[ah] = ref(cc,ac) + ref(cc,ac + ido);
+ ch[ah + ido*l1] = ref(cc,ac) - ref(cc,ac + ido);
+ ch[ah+1] = ref(cc,ac+1) + ref(cc,ac + ido + 1);
+ ch[ah + ido*l1 + 1] = ref(cc,ac+1) - ref(cc,ac + ido + 1);
+ }
+ } else {
+ for (k=0; k<l1; k++) {
+ for (i=0; i<ido-1; i+=2) {
+ ah = i + k*ido;
+ ac = i + 2*k*ido;
+ ch[ah] = ref(cc,ac) + ref(cc,ac + ido);
+ tr2 = ref(cc,ac) - ref(cc,ac + ido);
+ ch[ah+1] = ref(cc,ac+1) + ref(cc,ac + 1 + ido);
+ ti2 = ref(cc,ac+1) - ref(cc,ac + 1 + ido);
+ ch[ah+l1*ido+1] = wa1[i]*ti2 + isign*wa1[i+1]*tr2;
+ ch[ah+l1*ido] = wa1[i]*tr2 - isign*wa1[i+1]*ti2;
+ }
+ }
+ }
+ } /* passf2 */
+
+
+static void passf3(int ido, int l1, const Treal cc[], Treal ch[],
+ const Treal wa1[], const Treal wa2[], int isign)
+ /* isign==+1 for backward transform */
+ {
+ static const Treal taur = -0.5;
+ static const Treal taui = 0.866025403784439;
+ int i, k, ac, ah;
+ Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
+ if (ido == 2) {
+ for (k=1; k<=l1; k++) {
+ ac = (3*k - 2)*ido;
+ tr2 = ref(cc,ac) + ref(cc,ac + ido);
+ cr2 = ref(cc,ac - ido) + taur*tr2;
+ ah = (k - 1)*ido;
+ ch[ah] = ref(cc,ac - ido) + tr2;
+
+ ti2 = ref(cc,ac + 1) + ref(cc,ac + ido + 1);
+ ci2 = ref(cc,ac - ido + 1) + taur*ti2;
+ ch[ah + 1] = ref(cc,ac - ido + 1) + ti2;
+
+ cr3 = isign*taui*(ref(cc,ac) - ref(cc,ac + ido));
+ ci3 = isign*taui*(ref(cc,ac + 1) - ref(cc,ac + ido + 1));
+ ch[ah + l1*ido] = cr2 - ci3;
+ ch[ah + 2*l1*ido] = cr2 + ci3;
+ ch[ah + l1*ido + 1] = ci2 + cr3;
+ ch[ah + 2*l1*ido + 1] = ci2 - cr3;
+ }
+ } else {
+ for (k=1; k<=l1; k++) {
+ for (i=0; i<ido-1; i+=2) {
+ ac = i + (3*k - 2)*ido;
+ tr2 = ref(cc,ac) + ref(cc,ac + ido);
+ cr2 = ref(cc,ac - ido) + taur*tr2;
+ ah = i + (k-1)*ido;
+ ch[ah] = ref(cc,ac - ido) + tr2;
+ ti2 = ref(cc,ac + 1) + ref(cc,ac + ido + 1);
+ ci2 = ref(cc,ac - ido + 1) + taur*ti2;
+ ch[ah + 1] = ref(cc,ac - ido + 1) + ti2;
+ cr3 = isign*taui*(ref(cc,ac) - ref(cc,ac + ido));
+ ci3 = isign*taui*(ref(cc,ac + 1) - ref(cc,ac + ido + 1));
+ dr2 = cr2 - ci3;
+ dr3 = cr2 + ci3;
+ di2 = ci2 + cr3;
+ di3 = ci2 - cr3;
+ ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;
+ ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;
+ ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;
+ ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;
+ }
+ }
+ }
+ } /* passf3 */
+
+
+static void passf4(int ido, int l1, const Treal cc[], Treal ch[],
+ const Treal wa1[], const Treal wa2[], const Treal wa3[], int isign)
+ /* isign == -1 for forward transform and +1 for backward transform */
+ {
+ int i, k, ac, ah;
+ Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
+ if (ido == 2) {
+ for (k=0; k<l1; k++) {
+ ac = 4*k*ido + 1;
+ ti1 = ref(cc,ac) - ref(cc,ac + 2*ido);
+ ti2 = ref(cc,ac) + ref(cc,ac + 2*ido);
+ tr4 = ref(cc,ac + 3*ido) - ref(cc,ac + ido);
+ ti3 = ref(cc,ac + ido) + ref(cc,ac + 3*ido);
+ tr1 = ref(cc,ac - 1) - ref(cc,ac + 2*ido - 1);
+ tr2 = ref(cc,ac - 1) + ref(cc,ac + 2*ido - 1);
+ ti4 = ref(cc,ac + ido - 1) - ref(cc,ac + 3*ido - 1);
+ tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 3*ido - 1);
+ ah = k*ido;
+ ch[ah] = tr2 + tr3;
+ ch[ah + 2*l1*ido] = tr2 - tr3;
+ ch[ah + 1] = ti2 + ti3;
+ ch[ah + 2*l1*ido + 1] = ti2 - ti3;
+ ch[ah + l1*ido] = tr1 + isign*tr4;
+ ch[ah + 3*l1*ido] = tr1 - isign*tr4;
+ ch[ah + l1*ido + 1] = ti1 + isign*ti4;
+ ch[ah + 3*l1*ido + 1] = ti1 - isign*ti4;
+ }
+ } else {
+ for (k=0; k<l1; k++) {
+ for (i=0; i<ido-1; i+=2) {
+ ac = i + 1 + 4*k*ido;
+ ti1 = ref(cc,ac) - ref(cc,ac + 2*ido);
+ ti2 = ref(cc,ac) + ref(cc,ac + 2*ido);
+ ti3 = ref(cc,ac + ido) + ref(cc,ac + 3*ido);
+ tr4 = ref(cc,ac + 3*ido) - ref(cc,ac + ido);
+ tr1 = ref(cc,ac - 1) - ref(cc,ac + 2*ido - 1);
+ tr2 = ref(cc,ac - 1) + ref(cc,ac + 2*ido - 1);
+ ti4 = ref(cc,ac + ido - 1) - ref(cc,ac + 3*ido - 1);
+ tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 3*ido - 1);
+ ah = i + k*ido;
+ ch[ah] = tr2 + tr3;
+ cr3 = tr2 - tr3;
+ ch[ah + 1] = ti2 + ti3;
+ ci3 = ti2 - ti3;
+ cr2 = tr1 + isign*tr4;
+ cr4 = tr1 - isign*tr4;
+ ci2 = ti1 + isign*ti4;
+ ci4 = ti1 - isign*ti4;
+ ch[ah + l1*ido] = wa1[i]*cr2 - isign*wa1[i + 1]*ci2;
+ ch[ah + l1*ido + 1] = wa1[i]*ci2 + isign*wa1[i + 1]*cr2;
+ ch[ah + 2*l1*ido] = wa2[i]*cr3 - isign*wa2[i + 1]*ci3;
+ ch[ah + 2*l1*ido + 1] = wa2[i]*ci3 + isign*wa2[i + 1]*cr3;
+ ch[ah + 3*l1*ido] = wa3[i]*cr4 -isign*wa3[i + 1]*ci4;
+ ch[ah + 3*l1*ido + 1] = wa3[i]*ci4 + isign*wa3[i + 1]*cr4;
+ }
+ }
+ }
+ } /* passf4 */
+
+
+static void passf5(int ido, int l1, const Treal cc[], Treal ch[],
+ const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[], int isign)
+ /* isign == -1 for forward transform and +1 for backward transform */
+ {
+ static const Treal tr11 = 0.309016994374947;
+ static const Treal ti11 = 0.951056516295154;
+ static const Treal tr12 = -0.809016994374947;
+ static const Treal ti12 = 0.587785252292473;
+ int i, k, ac, ah;
+ Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
+ ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
+ if (ido == 2) {
+ for (k = 1; k <= l1; ++k) {
+ ac = (5*k - 4)*ido + 1;
+ ti5 = ref(cc,ac) - ref(cc,ac + 3*ido);
+ ti2 = ref(cc,ac) + ref(cc,ac + 3*ido);
+ ti4 = ref(cc,ac + ido) - ref(cc,ac + 2*ido);
+ ti3 = ref(cc,ac + ido) + ref(cc,ac + 2*ido);
+ tr5 = ref(cc,ac - 1) - ref(cc,ac + 3*ido - 1);
+ tr2 = ref(cc,ac - 1) + ref(cc,ac + 3*ido - 1);
+ tr4 = ref(cc,ac + ido - 1) - ref(cc,ac + 2*ido - 1);
+ tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 2*ido - 1);
+ ah = (k - 1)*ido;
+ ch[ah] = ref(cc,ac - ido - 1) + tr2 + tr3;
+ ch[ah + 1] = ref(cc,ac - ido) + ti2 + ti3;
+ cr2 = ref(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3;
+ ci2 = ref(cc,ac - ido) + tr11*ti2 + tr12*ti3;
+ cr3 = ref(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3;
+ ci3 = ref(cc,ac - ido) + tr12*ti2 + tr11*ti3;
+ cr5 = isign*(ti11*tr5 + ti12*tr4);
+ ci5 = isign*(ti11*ti5 + ti12*ti4);
+ cr4 = isign*(ti12*tr5 - ti11*tr4);
+ ci4 = isign*(ti12*ti5 - ti11*ti4);
+ ch[ah + l1*ido] = cr2 - ci5;
+ ch[ah + 4*l1*ido] = cr2 + ci5;
+ ch[ah + l1*ido + 1] = ci2 + cr5;
+ ch[ah + 2*l1*ido + 1] = ci3 + cr4;
+ ch[ah + 2*l1*ido] = cr3 - ci4;
+ ch[ah + 3*l1*ido] = cr3 + ci4;
+ ch[ah + 3*l1*ido + 1] = ci3 - cr4;
+ ch[ah + 4*l1*ido + 1] = ci2 - cr5;
+ }
+ } else {
+ for (k=1; k<=l1; k++) {
+ for (i=0; i<ido-1; i+=2) {
+ ac = i + 1 + (k*5 - 4)*ido;
+ ti5 = ref(cc,ac) - ref(cc,ac + 3*ido);
+ ti2 = ref(cc,ac) + ref(cc,ac + 3*ido);
+ ti4 = ref(cc,ac + ido) - ref(cc,ac + 2*ido);
+ ti3 = ref(cc,ac + ido) + ref(cc,ac + 2*ido);
+ tr5 = ref(cc,ac - 1) - ref(cc,ac + 3*ido - 1);
+ tr2 = ref(cc,ac - 1) + ref(cc,ac + 3*ido - 1);
+ tr4 = ref(cc,ac + ido - 1) - ref(cc,ac + 2*ido - 1);
+ tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 2*ido - 1);
+ ah = i + (k - 1)*ido;
+ ch[ah] = ref(cc,ac - ido - 1) + tr2 + tr3;
+ ch[ah + 1] = ref(cc,ac - ido) + ti2 + ti3;
+ cr2 = ref(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3;
+
+ ci2 = ref(cc,ac - ido) + tr11*ti2 + tr12*ti3;
+ cr3 = ref(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3;
+
+ ci3 = ref(cc,ac - ido) + tr12*ti2 + tr11*ti3;
+ cr5 = isign*(ti11*tr5 + ti12*tr4);
+ ci5 = isign*(ti11*ti5 + ti12*ti4);
+ cr4 = isign*(ti12*tr5 - ti11*tr4);
+ ci4 = isign*(ti12*ti5 - ti11*ti4);
+ dr3 = cr3 - ci4;
+ dr4 = cr3 + ci4;
+ di3 = ci3 + cr4;
+ di4 = ci3 - cr4;
+ dr5 = cr2 + ci5;
+ dr2 = cr2 - ci5;
+ di5 = ci2 - cr5;
+ di2 = ci2 + cr5;
+ ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;
+ ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;
+ ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;
+ ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;
+ ch[ah + 3*l1*ido] = wa3[i]*dr4 - isign*wa3[i+1]*di4;
+ ch[ah + 3*l1*ido + 1] = wa3[i]*di4 + isign*wa3[i+1]*dr4;
+ ch[ah + 4*l1*ido] = wa4[i]*dr5 - isign*wa4[i+1]*di5;
+ ch[ah + 4*l1*ido + 1] = wa4[i]*di5 + isign*wa4[i+1]*dr5;
+ }
+ }
+ }
+ } /* passf5 */
+
+
+static void passf(int *nac, int ido, int ip, int l1, int idl1,
+ Treal cc[], Treal ch[],
+ const Treal wa[], int isign)
+ /* isign is -1 for forward transform and +1 for backward transform */
+ {
+ int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, idj, idl, inc,idp;
+ Treal wai, war;
+
+ idot = ido / 2;
+ /* nt = ip*idl1;*/
+ ipph = (ip + 1) / 2;
+ idp = ip*ido;
+ if (ido >= l1) {
+ for (j=1; j<ipph; j++) {
+ jc = ip - j;
+ for (k=0; k<l1; k++) {
+ for (i=0; i<ido; i++) {
+ ch[i + (k + j*l1)*ido] =
+ ref(cc,i + (j + k*ip)*ido) + ref(cc,i + (jc + k*ip)*ido);
+ ch[i + (k + jc*l1)*ido] =
+ ref(cc,i + (j + k*ip)*ido) - ref(cc,i + (jc + k*ip)*ido);
+ }
+ }
+ }
+ for (k=0; k<l1; k++)
+ for (i=0; i<ido; i++)
+ ch[i + k*ido] = ref(cc,i + k*ip*ido);
+ } else {
+ for (j=1; j<ipph; j++) {
+ jc = ip - j;
+ for (i=0; i<ido; i++) {
+ for (k=0; k<l1; k++) {
+ ch[i + (k + j*l1)*ido] = ref(cc,i + (j + k*ip)*ido) + ref(cc,i + (jc + k*
+ ip)*ido);
+ ch[i + (k + jc*l1)*ido] = ref(cc,i + (j + k*ip)*ido) - ref(cc,i + (jc + k*
+ ip)*ido);
+ }
+ }
+ }
+ for (i=0; i<ido; i++)
+ for (k=0; k<l1; k++)
+ ch[i + k*ido] = ref(cc,i + k*ip*ido);
+ }
+
+ idl = 2 - ido;
+ inc = 0;
+ for (l=1; l<ipph; l++) {
+ lc = ip - l;
+ idl += ido;
+ for (ik=0; ik<idl1; ik++) {
+ cc[ik + l*idl1] = ch[ik] + wa[idl - 2]*ch[ik + idl1];
+ cc[ik + lc*idl1] = isign*wa[idl-1]*ch[ik + (ip-1)*idl1];
+ }
+ idlj = idl;
+ inc += ido;
+ for (j=2; j<ipph; j++) {
+ jc = ip - j;
+ idlj += inc;
+ if (idlj > idp) idlj -= idp;
+ war = wa[idlj - 2];
+ wai = wa[idlj-1];
+ for (ik=0; ik<idl1; ik++) {
+ cc[ik + l*idl1] += war*ch[ik + j*idl1];
+ cc[ik + lc*idl1] += isign*wai*ch[ik + jc*idl1];
+ }
+ }
+ }
+ for (j=1; j<ipph; j++)
+ for (ik=0; ik<idl1; ik++)
+ ch[ik] += ch[ik + j*idl1];
+ for (j=1; j<ipph; j++) {
+ jc = ip - j;
+ for (ik=1; ik<idl1; ik+=2) {
+ ch[ik - 1 + j*idl1] = cc[ik - 1 + j*idl1] - cc[ik + jc*idl1];
+ ch[ik - 1 + jc*idl1] = cc[ik - 1 + j*idl1] + cc[ik + jc*idl1];
+ ch[ik + j*idl1] = cc[ik + j*idl1] + cc[ik - 1 + jc*idl1];
+ ch[ik + jc*idl1] = cc[ik + j*idl1] - cc[ik - 1 + jc*idl1];
+ }
+ }
+ *nac = 1;
+ if (ido == 2) return;
+ *nac = 0;
+ for (ik=0; ik<idl1; ik++)
+ cc[ik] = ch[ik];
+ for (j=1; j<ip; j++) {
+ for (k=0; k<l1; k++) {
+ cc[(k + j*l1)*ido + 0] = ch[(k + j*l1)*ido + 0];
+ cc[(k + j*l1)*ido + 1] = ch[(k + j*l1)*ido + 1];
+ }
+ }
+ if (idot <= l1) {
+ idij = 0;
+ for (j=1; j<ip; j++) {
+ idij += 2;
+ for (i=3; i<ido; i+=2) {
+ idij += 2;
+ for (k=0; k<l1; k++) {
+ cc[i - 1 + (k + j*l1)*ido] =
+ wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] -
+ isign*wa[idij-1]*ch[i + (k + j*l1)*ido];
+ cc[i + (k + j*l1)*ido] =
+ wa[idij - 2]*ch[i + (k + j*l1)*ido] +
+ isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido];
+ }
+ }
+ }
+ } else {
+ idj = 2 - ido;
+ for (j=1; j<ip; j++) {
+ idj += ido;
+ for (k = 0; k < l1; k++) {
+ idij = idj;
+ for (i=3; i<ido; i+=2) {
+ idij += 2;
+ cc[i - 1 + (k + j*l1)*ido] =
+ wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] -
+ isign*wa[idij-1]*ch[i + (k + j*l1)*ido];
+ cc[i + (k + j*l1)*ido] =
+ wa[idij - 2]*ch[i + (k + j*l1)*ido] +
+ isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido];
+ }
+ }
+ }
+ }
+ } /* passf */
+
+
+ /* ----------------------------------------------------------------------
+radf2,radb2, radf3,radb3, radf4,radb4, radf5,radb5, radfg,radbg.
+Treal FFT passes fwd and bwd.
+---------------------------------------------------------------------- */
+
+static void radf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[])
+ {
+ int i, k, ic;
+ Treal ti2, tr2;
+ for (k=0; k<l1; k++) {
+ ch[2*k*ido] =
+ ref(cc,k*ido) + ref(cc,(k + l1)*ido);
+ ch[(2*k+1)*ido + ido-1] =
+ ref(cc,k*ido) - ref(cc,(k + l1)*ido);
+ }
+ if (ido < 2) return;
+ if (ido != 2) {
+ for (k=0; k<l1; k++) {
+ for (i=2; i<ido; i+=2) {
+ ic = ido - i;
+ tr2 = wa1[i - 2]*ref(cc, i-1 + (k + l1)*ido) + wa1[i - 1]*ref(cc, i + (k + l1)*ido);
+ ti2 = wa1[i - 2]*ref(cc, i + (k + l1)*ido) - wa1[i - 1]*ref(cc, i-1 + (k + l1)*ido);
+ ch[i + 2*k*ido] = ref(cc,i + k*ido) + ti2;
+ ch[ic + (2*k+1)*ido] = ti2 - ref(cc,i + k*ido);
+ ch[i - 1 + 2*k*ido] = ref(cc,i - 1 + k*ido) + tr2;
+ ch[ic - 1 + (2*k+1)*ido] = ref(cc,i - 1 + k*ido) - tr2;
+ }
+ }
+ if (ido % 2 == 1) return;
+ }
+ for (k=0; k<l1; k++) {
+ ch[(2*k+1)*ido] = -ref(cc,ido-1 + (k + l1)*ido);
+ ch[ido-1 + 2*k*ido] = ref(cc,ido-1 + k*ido);
+ }
+ } /* radf2 */
+
+
+static void radb2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[])
+ {
+ int i, k, ic;
+ Treal ti2, tr2;
+ for (k=0; k<l1; k++) {
+ ch[k*ido] =
+ ref(cc,2*k*ido) + ref(cc,ido-1 + (2*k+1)*ido);
+ ch[(k + l1)*ido] =
+ ref(cc,2*k*ido) - ref(cc,ido-1 + (2*k+1)*ido);
+ }
+ if (ido < 2) return;
+ if (ido != 2) {
+ for (k = 0; k < l1; ++k) {
+ for (i = 2; i < ido; i += 2) {
+ ic = ido - i;
+ ch[i-1 + k*ido] =
+ ref(cc,i-1 + 2*k*ido) + ref(cc,ic-1 + (2*k+1)*ido);
+ tr2 = ref(cc,i-1 + 2*k*ido) - ref(cc,ic-1 + (2*k+1)*ido);
+ ch[i + k*ido] =
+ ref(cc,i + 2*k*ido) - ref(cc,ic + (2*k+1)*ido);
+ ti2 = ref(cc,i + (2*k)*ido) + ref(cc,ic + (2*k+1)*ido);
+ ch[i-1 + (k + l1)*ido] =
+ wa1[i - 2]*tr2 - wa1[i - 1]*ti2;
+ ch[i + (k + l1)*ido] =
+ wa1[i - 2]*ti2 + wa1[i - 1]*tr2;
+ }
+ }
+ if (ido % 2 == 1) return;
+ }
+ for (k = 0; k < l1; k++) {
+ ch[ido-1 + k*ido] = 2*ref(cc,ido-1 + 2*k*ido);
+ ch[ido-1 + (k + l1)*ido] = -2*ref(cc,(2*k+1)*ido);
+ }
+ } /* radb2 */
+
+
+static void radf3(int ido, int l1, const Treal cc[], Treal ch[],
+ const Treal wa1[], const Treal wa2[])
+ {
+ static const Treal taur = -0.5;
+ static const Treal taui = 0.866025403784439;
+ int i, k, ic;
+ Treal ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3;
+ for (k=0; k<l1; k++) {
+ cr2 = ref(cc,(k + l1)*ido) + ref(cc,(k + 2*l1)*ido);
+ ch[3*k*ido] = ref(cc,k*ido) + cr2;
+ ch[(3*k+2)*ido] = taui*(ref(cc,(k + l1*2)*ido) - ref(cc,(k + l1)*ido));
+ ch[ido-1 + (3*k + 1)*ido] = ref(cc,k*ido) + taur*cr2;
+ }
+ if (ido == 1) return;
+ for (k=0; k<l1; k++) {
+ for (i=2; i<ido; i+=2) {
+ ic = ido - i;
+ dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) +
+ wa1[i - 1]*ref(cc,i + (k + l1)*ido);
+ di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
+ dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + l1*2)*ido) + wa2[i - 1]*ref(cc,i + (k + l1*2)*ido);
+ di3 = wa2[i - 2]*ref(cc,i + (k + l1*2)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + l1*2)*ido);
+ cr2 = dr2 + dr3;
+ ci2 = di2 + di3;
+ ch[i - 1 + 3*k*ido] = ref(cc,i - 1 + k*ido) + cr2;
+ ch[i + 3*k*ido] = ref(cc,i + k*ido) + ci2;
+ tr2 = ref(cc,i - 1 + k*ido) + taur*cr2;
+ ti2 = ref(cc,i + k*ido) + taur*ci2;
+ tr3 = taui*(di2 - di3);
+ ti3 = taui*(dr3 - dr2);
+ ch[i - 1 + (3*k + 2)*ido] = tr2 + tr3;
+ ch[ic - 1 + (3*k + 1)*ido] = tr2 - tr3;
+ ch[i + (3*k + 2)*ido] = ti2 + ti3;
+ ch[ic + (3*k + 1)*ido] = ti3 - ti2;
+ }
+ }
+ } /* radf3 */
+
+
+static void radb3(int ido, int l1, const Treal cc[], Treal ch[],
+ const Treal wa1[], const Treal wa2[])
+ {
+ static const Treal taur = -0.5;
+ static const Treal taui = 0.866025403784439;
+ int i, k, ic;
+ Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
+ for (k=0; k<l1; k++) {
+ tr2 = 2*ref(cc,ido-1 + (3*k + 1)*ido);
+ cr2 = ref(cc,3*k*ido) + taur*tr2;
+ ch[k*ido] = ref(cc,3*k*ido) + tr2;
+ ci3 = 2*taui*ref(cc,(3*k + 2)*ido);
+ ch[(k + l1)*ido] = cr2 - ci3;
+ ch[(k + 2*l1)*ido] = cr2 + ci3;
+ }
+ if (ido == 1) return;
+ for (k=0; k<l1; k++) {
+ for (i=2; i<ido; i+=2) {
+ ic = ido - i;
+ tr2 = ref(cc,i - 1 + (3*k + 2)*ido) + ref(cc,ic - 1 + (3*k + 1)*ido);
+ cr2 = ref(cc,i - 1 + 3*k*ido) + taur*tr2;
+ ch[i - 1 + k*ido] = ref(cc,i - 1 + 3*k*ido) + tr2;
+ ti2 = ref(cc,i + (3*k + 2)*ido) - ref(cc,ic + (3*k + 1)*ido);
+ ci2 = ref(cc,i + 3*k*ido) + taur*ti2;
+ ch[i + k*ido] = ref(cc,i + 3*k*ido) + ti2;
+ cr3 = taui*(ref(cc,i - 1 + (3*k + 2)*ido) - ref(cc,ic - 1 + (3*k + 1)*ido));
+ ci3 = taui*(ref(cc,i + (3*k + 2)*ido) + ref(cc,ic + (3*k + 1)*ido));
+ dr2 = cr2 - ci3;
+ dr3 = cr2 + ci3;
+ di2 = ci2 + cr3;
+ di3 = ci2 - cr3;
+ ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2;
+ ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2;
+ ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3;
+ ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3;
+ }
+ }
+ } /* radb3 */
+
+
+static void radf4(int ido, int l1, const Treal cc[], Treal ch[],
+ const Treal wa1[], const Treal wa2[], const Treal wa3[])
+ {
+ static const Treal hsqt2 = 0.7071067811865475;
+ int i, k, ic;
+ Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
+ for (k=0; k<l1; k++) {
+ tr1 = ref(cc,(k + l1)*ido) + ref(cc,(k + 3*l1)*ido);
+ tr2 = ref(cc,k*ido) + ref(cc,(k + 2*l1)*ido);
+ ch[4*k*ido] = tr1 + tr2;
+ ch[ido-1 + (4*k + 3)*ido] = tr2 - tr1;
+ ch[ido-1 + (4*k + 1)*ido] = ref(cc,k*ido) - ref(cc,(k + 2*l1)*ido);
+ ch[(4*k + 2)*ido] = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + l1)*ido);
+ }
+ if (ido < 2) return;
+ if (ido != 2) {
+ for (k=0; k<l1; k++) {
+ for (i=2; i<ido; i += 2) {
+ ic = ido - i;
+ cr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido);
+ ci2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
+ cr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*
+ ido);
+ ci3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*
+ ido);
+ cr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*
+ ido);
+ ci4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*
+ ido);
+ tr1 = cr2 + cr4;
+ tr4 = cr4 - cr2;
+ ti1 = ci2 + ci4;
+ ti4 = ci2 - ci4;
+ ti2 = ref(cc,i + k*ido) + ci3;
+ ti3 = ref(cc,i + k*ido) - ci3;
+ tr2 = ref(cc,i - 1 + k*ido) + cr3;
+ tr3 = ref(cc,i - 1 + k*ido) - cr3;
+ ch[i - 1 + 4*k*ido] = tr1 + tr2;
+ ch[ic - 1 + (4*k + 3)*ido] = tr2 - tr1;
+ ch[i + 4*k*ido] = ti1 + ti2;
+ ch[ic + (4*k + 3)*ido] = ti1 - ti2;
+ ch[i - 1 + (4*k + 2)*ido] = ti4 + tr3;
+ ch[ic - 1 + (4*k + 1)*ido] = tr3 - ti4;
+ ch[i + (4*k + 2)*ido] = tr4 + ti3;
+ ch[ic + (4*k + 1)*ido] = tr4 - ti3;
+ }
+ }
+ if (ido % 2 == 1) return;
+ }
+ for (k=0; k<l1; k++) {
+ ti1 = -hsqt2*(ref(cc,ido-1 + (k + l1)*ido) + ref(cc,ido-1 + (k + 3*l1)*ido));
+ tr1 = hsqt2*(ref(cc,ido-1 + (k + l1)*ido) - ref(cc,ido-1 + (k + 3*l1)*ido));
+ ch[ido-1 + 4*k*ido] = tr1 + ref(cc,ido-1 + k*ido);
+ ch[ido-1 + (4*k + 2)*ido] = ref(cc,ido-1 + k*ido) - tr1;
+ ch[(4*k + 1)*ido] = ti1 - ref(cc,ido-1 + (k + 2*l1)*ido);
+ ch[(4*k + 3)*ido] = ti1 + ref(cc,ido-1 + (k + 2*l1)*ido);
+ }
+ } /* radf4 */
+
+
+static void radb4(int ido, int l1, const Treal cc[], Treal ch[],
+ const Treal wa1[], const Treal wa2[], const Treal wa3[])
+ {
+ static const Treal sqrt2 = 1.414213562373095;
+ int i, k, ic;
+ Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
+ for (k = 0; k < l1; k++) {
+ tr1 = ref(cc,4*k*ido) - ref(cc,ido-1 + (4*k + 3)*ido);
+ tr2 = ref(cc,4*k*ido) + ref(cc,ido-1 + (4*k + 3)*ido);
+ tr3 = ref(cc,ido-1 + (4*k + 1)*ido) + ref(cc,ido-1 + (4*k + 1)*ido);
+ tr4 = ref(cc,(4*k + 2)*ido) + ref(cc,(4*k + 2)*ido);
+ ch[k*ido] = tr2 + tr3;
+ ch[(k + l1)*ido] = tr1 - tr4;
+ ch[(k + 2*l1)*ido] = tr2 - tr3;
+ ch[(k + 3*l1)*ido] = tr1 + tr4;
+ }
+ if (ido < 2) return;
+ if (ido != 2) {
+ for (k = 0; k < l1; ++k) {
+ for (i = 2; i < ido; i += 2) {
+ ic = ido - i;
+ ti1 = ref(cc,i + 4*k*ido) + ref(cc,ic + (4*k + 3)*ido);
+ ti2 = ref(cc,i + 4*k*ido) - ref(cc,ic + (4*k + 3)*ido);
+ ti3 = ref(cc,i + (4*k + 2)*ido) - ref(cc,ic + (4*k + 1)*ido);
+ tr4 = ref(cc,i + (4*k + 2)*ido) + ref(cc,ic + (4*k + 1)*ido);
+ tr1 = ref(cc,i - 1 + 4*k*ido) - ref(cc,ic - 1 + (4*k + 3)*ido);
+ tr2 = ref(cc,i - 1 + 4*k*ido) + ref(cc,ic - 1 + (4*k + 3)*ido);
+ ti4 = ref(cc,i - 1 + (4*k + 2)*ido) - ref(cc,ic - 1 + (4*k + 1)*ido);
+ tr3 = ref(cc,i - 1 + (4*k + 2)*ido) + ref(cc,ic - 1 + (4*k + 1)*ido);
+ ch[i - 1 + k*ido] = tr2 + tr3;
+ cr3 = tr2 - tr3;
+ ch[i + k*ido] = ti2 + ti3;
+ ci3 = ti2 - ti3;
+ cr2 = tr1 - tr4;
+ cr4 = tr1 + tr4;
+ ci2 = ti1 + ti4;
+ ci4 = ti1 - ti4;
+ ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*cr2 - wa1[i - 1]*ci2;
+ ch[i + (k + l1)*ido] = wa1[i - 2]*ci2 + wa1[i - 1]*cr2;
+ ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*cr3 - wa2[i - 1]*ci3;
+ ch[i + (k + 2*l1)*ido] = wa2[i - 2]*ci3 + wa2[i - 1]*cr3;
+ ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*cr4 - wa3[i - 1]*ci4;
+ ch[i + (k + 3*l1)*ido] = wa3[i - 2]*ci4 + wa3[i - 1]*cr4;
+ }
+ }
+ if (ido % 2 == 1) return;
+ }
+ for (k = 0; k < l1; k++) {
+ ti1 = ref(cc,(4*k + 1)*ido) + ref(cc,(4*k + 3)*ido);
+ ti2 = ref(cc,(4*k + 3)*ido) - ref(cc,(4*k + 1)*ido);
+ tr1 = ref(cc,ido-1 + 4*k*ido) - ref(cc,ido-1 + (4*k + 2)*ido);
+ tr2 = ref(cc,ido-1 + 4*k*ido) + ref(cc,ido-1 + (4*k + 2)*ido);
+ ch[ido-1 + k*ido] = tr2 + tr2;
+ ch[ido-1 + (k + l1)*ido] = sqrt2*(tr1 - ti1);
+ ch[ido-1 + (k + 2*l1)*ido] = ti2 + ti2;
+ ch[ido-1 + (k + 3*l1)*ido] = -sqrt2*(tr1 + ti1);
+ }
+ } /* radb4 */
+
+
+static void radf5(int ido, int l1, const Treal cc[], Treal ch[],
+ const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[])
+ {
+ static const Treal tr11 = 0.309016994374947;
+ static const Treal ti11 = 0.951056516295154;
+ static const Treal tr12 = -0.809016994374947;
+ static const Treal ti12 = 0.587785252292473;
+ int i, k, ic;
+ Treal ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3, dr4, dr5,
+ cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;
+ for (k = 0; k < l1; k++) {
+ cr2 = ref(cc,(k + 4*l1)*ido) + ref(cc,(k + l1)*ido);
+ ci5 = ref(cc,(k + 4*l1)*ido) - ref(cc,(k + l1)*ido);
+ cr3 = ref(cc,(k + 3*l1)*ido) + ref(cc,(k + 2*l1)*ido);
+ ci4 = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + 2*l1)*ido);
+ ch[5*k*ido] = ref(cc,k*ido) + cr2 + cr3;
+ ch[ido-1 + (5*k + 1)*ido] = ref(cc,k*ido) + tr11*cr2 + tr12*cr3;
+ ch[(5*k + 2)*ido] = ti11*ci5 + ti12*ci4;
+ ch[ido-1 + (5*k + 3)*ido] = ref(cc,k*ido) + tr12*cr2 + tr11*cr3;
+ ch[(5*k + 4)*ido] = ti12*ci5 - ti11*ci4;
+ }
+ if (ido == 1) return;
+ for (k = 0; k < l1; ++k) {
+ for (i = 2; i < ido; i += 2) {
+ ic = ido - i;
+ dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido);
+ di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
+ dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*ido);
+ di3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*ido);
+ dr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*ido);
+ di4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*ido);
+ dr5 = wa4[i - 2]*ref(cc,i - 1 + (k + 4*l1)*ido) + wa4[i - 1]*ref(cc,i + (k + 4*l1)*ido);
+ di5 = wa4[i - 2]*ref(cc,i + (k + 4*l1)*ido) - wa4[i - 1]*ref(cc,i - 1 + (k + 4*l1)*ido);
+ cr2 = dr2 + dr5;
+ ci5 = dr5 - dr2;
+ cr5 = di2 - di5;
+ ci2 = di2 + di5;
+ cr3 = dr3 + dr4;
+ ci4 = dr4 - dr3;
+ cr4 = di3 - di4;
+ ci3 = di3 + di4;
+ ch[i - 1 + 5*k*ido] = ref(cc,i - 1 + k*ido) + cr2 + cr3;
+ ch[i + 5*k*ido] = ref(cc,i + k*ido) + ci2 + ci3;
+ tr2 = ref(cc,i - 1 + k*ido) + tr11*cr2 + tr12*cr3;
+ ti2 = ref(cc,i + k*ido) + tr11*ci2 + tr12*ci3;
+ tr3 = ref(cc,i - 1 + k*ido) + tr12*cr2 + tr11*cr3;
+ ti3 = ref(cc,i + k*ido) + tr12*ci2 + tr11*ci3;
+ tr5 = ti11*cr5 + ti12*cr4;
+ ti5 = ti11*ci5 + ti12*ci4;
+ tr4 = ti12*cr5 - ti11*cr4;
+ ti4 = ti12*ci5 - ti11*ci4;
+ ch[i - 1 + (5*k + 2)*ido] = tr2 + tr5;
+ ch[ic - 1 + (5*k + 1)*ido] = tr2 - tr5;
+ ch[i + (5*k + 2)*ido] = ti2 + ti5;
+ ch[ic + (5*k + 1)*ido] = ti5 - ti2;
+ ch[i - 1 + (5*k + 4)*ido] = tr3 + tr4;
+ ch[ic - 1 + (5*k + 3)*ido] = tr3 - tr4;
+ ch[i + (5*k + 4)*ido] = ti3 + ti4;
+ ch[ic + (5*k + 3)*ido] = ti4 - ti3;
+ }
+ }
+ } /* radf5 */
+
+
+static void radb5(int ido, int l1, const Treal cc[], Treal ch[],
+ const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[])
+ {
+ static const Treal tr11 = 0.309016994374947;
+ static const Treal ti11 = 0.951056516295154;
+ static const Treal tr12 = -0.809016994374947;
+ static const Treal ti12 = 0.587785252292473;
+ int i, k, ic;
+ Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
+ ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
+ for (k = 0; k < l1; k++) {
+ ti5 = 2*ref(cc,(5*k + 2)*ido);
+ ti4 = 2*ref(cc,(5*k + 4)*ido);
+ tr2 = 2*ref(cc,ido-1 + (5*k + 1)*ido);
+ tr3 = 2*ref(cc,ido-1 + (5*k + 3)*ido);
+ ch[k*ido] = ref(cc,5*k*ido) + tr2 + tr3;
+ cr2 = ref(cc,5*k*ido) + tr11*tr2 + tr12*tr3;
+ cr3 = ref(cc,5*k*ido) + tr12*tr2 + tr11*tr3;
+ ci5 = ti11*ti5 + ti12*ti4;
+ ci4 = ti12*ti5 - ti11*ti4;
+ ch[(k + l1)*ido] = cr2 - ci5;
+ ch[(k + 2*l1)*ido] = cr3 - ci4;
+ ch[(k + 3*l1)*ido] = cr3 + ci4;
+ ch[(k + 4*l1)*ido] = cr2 + ci5;
+ }
+ if (ido == 1) return;
+ for (k = 0; k < l1; ++k) {
+ for (i = 2; i < ido; i += 2) {
+ ic = ido - i;
+ ti5 = ref(cc,i + (5*k + 2)*ido) + ref(cc,ic + (5*k + 1)*ido);
+ ti2 = ref(cc,i + (5*k + 2)*ido) - ref(cc,ic + (5*k + 1)*ido);
+ ti4 = ref(cc,i + (5*k + 4)*ido) + ref(cc,ic + (5*k + 3)*ido);
+ ti3 = ref(cc,i + (5*k + 4)*ido) - ref(cc,ic + (5*k + 3)*ido);
+ tr5 = ref(cc,i - 1 + (5*k + 2)*ido) - ref(cc,ic - 1 + (5*k + 1)*ido);
+ tr2 = ref(cc,i - 1 + (5*k + 2)*ido) + ref(cc,ic - 1 + (5*k + 1)*ido);
+ tr4 = ref(cc,i - 1 + (5*k + 4)*ido) - ref(cc,ic - 1 + (5*k + 3)*ido);
+ tr3 = ref(cc,i - 1 + (5*k + 4)*ido) + ref(cc,ic - 1 + (5*k + 3)*ido);
+ ch[i - 1 + k*ido] = ref(cc,i - 1 + 5*k*ido) + tr2 + tr3;
+ ch[i + k*ido] = ref(cc,i + 5*k*ido) + ti2 + ti3;
+ cr2 = ref(cc,i - 1 + 5*k*ido) + tr11*tr2 + tr12*tr3;
+
+ ci2 = ref(cc,i + 5*k*ido) + tr11*ti2 + tr12*ti3;
+ cr3 = ref(cc,i - 1 + 5*k*ido) + tr12*tr2 + tr11*tr3;
+
+ ci3 = ref(cc,i + 5*k*ido) + tr12*ti2 + tr11*ti3;
+ cr5 = ti11*tr5 + ti12*tr4;
+ ci5 = ti11*ti5 + ti12*ti4;
+ cr4 = ti12*tr5 - ti11*tr4;
+ ci4 = ti12*ti5 - ti11*ti4;
+ dr3 = cr3 - ci4;
+ dr4 = cr3 + ci4;
+ di3 = ci3 + cr4;
+ di4 = ci3 - cr4;
+ dr5 = cr2 + ci5;
+ dr2 = cr2 - ci5;
+ di5 = ci2 - cr5;
+ di2 = ci2 + cr5;
+ ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2;
+ ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2;
+ ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3;
+ ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3;
+ ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*dr4 - wa3[i - 1]*di4;
+ ch[i + (k + 3*l1)*ido] = wa3[i - 2]*di4 + wa3[i - 1]*dr4;
+ ch[i - 1 + (k + 4*l1)*ido] = wa4[i - 2]*dr5 - wa4[i - 1]*di5;
+ ch[i + (k + 4*l1)*ido] = wa4[i - 2]*di5 + wa4[i - 1]*dr5;
+ }
+ }
+ } /* radb5 */
+
+
+static void radfg(int ido, int ip, int l1, int idl1,
+ Treal cc[], Treal ch[], const Treal wa[])
+ {
+ static const Treal twopi = 6.28318530717959;
+ int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is, nbd;
+ Treal dc2, ai1, ai2, ar1, ar2, ds2, dcp, arg, dsp, ar1h, ar2h;
+ arg = twopi / ip;
+ dcp = cos(arg);
+ dsp = sin(arg);
+ ipph = (ip + 1) / 2;
+ nbd = (ido - 1) / 2;
+ if (ido != 1) {
+ for (ik=0; ik<idl1; ik++) ch[ik] = cc[ik];
+ for (j=1; j<ip; j++)
+ for (k=0; k<l1; k++)
+ ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido];
+ if (nbd <= l1) {
+ is = -ido;
+ for (j=1; j<ip; j++) {
+ is += ido;
+ idij = is-1;
+ for (i=2; i<ido; i+=2) {
+ idij += 2;
+ for (k=0; k<l1; k++) {
+ ch[i - 1 + (k + j*l1)*ido] =
+ wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido];
+ ch[i + (k + j*l1)*ido] =
+ wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido];
+ }
+ }
+ }
+ } else {
+ is = -ido;
+ for (j=1; j<ip; j++) {
+ is += ido;
+ for (k=0; k<l1; k++) {
+ idij = is-1;
+ for (i=2; i<ido; i+=2) {
+ idij += 2;
+ ch[i - 1 + (k + j*l1)*ido] =
+ wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido];
+ ch[i + (k + j*l1)*ido] =
+ wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido];
+ }
+ }
+ }
+ }
+ if (nbd >= l1) {
+ for (j=1; j<ipph; j++) {
+ jc = ip - j;
+ for (k=0; k<l1; k++) {
+ for (i=2; i<ido; i+=2) {
+ cc[i - 1 + (k + j*l1)*ido] = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
+ cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido];
+ cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
+ cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido];
+ }
+ }
+ }
+ } else {
+ for (j=1; j<ipph; j++) {
+ jc = ip - j;
+ for (i=2; i<ido; i+=2) {
+ for (k=0; k<l1; k++) {
+ cc[i - 1 + (k + j*l1)*ido] =
+ ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
+ cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido];
+ cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
+ cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido];
+ }
+ }
+ }
+ }
+ } else { /* now ido == 1 */
+ for (ik=0; ik<idl1; ik++) cc[ik] = ch[ik];
+ }
+ for (j=1; j<ipph; j++) {
+ jc = ip - j;
+ for (k=0; k<l1; k++) {
+ cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido] + ch[(k + jc*l1)*ido];
+ cc[(k + jc*l1)*ido] = ch[(k + jc*l1)*ido] - ch[(k + j*l1)*ido];
+ }
+ }
+
+ ar1 = 1;
+ ai1 = 0;
+ for (l=1; l<ipph; l++) {
+ lc = ip - l;
+ ar1h = dcp*ar1 - dsp*ai1;
+ ai1 = dcp*ai1 + dsp*ar1;
+ ar1 = ar1h;
+ for (ik=0; ik<idl1; ik++) {
+ ch[ik + l*idl1] = cc[ik] + ar1*cc[ik + idl1];
+ ch[ik + lc*idl1] = ai1*cc[ik + (ip-1)*idl1];
+ }
+ dc2 = ar1;
+ ds2 = ai1;
+ ar2 = ar1;
+ ai2 = ai1;
+ for (j=2; j<ipph; j++) {
+ jc = ip - j;
+ ar2h = dc2*ar2 - ds2*ai2;
+ ai2 = dc2*ai2 + ds2*ar2;
+ ar2 = ar2h;
+ for (ik=0; ik<idl1; ik++) {
+ ch[ik + l*idl1] += ar2*cc[ik + j*idl1];
+ ch[ik + lc*idl1] += ai2*cc[ik + jc*idl1];
+ }
+ }
+ }
+ for (j=1; j<ipph; j++)
+ for (ik=0; ik<idl1; ik++)
+ ch[ik] += cc[ik + j*idl1];
+
+ if (ido >= l1) {
+ for (k=0; k<l1; k++) {
+ for (i=0; i<ido; i++) {
+ ref(cc,i + k*ip*ido) = ch[i + k*ido];
+ }
+ }
+ } else {
+ for (i=0; i<ido; i++) {
+ for (k=0; k<l1; k++) {
+ ref(cc,i + k*ip*ido) = ch[i + k*ido];
+ }
+ }
+ }
+ for (j=1; j<ipph; j++) {
+ jc = ip - j;
+ j2 = 2*j;
+ for (k=0; k<l1; k++) {
+ ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) =
+ ch[(k + j*l1)*ido];
+ ref(cc,(j2 + k*ip)*ido) =
+ ch[(k + jc*l1)*ido];
+ }
+ }
+ if (ido == 1) return;
+ if (nbd >= l1) {
+ for (j=1; j<ipph; j++) {
+ jc = ip - j;
+ j2 = 2*j;
+ for (k=0; k<l1; k++) {
+ for (i=2; i<ido; i+=2) {
+ ic = ido - i;
+ ref(cc,i - 1 + (j2 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
+ ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido];
+ ref(cc,i + (j2 + k*ip)*ido) = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
+ ref(cc,ic + (j2 - 1 + k*ip)*ido) = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido];
+ }
+ }
+ }
+ } else {
+ for (j=1; j<ipph; j++) {
+ jc = ip - j;
+ j2 = 2*j;
+ for (i=2; i<ido; i+=2) {
+ ic = ido - i;
+ for (k=0; k<l1; k++) {
+ ref(cc,i - 1 + (j2 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
+ ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido];
+ ref(cc,i + (j2 + k*ip)*ido) = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
+ ref(cc,ic + (j2 - 1 + k*ip)*ido) = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido];
+ }
+ }
+ }
+ }
+ } /* radfg */
+
+
+static void radbg(int ido, int ip, int l1, int idl1,
+ Treal cc[], Treal ch[], const Treal wa[])
+ {
+ static const Treal twopi = 6.28318530717959;
+ int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is;
+ Treal dc2, ai1, ai2, ar1, ar2, ds2;
+ int nbd;
+ Treal dcp, arg, dsp, ar1h, ar2h;
+ arg = twopi / ip;
+ dcp = cos(arg);
+ dsp = sin(arg);
+ nbd = (ido - 1) / 2;
+ ipph = (ip + 1) / 2;
+ if (ido >= l1) {
+ for (k=0; k<l1; k++) {
+ for (i=0; i<ido; i++) {
+ ch[i + k*ido] = ref(cc,i + k*ip*ido);
+ }
+ }
+ } else {
+ for (i=0; i<ido; i++) {
+ for (k=0; k<l1; k++) {
+ ch[i + k*ido] = ref(cc,i + k*ip*ido);
+ }
+ }
+ }
+ for (j=1; j<ipph; j++) {
+ jc = ip - j;
+ j2 = 2*j;
+ for (k=0; k<l1; k++) {
+ ch[(k + j*l1)*ido] = ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) + ref(cc,ido-1 + (j2 - 1 + k*ip)*
+ ido);
+ ch[(k + jc*l1)*ido] = ref(cc,(j2 + k*ip)*ido) + ref(cc,(j2 + k*ip)*ido);
+ }
+ }
+
+ if (ido != 1) {
+ if (nbd >= l1) {
+ for (j=1; j<ipph; j++) {
+ jc = ip - j;
+ for (k=0; k<l1; k++) {
+ for (i=2; i<ido; i+=2) {
+ ic = ido - i;
+ ch[i - 1 + (k + j*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc,
+ ic - 1 + (2*j - 1 + k*ip)*ido);
+ ch[i - 1 + (k + jc*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) -
+ ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);
+ ch[i + (k + j*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic
+ + (2*j - 1 + k*ip)*ido);
+ ch[i + (k + jc*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic
+ + (2*j - 1 + k*ip)*ido);
+ }
+ }
+ }
+ } else {
+ for (j=1; j<ipph; j++) {
+ jc = ip - j;
+ for (i=2; i<ido; i+=2) {
+ ic = ido - i;
+ for (k=0; k<l1; k++) {
+ ch[i - 1 + (k + j*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc,
+ ic - 1 + (2*j - 1 + k*ip)*ido);
+ ch[i - 1 + (k + jc*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) -
+ ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);
+ ch[i + (k + j*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic
+ + (2*j - 1 + k*ip)*ido);
+ ch[i + (k + jc*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic
+ + (2*j - 1 + k*ip)*ido);
+ }
+ }
+ }
+ }
+ }
+
+ ar1 = 1;
+ ai1 = 0;
+ for (l=1; l<ipph; l++) {
+ lc = ip - l;
+ ar1h = dcp*ar1 - dsp*ai1;
+ ai1 = dcp*ai1 + dsp*ar1;
+ ar1 = ar1h;
+ for (ik=0; ik<idl1; ik++) {
+ cc[ik + l*idl1] = ch[ik] + ar1*ch[ik + idl1];
+ cc[ik + lc*idl1] = ai1*ch[ik + (ip-1)*idl1];
+ }
+ dc2 = ar1;
+ ds2 = ai1;
+ ar2 = ar1;
+ ai2 = ai1;
+ for (j=2; j<ipph; j++) {
+ jc = ip - j;
+ ar2h = dc2*ar2 - ds2*ai2;
+ ai2 = dc2*ai2 + ds2*ar2;
+ ar2 = ar2h;
+ for (ik=0; ik<idl1; ik++) {
+ cc[ik + l*idl1] += ar2*ch[ik + j*idl1];
+ cc[ik + lc*idl1] += ai2*ch[ik + jc*idl1];
+ }
+ }
+ }
+ for (j=1; j<ipph; j++) {
+ for (ik=0; ik<idl1; ik++) {
+ ch[ik] += ch[ik + j*idl1];
+ }
+ }
+ for (j=1; j<ipph; j++) {
+ jc = ip - j;
+ for (k=0; k<l1; k++) {
+ ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido] - cc[(k + jc*l1)*ido];
+ ch[(k + jc*l1)*ido] = cc[(k + j*l1)*ido] + cc[(k + jc*l1)*ido];
+ }
+ }
+
+ if (ido == 1) return;
+ if (nbd >= l1) {
+ for (j=1; j<ipph; j++) {
+ jc = ip - j;
+ for (k=0; k<l1; k++) {
+ for (i=2; i<ido; i+=2) {
+ ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido];
+ ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] + cc[i + (k + jc*l1)*ido];
+ ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido];
+ ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido];
+ }
+ }
+ }
+ } else {
+ for (j=1; j<ipph; j++) {
+ jc = ip - j;
+ for (i=2; i<ido; i+=2) {
+ for (k=0; k<l1; k++) {
+ ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido];
+ ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j *l1)*ido] + cc[i + (k + jc*l1)*ido];
+ ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido];
+ ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido];
+ }
+ }
+ }
+ }
+ for (ik=0; ik<idl1; ik++) cc[ik] = ch[ik];
+ for (j=1; j<ip; j++)
+ for (k=0; k<l1; k++)
+ cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido];
+ if (nbd <= l1) {
+ is = -ido;
+ for (j=1; j<ip; j++) {
+ is += ido;
+ idij = is-1;
+ for (i=2; i<ido; i+=2) {
+ idij += 2;
+ for (k=0; k<l1; k++) {
+ cc[i - 1 + (k + j*l1)*ido] = wa[idij - 1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]*
+ ch[i + (k + j*l1)*ido];
+ cc[i + (k + j*l1)*ido] = wa[idij - 1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido];
+ }
+ }
+ }
+ } else {
+ is = -ido;
+ for (j=1; j<ip; j++) {
+ is += ido;
+ for (k=0; k<l1; k++) {
+ idij = is - 1;
+ for (i=2; i<ido; i+=2) {
+ idij += 2;
+ cc[i - 1 + (k + j*l1)*ido] = wa[idij-1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]*
+ ch[i + (k + j*l1)*ido];
+ cc[i + (k + j*l1)*ido] = wa[idij-1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido];
+ }
+ }
+ }
+ }
+ } /* radbg */
+
+ /* ----------------------------------------------------------------------
+cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs.
+---------------------------------------------------------------------- */
+
+static void cfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2], int isign)
+ {
+ int idot, i;
+ int k1, l1, l2;
+ int na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1;
+ Treal *cinput, *coutput;
+ nf = ifac[1];
+ na = 0;
+ l1 = 1;
+ iw = 0;
+ for (k1=2; k1<=nf+1; k1++) {
+ ip = ifac[k1];
+ l2 = ip*l1;
+ ido = n / l2;
+ idot = ido + ido;
+ idl1 = idot*l1;
+ if (na) {
+ cinput = ch;
+ coutput = c;
+ } else {
+ cinput = c;
+ coutput = ch;
+ }
+ switch (ip) {
+ case 4:
+ ix2 = iw + idot;
+ ix3 = ix2 + idot;
+ passf4(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], isign);
+ na = !na;
+ break;
+ case 2:
+ passf2(idot, l1, cinput, coutput, &wa[iw], isign);
+ na = !na;
+ break;
+ case 3:
+ ix2 = iw + idot;
+ passf3(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], isign);
+ na = !na;
+ break;
+ case 5:
+ ix2 = iw + idot;
+ ix3 = ix2 + idot;
+ ix4 = ix3 + idot;
+ passf5(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
+ na = !na;
+ break;
+ default:
+ passf(&nac, idot, ip, l1, idl1, cinput, coutput, &wa[iw], isign);
+ if (nac != 0) na = !na;
+ }
+ l1 = l2;
+ iw += (ip - 1)*idot;
+ }
+ if (na == 0) return;
+ for (i=0; i<2*n; i++) c[i] = ch[i];
+ } /* cfftf1 */
+
+
+void cfftf(int n, Treal c[], Treal wsave[])
+ {
+ int iw1, iw2;
+ if (n == 1) return;
+ iw1 = 2*n;
+ iw2 = iw1 + 2*n;
+ cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), -1);
+ } /* cfftf */
+
+
+void cfftb(int n, Treal c[], Treal wsave[])
+ {
+ int iw1, iw2;
+ if (n == 1) return;
+ iw1 = 2*n;
+ iw2 = iw1 + 2*n;
+ cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), +1);
+ } /* cfftb */
+
+
+static void factorize(int n, int ifac[MAXFAC+2], const int ntryh[NSPECIAL])
+ /* Factorize n in factors in ntryh and rest. On exit,
+ifac[0] contains n and ifac[1] contains number of factors,
+the factors start from ifac[2]. */
+ {
+ int ntry=3, i, j=0, ib, nf=0, nl=n, nq, nr;
+startloop:
+ if (j < NSPECIAL)
+ ntry = ntryh[j];
+ else
+ ntry+= 2;
+ j++;
+ do {
+ nq = nl / ntry;
+ nr = nl - ntry*nq;
+ if (nr != 0) goto startloop;
+ nf++;
+ ifac[nf + 1] = ntry;
+ nl = nq;
+ if (ntry == 2 && nf != 1) {
+ for (i=2; i<=nf; i++) {
+ ib = nf - i + 2;
+ ifac[ib + 1] = ifac[ib];
+ }
+ ifac[2] = 2;
+ }
+ } while (nl != 1);
+ ifac[0] = n;
+ ifac[1] = nf;
+ }
+
+
+static void cffti1(int n, Treal wa[], int ifac[MAXFAC+2])
+ {
+ static const Treal twopi = 6.28318530717959;
+ Treal arg, argh, argld, fi;
+ int idot, i, j;
+ int i1, k1, l1, l2;
+ int ld, ii, nf, ip;
+ int ido, ipm;
+
+ static const int ntryh[NSPECIAL] = {
+ 3,4,2,5 }; /* Do not change the order of these. */
+
+ factorize(n,ifac,ntryh);
+ nf = ifac[1];
+ argh = twopi/(Treal)n;
+ i = 1;
+ l1 = 1;
+ for (k1=1; k1<=nf; k1++) {
+ ip = ifac[k1+1];
+ ld = 0;
+ l2 = l1*ip;
+ ido = n / l2;
+ idot = ido + ido + 2;
+ ipm = ip - 1;
+ for (j=1; j<=ipm; j++) {
+ i1 = i;
+ wa[i-1] = 1;
+ wa[i] = 0;
+ ld += l1;
+ fi = 0;
+ argld = ld*argh;
+ for (ii=4; ii<=idot; ii+=2) {
+ i+= 2;
+ fi+= 1;
+ arg = fi*argld;
+ wa[i-1] = cos(arg);
+ wa[i] = sin(arg);
+ }
+ if (ip > 5) {
+ wa[i1-1] = wa[i-1];
+ wa[i1] = wa[i];
+ }
+ }
+ l1 = l2;
+ }
+ } /* cffti1 */
+
+
+void cffti(int n, Treal wsave[])
+ {
+ int iw1, iw2;
+ if (n == 1) return;
+ iw1 = 2*n;
+ iw2 = iw1 + 2*n;
+ cffti1(n, wsave+iw1, (int*)(wsave+iw2));
+ } /* cffti */
+
+ /* ----------------------------------------------------------------------
+rfftf1, rfftb1, rfftf, rfftb, rffti1, rffti. Treal FFTs.
+---------------------------------------------------------------------- */
+
+static void rfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2])
+ {
+ int i;
+ int k1, l1, l2, na, kh, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
+ Treal *cinput, *coutput;
+ nf = ifac[1];
+ na = 1;
+ l2 = n;
+ iw = n-1;
+ for (k1 = 1; k1 <= nf; ++k1) {
+ kh = nf - k1;
+ ip = ifac[kh + 2];
+ l1 = l2 / ip;
+ ido = n / l2;
+ idl1 = ido*l1;
+ iw -= (ip - 1)*ido;
+ na = !na;
+ if (na) {
+ cinput = ch;
+ coutput = c;
+ } else {
+ cinput = c;
+ coutput = ch;
+ }
+ switch (ip) {
+ case 4:
+ ix2 = iw + ido;
+ ix3 = ix2 + ido;
+ radf4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]);
+ break;
+ case 2:
+ radf2(ido, l1, cinput, coutput, &wa[iw]);
+ break;
+ case 3:
+ ix2 = iw + ido;
+ radf3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]);
+ break;
+ case 5:
+ ix2 = iw + ido;
+ ix3 = ix2 + ido;
+ ix4 = ix3 + ido;
+ radf5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
+ break;
+ default:
+ if (ido == 1)
+ na = !na;
+ if (na == 0) {
+ radfg(ido, ip, l1, idl1, c, ch, &wa[iw]);
+ na = 1;
+ } else {
+ radfg(ido, ip, l1, idl1, ch, c, &wa[iw]);
+ na = 0;
+ }
+ }
+ l2 = l1;
+ }
+ if (na == 1) return;
+ for (i = 0; i < n; i++) c[i] = ch[i];
+ } /* rfftf1 */
+
+
+void rfftb1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2])
+ {
+ int i;
+ int k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
+ Treal *cinput, *coutput;
+ nf = ifac[1];
+ na = 0;
+ l1 = 1;
+ iw = 0;
+ for (k1=1; k1<=nf; k1++) {
+ ip = ifac[k1 + 1];
+ l2 = ip*l1;
+ ido = n / l2;
+ idl1 = ido*l1;
+ if (na) {
+ cinput = ch;
+ coutput = c;
+ } else {
+ cinput = c;
+ coutput = ch;
+ }
+ switch (ip) {
+ case 4:
+ ix2 = iw + ido;
+ ix3 = ix2 + ido;
+ radb4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]);
+ na = !na;
+ break;
+ case 2:
+ radb2(ido, l1, cinput, coutput, &wa[iw]);
+ na = !na;
+ break;
+ case 3:
+ ix2 = iw + ido;
+ radb3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]);
+ na = !na;
+ break;
+ case 5:
+ ix2 = iw + ido;
+ ix3 = ix2 + ido;
+ ix4 = ix3 + ido;
+ radb5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
+ na = !na;
+ break;
+ default:
+ radbg(ido, ip, l1, idl1, cinput, coutput, &wa[iw]);
+ if (ido == 1) na = !na;
+ }
+ l1 = l2;
+ iw += (ip - 1)*ido;
+ }
+ if (na == 0) return;
+ for (i=0; i<n; i++) c[i] = ch[i];
+ } /* rfftb1 */
+
+
+void rfftf(int n, Treal r[], Treal wsave[])
+ {
+ if (n == 1) return;
+ rfftf1(n, r, wsave, wsave+n, (int*)(wsave+2*n));
+ } /* rfftf */
+
+
+void rfftb(int n, Treal r[], Treal wsave[])
+ {
+ if (n == 1) return;
+ rfftb1(n, r, wsave, wsave+n, (int*)(wsave+2*n));
+ } /* rfftb */
+
+
+static void rffti1(int n, Treal wa[], int ifac[MAXFAC+2])
+ {
+ static const Treal twopi = 6.28318530717959;
+ Treal arg, argh, argld, fi;
+ int i, j;
+ int k1, l1, l2;
+ int ld, ii, nf, ip, is;
+ int ido, ipm, nfm1;
+ static const int ntryh[NSPECIAL] = {
+ 4,2,3,5 }; /* Do not change the order of these. */
+ factorize(n,ifac,ntryh);
+ nf = ifac[1];
+ argh = twopi / n;
+ is = 0;
+ nfm1 = nf - 1;
+ l1 = 1;
+ if (nfm1 == 0) return;
+ for (k1 = 1; k1 <= nfm1; k1++) {
+ ip = ifac[k1 + 1];
+ ld = 0;
+ l2 = l1*ip;
+ ido = n / l2;
+ ipm = ip - 1;
+ for (j = 1; j <= ipm; ++j) {
+ ld += l1;
+ i = is;
+ argld = (Treal) ld*argh;
+ fi = 0;
+ for (ii = 3; ii <= ido; ii += 2) {
+ i += 2;
+ fi += 1;
+ arg = fi*argld;
+ wa[i - 2] = cos(arg);
+ wa[i - 1] = sin(arg);
+ }
+ is += ido;
+ }
+ l1 = l2;
+ }
+ } /* rffti1 */
+
+
+void rffti(int n, Treal wsave[])
+ {
+ if (n == 1) return;
+ rffti1(n, wsave+n, (int*)(wsave+2*n));
+ } /* rffti */
+
+#ifdef __cplusplus
+}
+#endif
diff --git a/numpy/fft/fftpack.h b/numpy/fft/fftpack.h
new file mode 100644
index 000000000..d134784a2
--- /dev/null
+++ b/numpy/fft/fftpack.h
@@ -0,0 +1,28 @@
+/*
+ * This file is part of tela the Tensor Language.
+ * Copyright (c) 1994-1995 Pekka Janhunen
+ */
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#define DOUBLE
+
+#ifdef DOUBLE
+#define Treal double
+#else
+#define Treal float
+#endif
+
+extern void cfftf(int N, Treal data[], const Treal wrk[]);
+extern void cfftb(int N, Treal data[], const Treal wrk[]);
+extern void cffti(int N, Treal wrk[]);
+
+extern void rfftf(int N, Treal data[], const Treal wrk[]);
+extern void rfftb(int N, Treal data[], const Treal wrk[]);
+extern void rffti(int N, Treal wrk[]);
+
+#ifdef __cplusplus
+}
+#endif
diff --git a/numpy/fft/fftpack.py b/numpy/fft/fftpack.py
new file mode 100644
index 000000000..0af1cdc79
--- /dev/null
+++ b/numpy/fft/fftpack.py
@@ -0,0 +1,550 @@
+"""
+Discrete Fourier Transforms - FFT.py
+
+The underlying code for these functions is an f2c translated and modified
+version of the FFTPACK routines.
+
+fft(a, n=None, axis=-1)
+ifft(a, n=None, axis=-1)
+rfft(a, n=None, axis=-1)
+irfft(a, n=None, axis=-1)
+hfft(a, n=None, axis=-1)
+ihfft(a, n=None, axis=-1)
+fftn(a, s=None, axes=None)
+ifftn(a, s=None, axes=None)
+rfftn(a, s=None, axes=None)
+irfftn(a, s=None, axes=None)
+fft2(a, s=None, axes=(-2,-1))
+ifft2(a, s=None, axes=(-2, -1))
+rfft2(a, s=None, axes=(-2,-1))
+irfft2(a, s=None, axes=(-2, -1))
+"""
+__all__ = ['fft','ifft', 'rfft', 'irfft', 'hfft', 'ihfft', 'rfftn',
+ 'irfftn', 'rfft2', 'irfft2', 'fft2', 'ifft2', 'fftn', 'ifftn',
+ 'refft', 'irefft','refftn','irefftn', 'refft2', 'irefft2']
+
+from numpy.core import asarray, zeros, swapaxes, shape, conjugate, \
+ take
+import fftpack_lite as fftpack
+from helper import *
+
+_fft_cache = {}
+_real_fft_cache = {}
+
+def _raw_fft(a, n=None, axis=-1, init_function=fftpack.cffti,
+ work_function=fftpack.cfftf, fft_cache = _fft_cache ):
+ a = asarray(a)
+
+ if n is None:
+ n = a.shape[axis]
+
+ if n < 1:
+ raise ValueError("Invalid number of FFT data points (%d) specified." % n)
+
+ try:
+ wsave = fft_cache[n]
+ except(KeyError):
+ wsave = init_function(n)
+ fft_cache[n] = wsave
+
+ if a.shape[axis] != n:
+ s = list(a.shape)
+ if s[axis] > n:
+ index = [slice(None)]*len(s)
+ index[axis] = slice(0,n)
+ a = a[index]
+ else:
+ index = [slice(None)]*len(s)
+ index[axis] = slice(0,s[axis])
+ s[axis] = n
+ z = zeros(s, a.dtype.char)
+ z[index] = a
+ a = z
+
+ if axis != -1:
+ a = swapaxes(a, axis, -1)
+ r = work_function(a, wsave)
+ if axis != -1:
+ r = swapaxes(r, axis, -1)
+ return r
+
+
+def fft(a, n=None, axis=-1):
+ """
+ Compute the one dimensional fft on a given axis.
+
+ Return the n point discrete Fourier transform of a. n defaults to the
+ length of a. If n is larger than the length of a, then a will be
+ zero-padded to make up the difference. If n is smaller than the length of
+ a, only the first n items in a will be used.
+
+ Parameters
+ ----------
+
+ a : array
+ input array
+ n : int
+ length of the fft
+ axis : int
+ axis over which to compute the fft
+
+ Notes
+ -----
+
+ The packing of the result is "standard": If A = fft(a, n), then A[0]
+ contains the zero-frequency term, A[1:n/2+1] contains the
+ positive-frequency terms, and A[n/2+1:] contains the negative-frequency
+ terms, in order of decreasingly negative frequency. So for an 8-point
+ transform, the frequencies of the result are [ 0, 1, 2, 3, 4, -3, -2, -1].
+
+ This is most efficient for n a power of two. This also stores a cache of
+ working memory for different sizes of fft's, so you could theoretically
+ run into memory problems if you call this too many times with too many
+ different n's.
+
+ """
+
+ return _raw_fft(a, n, axis, fftpack.cffti, fftpack.cfftf, _fft_cache)
+
+
+def ifft(a, n=None, axis=-1):
+ """
+ Compute the one-dimensonal inverse fft along an axis.
+
+ Return the `n` point inverse discrete Fourier transform of `a`. The length
+ `n` defaults to the length of `a`. If `n` is larger than the length of `a`,
+ then `a` will be zero-padded to make up the difference. If `n` is smaller
+ than the length of `a`, then `a` will be truncated to reduce its size.
+
+ Parameters
+ ----------
+
+ a : array_like
+ Input array.
+ n : int, optional
+ Length of the fft.
+ axis : int, optional
+ Axis over which to compute the inverse fft.
+
+ See Also
+ --------
+ fft
+
+ Notes
+ -----
+ The input array is expected to be packed the same way as the output of
+ fft, as discussed in the fft documentation.
+
+ This is the inverse of fft: ifft(fft(a)) == a within numerical
+ accuracy.
+
+ This is most efficient for `n` a power of two. This also stores a cache of
+ working memory for different sizes of fft's, so you could theoretically
+ run into memory problems if you call this too many times with too many
+ different `n` values.
+
+ """
+
+ a = asarray(a).astype(complex)
+ if n is None:
+ n = shape(a)[axis]
+ return _raw_fft(a, n, axis, fftpack.cffti, fftpack.cfftb, _fft_cache) / n
+
+
+def rfft(a, n=None, axis=-1):
+ """
+ Compute the one-dimensional fft for real input.
+
+ Return the n point discrete Fourier transform of the real valued
+ array a. n defaults to the length of a. n is the length of the
+ input, not the output.
+
+ Parameters
+ ----------
+
+ a : array
+ input array with real data type
+ n : int
+ length of the fft
+ axis : int
+ axis over which to compute the fft
+
+ Notes
+ -----
+
+ The returned array will be the nonnegative frequency terms of the
+ Hermite-symmetric, complex transform of the real array. So for an 8-point
+ transform, the frequencies in the result are [ 0, 1, 2, 3, 4]. The first
+ term will be real, as will the last if n is even. The negative frequency
+ terms are not needed because they are the complex conjugates of the
+ positive frequency terms. (This is what I mean when I say
+ Hermite-symmetric.)
+
+ This is most efficient for n a power of two.
+
+ """
+
+ a = asarray(a).astype(float)
+ return _raw_fft(a, n, axis, fftpack.rffti, fftpack.rfftf, _real_fft_cache)
+
+
+def irfft(a, n=None, axis=-1):
+ """
+ Compute the one-dimensional inverse fft for real input.
+
+ Parameters
+ ----------
+ a : array
+ Input array with real data type.
+ n : int
+ Length of the fft.
+ axis : int
+ Axis over which to compute the fft.
+
+ See Also
+ --------
+ rfft
+
+ Notes
+ -----
+ Return the real valued `n` point inverse discrete Fourier transform
+ of `a`, where `a` contains the nonnegative frequency terms of a
+ Hermite-symmetric sequence. `n` is the length of the result, not the
+ input. If `n` is not supplied, the default is 2*(len(`a`)-1). If you
+ want the length of the result to be odd, you have to say so.
+
+ If you specify an `n` such that `a` must be zero-padded or truncated, the
+ extra/removed values will be added/removed at high frequencies. One can
+ thus resample a series to m points via Fourier interpolation by: a_resamp
+ = irfft(rfft(a), m).
+
+ Within numerical accuracy ``irfft`` is the inverse of ``rfft``::
+
+ irfft(rfft(a), len(a)) == a
+
+ """
+
+ a = asarray(a).astype(complex)
+ if n is None:
+ n = (shape(a)[axis] - 1) * 2
+ return _raw_fft(a, n, axis, fftpack.rffti, fftpack.rfftb,
+ _real_fft_cache) / n
+
+
+def hfft(a, n=None, axis=-1):
+ """
+ Compute the fft of a signal which spectrum has Hermitian symmetry.
+
+ Parameters
+ ----------
+ a : array
+ input array
+ n : int
+ length of the hfft
+ axis : int
+ axis over which to compute the hfft
+
+ See also
+ --------
+ rfft
+ ihfft
+
+ Notes
+ -----
+ These are a pair analogous to rfft/irfft, but for the
+ opposite case: here the signal is real in the frequency domain and has
+ Hermite symmetry in the time domain. So here it's hermite_fft for which
+ you must supply the length of the result if it is to be odd.
+
+ ihfft(hfft(a), len(a)) == a
+ within numerical accuracy.
+
+ """
+
+ a = asarray(a).astype(complex)
+ if n is None:
+ n = (shape(a)[axis] - 1) * 2
+ return irfft(conjugate(a), n, axis) * n
+
+
+def ihfft(a, n=None, axis=-1):
+ """
+ Compute the inverse fft of a signal whose spectrum has Hermitian symmetry.
+
+ Parameters
+ ----------
+ a : array_like
+ Input array.
+ n : int, optional
+ Length of the ihfft.
+ axis : int, optional
+ Axis over which to compute the ihfft.
+
+ See also
+ --------
+ rfft, hfft
+
+ Notes
+ -----
+ These are a pair analogous to rfft/irfft, but for the
+ opposite case: here the signal is real in the frequency domain and has
+ Hermite symmetry in the time domain. So here it's hermite_fft for which
+ you must supply the length of the result if it is to be odd.
+
+ ihfft(hfft(a), len(a)) == a
+ within numerical accuracy.
+
+ """
+
+ a = asarray(a).astype(float)
+ if n is None:
+ n = shape(a)[axis]
+ return conjugate(rfft(a, n, axis))/n
+
+
+def _cook_nd_args(a, s=None, axes=None, invreal=0):
+ if s is None:
+ shapeless = 1
+ if axes is None:
+ s = list(a.shape)
+ else:
+ s = take(a.shape, axes)
+ else:
+ shapeless = 0
+ s = list(s)
+ if axes is None:
+ axes = range(-len(s), 0)
+ if len(s) != len(axes):
+ raise ValueError, "Shape and axes have different lengths."
+ if invreal and shapeless:
+ s[axes[-1]] = (s[axes[-1]] - 1) * 2
+ return s, axes
+
+
+def _raw_fftnd(a, s=None, axes=None, function=fft):
+ a = asarray(a)
+ s, axes = _cook_nd_args(a, s, axes)
+ itl = range(len(axes))
+ itl.reverse()
+ for ii in itl:
+ a = function(a, n=s[ii], axis=axes[ii])
+ return a
+
+
+def fftn(a, s=None, axes=None):
+ """
+ Compute the N-dimensional Fast Fourier Transform.
+
+ Parameters
+ ----------
+ a : array_like
+ Input array.
+ s : sequence of ints
+ Shape of each axis of the input (s[0] refers to axis 0, s[1] to
+ axis 1, etc.). This corresponds to `n` for `fft(x, n)`.
+ Along any axis, if the given shape is smaller than that of the input,
+ the input is cropped. If it is larger, the input is padded with zeros.
+ axes : tuple of int
+ Axes over which to compute the FFT.
+
+ Notes
+ -----
+ Analogously to `fft`, the term for zero frequency in all axes is in the
+ low-order corner, while the term for the Nyquist frequency in all axes is
+ in the middle.
+
+ If neither `s` nor `axes` is specified, the transform is taken along all
+ axes. If `s` is specified and `axes` is not, the last ``len(s)`` axes are
+ used. If `axes` is specified and `s` is not, the input shape along the
+ specified axes is used. If `s` and `axes` are both specified and are not
+ the same length, an exception is raised.
+
+ """
+
+ return _raw_fftnd(a,s,axes,fft)
+
+def ifftn(a, s=None, axes=None):
+ """
+ Compute the inverse of fftn.
+
+ Parameters
+ ----------
+ a : array
+ input array
+ s : sequence (int)
+ shape of the ifft
+ axis : int
+ axis over which to compute the ifft
+
+ Notes
+ -----
+ The n-dimensional ifft of a. s is a sequence giving the shape of the input
+ an result along the transformed axes, as n for fft. Results are packed
+ analogously to fft: the term for zero frequency in all axes is in the
+ low-order corner, while the term for the Nyquist frequency in all axes is
+ in the middle.
+
+ If neither s nor axes is specified, the transform is taken along all
+ axes. If s is specified and axes is not, the last len(s) axes are used.
+ If axes are specified and s is not, the input shape along the specified
+ axes is used. If s and axes are both specified and are not the same
+ length, an exception is raised.
+
+ """
+
+ return _raw_fftnd(a, s, axes, ifft)
+
+
+def fft2(a, s=None, axes=(-2,-1)):
+ """
+ Compute the 2-D FFT of an array.
+
+ Parameters
+ ----------
+ a : array_like
+ Input array. The rank (dimensions) of `a` must be 2 or greater.
+ s : shape tuple
+ Shape of the FFT.
+ axes : sequence of 2 ints
+ The 2 axes over which to compute the FFT. The default is the last two
+ axes (-2, -1).
+
+ Notes
+ -----
+ This is really just ``fftn`` with different default behavior.
+
+ """
+
+ return _raw_fftnd(a,s,axes,fft)
+
+
+def ifft2(a, s=None, axes=(-2,-1)):
+ """
+ Compute the inverse 2d fft of an array.
+
+ Parameters
+ ----------
+ a : array
+ input array
+ s : sequence (int)
+ shape of the ifft
+ axis : int
+ axis over which to compute the ifft
+
+ Notes
+ -----
+ This is really just ifftn with different default behavior.
+
+ """
+
+ return _raw_fftnd(a, s, axes, ifft)
+
+
+def rfftn(a, s=None, axes=None):
+ """
+ Compute the n-dimensional fft of a real array.
+
+ Parameters
+ ----------
+ a : array (real)
+ input array
+ s : sequence (int)
+ shape of the fft
+ axis : int
+ axis over which to compute the fft
+
+ Notes
+ -----
+ A real transform as rfft is performed along the axis specified by the last
+ element of axes, then complex transforms as fft are performed along the
+ other axes.
+
+ """
+
+ a = asarray(a).astype(float)
+ s, axes = _cook_nd_args(a, s, axes)
+ a = rfft(a, s[-1], axes[-1])
+ for ii in range(len(axes)-1):
+ a = fft(a, s[ii], axes[ii])
+ return a
+
+def rfft2(a, s=None, axes=(-2,-1)):
+ """
+ Compute the 2-dimensional fft of a real array.
+
+ Parameters
+ ----------
+ a : array (real)
+ input array
+ s : sequence (int)
+ shape of the fft
+ axis : int
+ axis over which to compute the fft
+
+ Notes
+ -----
+ The 2-D fft of the real valued array a. This is really just rfftn with
+ different default behavior.
+
+ """
+
+ return rfftn(a, s, axes)
+
+def irfftn(a, s=None, axes=None):
+ """
+ Compute the n-dimensional inverse fft of a real array.
+
+ Parameters
+ ----------
+ a : array (real)
+ input array
+ s : sequence (int)
+ shape of the inverse fft
+ axis : int
+ axis over which to compute the inverse fft
+
+ Notes
+ -----
+ The transform implemented in ifftn is applied along
+ all axes but the last, then the transform implemented in irfft is performed
+ along the last axis. As with irfft, the length of the result along that
+ axis must be specified if it is to be odd.
+
+ """
+
+ a = asarray(a).astype(complex)
+ s, axes = _cook_nd_args(a, s, axes, invreal=1)
+ for ii in range(len(axes)-1):
+ a = ifft(a, s[ii], axes[ii])
+ a = irfft(a, s[-1], axes[-1])
+ return a
+
+def irfft2(a, s=None, axes=(-2,-1)):
+ """
+ Compute the 2-dimensional inverse fft of a real array.
+
+ Parameters
+ ----------
+ a : array (real)
+ input array
+ s : sequence (int)
+ shape of the inverse fft
+ axis : int
+ axis over which to compute the inverse fft
+
+ Notes
+ -----
+ This is really irfftn with different default.
+
+ """
+
+ return irfftn(a, s, axes)
+
+# Deprecated names
+from numpy import deprecate
+refft = deprecate(rfft, 'refft', 'rfft')
+irefft = deprecate(irfft, 'irefft', 'irfft')
+refft2 = deprecate(rfft2, 'refft2', 'rfft2')
+irefft2 = deprecate(irfft2, 'irefft2', 'irfft2')
+refftn = deprecate(rfftn, 'refftn', 'rfftn')
+irefftn = deprecate(irfftn, 'irefftn', 'irfftn')
diff --git a/numpy/fft/fftpack_litemodule.c b/numpy/fft/fftpack_litemodule.c
new file mode 100644
index 000000000..3375976d6
--- /dev/null
+++ b/numpy/fft/fftpack_litemodule.c
@@ -0,0 +1,330 @@
+#include "fftpack.h"
+#include "Python.h"
+#include "numpy/arrayobject.h"
+
+static PyObject *ErrorObject;
+
+/* ----------------------------------------------------- */
+
+static char fftpack_cfftf__doc__[] = "";
+
+PyObject *
+fftpack_cfftf(PyObject *NPY_UNUSED(self), PyObject *args)
+{
+ PyObject *op1, *op2;
+ PyArrayObject *data;
+ PyArray_Descr *descr;
+ double *wsave, *dptr;
+ npy_intp nsave;
+ int npts, nrepeats, i;
+
+ if(!PyArg_ParseTuple(args, "OO", &op1, &op2)) {
+ return NULL;
+ }
+ data = (PyArrayObject *)PyArray_CopyFromObject(op1,
+ PyArray_CDOUBLE, 1, 0);
+ if (data == NULL) {
+ return NULL;
+ }
+ descr = PyArray_DescrFromType(PyArray_DOUBLE);
+ if (PyArray_AsCArray(&op2, (void *)&wsave, &nsave, 1, descr) == -1) {
+ goto fail;
+ }
+ if (data == NULL) {
+ goto fail;
+ }
+
+ npts = data->dimensions[data->nd - 1];
+ if (nsave != npts*4 + 15) {
+ PyErr_SetString(ErrorObject, "invalid work array for fft size");
+ goto fail;
+ }
+
+ nrepeats = PyArray_SIZE(data)/npts;
+ dptr = (double *)data->data;
+ NPY_SIGINT_ON;
+ for (i = 0; i < nrepeats; i++) {
+ cfftf(npts, dptr, wsave);
+ dptr += npts*2;
+ }
+ NPY_SIGINT_OFF;
+ PyArray_Free(op2, (char *)wsave);
+ return (PyObject *)data;
+
+fail:
+ PyArray_Free(op2, (char *)wsave);
+ Py_DECREF(data);
+ return NULL;
+}
+
+static char fftpack_cfftb__doc__[] = "";
+
+PyObject *
+fftpack_cfftb(PyObject *NPY_UNUSED(self), PyObject *args)
+{
+ PyObject *op1, *op2;
+ PyArrayObject *data;
+ PyArray_Descr *descr;
+ double *wsave, *dptr;
+ npy_intp nsave;
+ int npts, nrepeats, i;
+
+ if(!PyArg_ParseTuple(args, "OO", &op1, &op2)) {
+ return NULL;
+ }
+ data = (PyArrayObject *)PyArray_CopyFromObject(op1,
+ PyArray_CDOUBLE, 1, 0);
+ if (data == NULL) {
+ return NULL;
+ }
+ descr = PyArray_DescrFromType(PyArray_DOUBLE);
+ if (PyArray_AsCArray(&op2, (void *)&wsave, &nsave, 1, descr) == -1) {
+ goto fail;
+ }
+ if (data == NULL) {
+ goto fail;
+ }
+
+ npts = data->dimensions[data->nd - 1];
+ if (nsave != npts*4 + 15) {
+ PyErr_SetString(ErrorObject, "invalid work array for fft size");
+ goto fail;
+ }
+
+ nrepeats = PyArray_SIZE(data)/npts;
+ dptr = (double *)data->data;
+ NPY_SIGINT_ON;
+ for (i = 0; i < nrepeats; i++) {
+ cfftb(npts, dptr, wsave);
+ dptr += npts*2;
+ }
+ NPY_SIGINT_OFF;
+ PyArray_Free(op2, (char *)wsave);
+ return (PyObject *)data;
+
+fail:
+ PyArray_Free(op2, (char *)wsave);
+ Py_DECREF(data);
+ return NULL;
+}
+
+static char fftpack_cffti__doc__[] ="";
+
+static PyObject *
+fftpack_cffti(PyObject *NPY_UNUSED(self), PyObject *args)
+{
+ PyArrayObject *op;
+ npy_intp dim;
+ long n;
+
+ if (!PyArg_ParseTuple(args, "l", &n)) {
+ return NULL;
+ }
+ /*Magic size needed by cffti*/
+ dim = 4*n + 15;
+ /*Create a 1 dimensional array of dimensions of type double*/
+ op = (PyArrayObject *)PyArray_SimpleNew(1, &dim, PyArray_DOUBLE);
+ if (op == NULL) {
+ return NULL;
+ }
+
+ NPY_SIGINT_ON;
+ cffti(n, (double *)((PyArrayObject*)op)->data);
+ NPY_SIGINT_OFF;
+
+ return (PyObject *)op;
+}
+
+static char fftpack_rfftf__doc__[] ="";
+
+PyObject *
+fftpack_rfftf(PyObject *NPY_UNUSED(self), PyObject *args)
+{
+ PyObject *op1, *op2;
+ PyArrayObject *data, *ret;
+ PyArray_Descr *descr;
+ double *wsave, *dptr, *rptr;
+ npy_intp nsave;
+ int npts, nrepeats, i, rstep;
+
+ if(!PyArg_ParseTuple(args, "OO", &op1, &op2)) {
+ return NULL;
+ }
+ data = (PyArrayObject *)PyArray_ContiguousFromObject(op1,
+ PyArray_DOUBLE, 1, 0);
+ if (data == NULL) {
+ return NULL;
+ }
+ npts = data->dimensions[data->nd-1];
+ data->dimensions[data->nd - 1] = npts/2 + 1;
+ ret = (PyArrayObject *)PyArray_Zeros(data->nd, data->dimensions,
+ PyArray_DescrFromType(PyArray_CDOUBLE), 0);
+ data->dimensions[data->nd - 1] = npts;
+ rstep = (ret->dimensions[ret->nd - 1])*2;
+
+ descr = PyArray_DescrFromType(PyArray_DOUBLE);
+ if (PyArray_AsCArray(&op2, (void *)&wsave, &nsave, 1, descr) == -1) {
+ goto fail;
+ }
+ if (data == NULL || ret == NULL) {
+ goto fail;
+ }
+ if (nsave != npts*2+15) {
+ PyErr_SetString(ErrorObject, "invalid work array for fft size");
+ goto fail;
+ }
+
+ nrepeats = PyArray_SIZE(data)/npts;
+ rptr = (double *)ret->data;
+ dptr = (double *)data->data;
+
+
+ NPY_SIGINT_ON;
+ for (i = 0; i < nrepeats; i++) {
+ memcpy((char *)(rptr+1), dptr, npts*sizeof(double));
+ rfftf(npts, rptr+1, wsave);
+ rptr[0] = rptr[1];
+ rptr[1] = 0.0;
+ rptr += rstep;
+ dptr += npts;
+ }
+ NPY_SIGINT_OFF;
+ PyArray_Free(op2, (char *)wsave);
+ Py_DECREF(data);
+ return (PyObject *)ret;
+
+fail:
+ PyArray_Free(op2, (char *)wsave);
+ Py_XDECREF(data);
+ Py_XDECREF(ret);
+ return NULL;
+}
+
+static char fftpack_rfftb__doc__[] ="";
+
+
+PyObject *
+fftpack_rfftb(PyObject *NPY_UNUSED(self), PyObject *args)
+{
+ PyObject *op1, *op2;
+ PyArrayObject *data, *ret;
+ PyArray_Descr *descr;
+ double *wsave, *dptr, *rptr;
+ npy_intp nsave;
+ int npts, nrepeats, i;
+
+ if(!PyArg_ParseTuple(args, "OO", &op1, &op2)) {
+ return NULL;
+ }
+ data = (PyArrayObject *)PyArray_ContiguousFromObject(op1,
+ PyArray_CDOUBLE, 1, 0);
+ if (data == NULL) {
+ return NULL;
+ }
+ npts = data->dimensions[data->nd - 1];
+ ret = (PyArrayObject *)PyArray_Zeros(data->nd, data->dimensions,
+ PyArray_DescrFromType(PyArray_DOUBLE), 0);
+
+ descr = PyArray_DescrFromType(PyArray_DOUBLE);
+ if (PyArray_AsCArray(&op2, (void *)&wsave, &nsave, 1, descr) == -1) {
+ goto fail;
+ }
+ if (data == NULL || ret == NULL) {
+ goto fail;
+ }
+ if (nsave != npts*2 + 15) {
+ PyErr_SetString(ErrorObject, "invalid work array for fft size");
+ goto fail;
+ }
+
+ nrepeats = PyArray_SIZE(ret)/npts;
+ rptr = (double *)ret->data;
+ dptr = (double *)data->data;
+
+ NPY_SIGINT_ON;
+ for (i = 0; i < nrepeats; i++) {
+ memcpy((char *)(rptr + 1), (dptr + 2), (npts - 1)*sizeof(double));
+ rptr[0] = dptr[0];
+ rfftb(npts, rptr, wsave);
+ rptr += npts;
+ dptr += npts*2;
+ }
+ NPY_SIGINT_OFF;
+ PyArray_Free(op2, (char *)wsave);
+ Py_DECREF(data);
+ return (PyObject *)ret;
+
+fail:
+ PyArray_Free(op2, (char *)wsave);
+ Py_XDECREF(data);
+ Py_XDECREF(ret);
+ return NULL;
+}
+
+
+static char fftpack_rffti__doc__[] ="";
+
+static PyObject *
+fftpack_rffti(PyObject *NPY_UNUSED(self), PyObject *args)
+{
+ PyArrayObject *op;
+ npy_intp dim;
+ long n;
+
+ if (!PyArg_ParseTuple(args, "l", &n)) {
+ return NULL;
+ }
+ /*Magic size needed by rffti*/
+ dim = 2*n + 15;
+ /*Create a 1 dimensional array of dimensions of type double*/
+ op = (PyArrayObject *)PyArray_SimpleNew(1, &dim, PyArray_DOUBLE);
+ if (op == NULL) {
+ return NULL;
+ }
+ NPY_SIGINT_ON;
+ rffti(n, (double *)((PyArrayObject*)op)->data);
+ NPY_SIGINT_OFF;
+
+ return (PyObject *)op;
+}
+
+
+/* List of methods defined in the module */
+
+static struct PyMethodDef fftpack_methods[] = {
+ {"cfftf", fftpack_cfftf, 1, fftpack_cfftf__doc__},
+ {"cfftb", fftpack_cfftb, 1, fftpack_cfftb__doc__},
+ {"cffti", fftpack_cffti, 1, fftpack_cffti__doc__},
+ {"rfftf", fftpack_rfftf, 1, fftpack_rfftf__doc__},
+ {"rfftb", fftpack_rfftb, 1, fftpack_rfftb__doc__},
+ {"rffti", fftpack_rffti, 1, fftpack_rffti__doc__},
+ {NULL, NULL} /* sentinel */
+};
+
+
+/* Initialization function for the module (*must* be called initfftpack) */
+
+static char fftpack_module_documentation[] =
+""
+;
+
+PyMODINIT_FUNC initfftpack_lite(void)
+{
+ PyObject *m, *d;
+
+ /* Create the module and add the functions */
+ m = Py_InitModule4("fftpack_lite", fftpack_methods,
+ fftpack_module_documentation,
+ (PyObject*)NULL,PYTHON_API_VERSION);
+
+ /* Import the array object */
+ import_array();
+
+ /* Add some symbolic constants to the module */
+ d = PyModule_GetDict(m);
+ ErrorObject = PyErr_NewException("fftpack.error", NULL, NULL);
+ PyDict_SetItemString(d, "error", ErrorObject);
+
+ /* XXXX Add constants here */
+
+}
diff --git a/numpy/fft/helper.py b/numpy/fft/helper.py
new file mode 100644
index 000000000..e2e36c323
--- /dev/null
+++ b/numpy/fft/helper.py
@@ -0,0 +1,116 @@
+"""
+Discrete Fourier Transforms - helper.py
+"""
+# Created by Pearu Peterson, September 2002
+
+__all__ = ['fftshift','ifftshift','fftfreq']
+
+from numpy.core import asarray, concatenate, arange, take, \
+ integer, empty
+import types
+
+def fftshift(x,axes=None):
+ """
+ Shift zero-frequency component to center of spectrum.
+
+ This function swaps half-spaces for all axes listed (defaults to all).
+ If len(x) is even then the Nyquist component is y[0].
+
+ Parameters
+ ----------
+ x : array_like
+ Input array.
+ axes : int or shape tuple, optional
+ Axes over which to shift. Default is None which shifts all axes.
+
+ See Also
+ --------
+ ifftshift
+
+ """
+ tmp = asarray(x)
+ ndim = len(tmp.shape)
+ if axes is None:
+ axes = range(ndim)
+ y = tmp
+ for k in axes:
+ n = tmp.shape[k]
+ p2 = (n+1)/2
+ mylist = concatenate((arange(p2,n),arange(p2)))
+ y = take(y,mylist,k)
+ return y
+
+
+def ifftshift(x,axes=None):
+ """
+ Inverse of fftshift.
+
+ Parameters
+ ----------
+ x : array_like
+ Input array.
+ axes : int or shape tuple, optional
+ Axes over which to calculate. Defaults to None which is over all axes.
+
+ See Also
+ --------
+ fftshift
+
+ """
+ tmp = asarray(x)
+ ndim = len(tmp.shape)
+ if axes is None:
+ axes = range(ndim)
+ y = tmp
+ for k in axes:
+ n = tmp.shape[k]
+ p2 = n-(n+1)/2
+ mylist = concatenate((arange(p2,n),arange(p2)))
+ y = take(y,mylist,k)
+ return y
+
+def fftfreq(n,d=1.0):
+ """
+ Discrete Fourier Transform sample frequencies.
+
+ The returned float array contains the frequency bins in
+ cycles/unit (with zero at the start) given a window length `n` and a
+ sample spacing `d`.
+ ::
+
+ f = [0,1,...,n/2-1,-n/2,...,-1]/(d*n) if n is even
+ f = [0,1,...,(n-1)/2,-(n-1)/2,...,-1]/(d*n) if n is odd
+
+ Parameters
+ ----------
+ n : int
+ Window length.
+ d : scalar
+ Sample spacing.
+
+ Returns
+ -------
+ out : ndarray, shape(`n`,)
+ Sample frequencies.
+
+ Examples
+ --------
+ >>> signal = np.array([-2., 8., -6., 4., 1., 0., 3., 5.])
+ >>> fourier = np.fft.fft(signal)
+ >>> n = len(signal)
+ >>> timestep = 0.1
+ >>> freq = np.fft.fftfreq(n, d=timestep)
+ >>> freq
+ array([ 0. , 1.25, 2.5 , 3.75, -5. , -3.75, -2.5 , -1.25])
+
+ """
+ assert isinstance(n,types.IntType) or isinstance(n, integer)
+ val = 1.0/(n*d)
+ results = empty(n, int)
+ N = (n-1)//2 + 1
+ p1 = arange(0,N,dtype=int)
+ results[:N] = p1
+ p2 = arange(-(n//2),0,dtype=int)
+ results[N:] = p2
+ return results * val
+ #return hstack((arange(0,(n-1)/2 + 1), arange(-(n/2),0))) / (n*d)
diff --git a/numpy/fft/info.py b/numpy/fft/info.py
new file mode 100644
index 000000000..c9ab599a4
--- /dev/null
+++ b/numpy/fft/info.py
@@ -0,0 +1,29 @@
+"""\
+Core FFT routines
+==================
+
+ Standard FFTs
+
+ fft
+ ifft
+ fft2
+ ifft2
+ fftn
+ ifftn
+
+ Real FFTs
+
+ rfft
+ irfft
+ rfft2
+ irfft2
+ rfftn
+ irfftn
+
+ Hermite FFTs
+
+ hfft
+ ihfft
+"""
+
+depends = ['core']
diff --git a/numpy/fft/setup.py b/numpy/fft/setup.py
new file mode 100644
index 000000000..6acad7c9a
--- /dev/null
+++ b/numpy/fft/setup.py
@@ -0,0 +1,19 @@
+
+
+def configuration(parent_package='',top_path=None):
+ from numpy.distutils.misc_util import Configuration
+ config = Configuration('fft',parent_package,top_path)
+
+ config.add_data_dir('tests')
+
+ # Configure fftpack_lite
+ config.add_extension('fftpack_lite',
+ sources=['fftpack_litemodule.c', 'fftpack.c']
+ )
+
+
+ return config
+
+if __name__ == '__main__':
+ from numpy.distutils.core import setup
+ setup(configuration=configuration)
diff --git a/numpy/fft/setupscons.py b/numpy/fft/setupscons.py
new file mode 100644
index 000000000..54551b0a3
--- /dev/null
+++ b/numpy/fft/setupscons.py
@@ -0,0 +1,15 @@
+def configuration(parent_package = '', top_path = None):
+ from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs
+ config = Configuration('fft', parent_package, top_path)
+
+ config.add_data_dir('tests')
+
+ config.add_sconscript('SConstruct',
+ source_files = ['fftpack_litemodule.c', 'fftpack.c',
+ 'fftpack.h'])
+
+ return config
+
+if __name__ == '__main__':
+ from numpy.distutils.core import setup
+ setup(configuration=configuration)
diff --git a/numpy/fft/tests/test_fftpack.py b/numpy/fft/tests/test_fftpack.py
new file mode 100644
index 000000000..1095de4bc
--- /dev/null
+++ b/numpy/fft/tests/test_fftpack.py
@@ -0,0 +1,23 @@
+from numpy.testing import *
+import numpy as np
+
+def fft1(x):
+ L = len(x)
+ phase = -2j*np.pi*(np.arange(L)/float(L))
+ phase = np.arange(L).reshape(-1,1) * phase
+ return np.sum(x*np.exp(phase),axis=1)
+
+class TestFFTShift(TestCase):
+ def test_fft_n(self):
+ self.failUnlessRaises(ValueError,np.fft.fft,[1,2,3],0)
+
+
+class TestFFT1D(TestCase):
+ def test_basic(self):
+ rand = np.random.random
+ x = rand(30) + 1j*rand(30)
+ assert_array_almost_equal(fft1(x), np.fft.fft(x))
+
+
+if __name__ == "__main__":
+ run_module_suite()
diff --git a/numpy/fft/tests/test_helper.py b/numpy/fft/tests/test_helper.py
new file mode 100644
index 000000000..f757b6032
--- /dev/null
+++ b/numpy/fft/tests/test_helper.py
@@ -0,0 +1,42 @@
+#!/usr/bin/env python
+# Copied from fftpack.helper by Pearu Peterson, October 2005
+""" Test functions for fftpack.helper module
+"""
+
+from numpy.testing import *
+from numpy.fft import fftshift,ifftshift,fftfreq
+
+from numpy import pi
+
+def random(size):
+ return rand(*size)
+
+class TestFFTShift(TestCase):
+ def test_definition(self):
+ x = [0,1,2,3,4,-4,-3,-2,-1]
+ y = [-4,-3,-2,-1,0,1,2,3,4]
+ assert_array_almost_equal(fftshift(x),y)
+ assert_array_almost_equal(ifftshift(y),x)
+ x = [0,1,2,3,4,-5,-4,-3,-2,-1]
+ y = [-5,-4,-3,-2,-1,0,1,2,3,4]
+ assert_array_almost_equal(fftshift(x),y)
+ assert_array_almost_equal(ifftshift(y),x)
+
+ def test_inverse(self):
+ for n in [1,4,9,100,211]:
+ x = random((n,))
+ assert_array_almost_equal(ifftshift(fftshift(x)),x)
+
+
+class TestFFTFreq(TestCase):
+ def test_definition(self):
+ x = [0,1,2,3,4,-4,-3,-2,-1]
+ assert_array_almost_equal(9*fftfreq(9),x)
+ assert_array_almost_equal(9*pi*fftfreq(9,pi),x)
+ x = [0,1,2,3,4,-5,-4,-3,-2,-1]
+ assert_array_almost_equal(10*fftfreq(10),x)
+ assert_array_almost_equal(10*pi*fftfreq(10,pi),x)
+
+
+if __name__ == "__main__":
+ run_module_suite()
View Lorry Controller Status