diff options
Diffstat (limited to 'numpy/fft')
-rw-r--r-- | numpy/fft/__init__.py | 9 | ||||
-rw-r--r-- | numpy/fft/fftpack.c | 1501 | ||||
-rw-r--r-- | numpy/fft/fftpack.h | 28 | ||||
-rw-r--r-- | numpy/fft/fftpack.py | 323 | ||||
-rw-r--r-- | numpy/fft/fftpack_litemodule.c | 262 | ||||
-rw-r--r-- | numpy/fft/helper.py | 66 | ||||
-rw-r--r-- | numpy/fft/info.py | 29 | ||||
-rw-r--r-- | numpy/fft/old.py | 19 | ||||
-rw-r--r-- | numpy/fft/setup.py | 20 | ||||
-rw-r--r-- | numpy/fft/tests/test_helper.py | 45 |
10 files changed, 2302 insertions, 0 deletions
diff --git a/numpy/fft/__init__.py b/numpy/fft/__init__.py new file mode 100644 index 000000000..0c5a7f2ed --- /dev/null +++ b/numpy/fft/__init__.py @@ -0,0 +1,9 @@ +# To get sub-modules +from info import __doc__ + +from fftpack import * +from helper import * + +def test(level=1, verbosity=1): + from numpy.testing import NumpyTest + return NumpyTest().test(level, verbosity) diff --git a/numpy/fft/fftpack.c b/numpy/fft/fftpack.c new file mode 100644 index 000000000..3e5e7d2ed --- /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, nt, 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..1cb24f2b7 --- /dev/null +++ b/numpy/fft/fftpack.py @@ -0,0 +1,323 @@ +""" +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 == None: n = a.shape[axis] + + 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): + """fft(a, n=None, axis=-1) + + Will return the n point discrete Fourier transform of a. n defaults to the + length of a. If n is larger than a, then a will be zero-padded to make up + the difference. If n is smaller than a, the first n items in a will be + used. + + 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): + """ifft(a, n=None, axis=-1) + + Will return the n point inverse discrete Fourier transform of a. n + defaults to the length of a. If n is larger than a, then a will be + zero-padded to make up the difference. If n is smaller than a, then a will + be truncated to reduce its size. + + The input array is expected to be packed the same way as the output of + fft, as discussed in it's 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's.""" + + a = asarray(a).astype(complex) + if n == 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): + """rfft(a, n=None, axis=-1) + + Will 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. + + 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): + """irfft(a, n=None, axis=-1) + + Will 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). + + This is the inverse of rfft: + irfft(rfft(a), len(a)) == a + within numerical accuracy.""" + + a = asarray(a).astype(complex) + if n == 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): + """hfft(a, n=None, axis=-1) + ihfft(a, n=None, axis=-1) + + 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 == None: + n = (shape(a)[axis] - 1) * 2 + return irfft(conjugate(a), n, axis) * n + + +def ihfft(a, n=None, axis=-1): + """hfft(a, n=None, axis=-1) + ihfft(a, n=None, axis=-1) + + 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 hfft 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 == 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 == None: + s = list(a.shape) + else: + s = take(a.shape, axes) + else: + shapeless = 0 + s = list(s) + if axes == 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): + """fftn(a, s=None, axes=None) + + The n-dimensional fft 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,fft) + +def ifftn(a, s=None, axes=None): + """ifftn(a, s=None, axes=None) + + The inverse of fftn.""" + + return _raw_fftnd(a, s, axes, ifft) + + +def fft2(a, s=None, axes=(-2,-1)): + """fft2(a, s=None, axes=(-2,-1)) + + The 2d fft of a. This is really just fftn with different default + behavior.""" + + return _raw_fftnd(a,s,axes,fft) + + +def ifft2(a, s=None, axes=(-2,-1)): + """ifft2(a, s=None, axes=(-2, -1)) + + The inverse of fft2d. This is really just ifftn with different + default behavior.""" + + return _raw_fftnd(a, s, axes, ifft) + + +def rfftn(a, s=None, axes=None): + """rfftn(a, s=None, axes=None) + + The n-dimensional discrete Fourier transform of a real array a. 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)): + """rfft2(a, s=None, axes=(-2,-1)) + + The 2d 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): + """irfftn(a, s=None, axes=None) + + The inverse of rfftn. The transform implemented in ifft 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)): + """irfft2(a, s=None, axes=(-2, -1)) + + The inverse of rfft2. This is really just irfftn with + different default behavior.""" + + 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..f6d9aaf6b --- /dev/null +++ b/numpy/fft/fftpack_litemodule.c @@ -0,0 +1,262 @@ +#include "fftpack.h" +#include "Python.h" +#include "numpy/arrayobject.h" + +static PyObject *ErrorObject; + +/* ----------------------------------------------------- */ + +static char fftpack_cfftf__doc__[] =""; + +PyObject * +fftpack_cfftf(PyObject *self, PyObject *args) +{ + PyObject *op1, *op2; + PyArrayObject *data; + double *wsave, *dptr; + int npts, nsave, 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; + if (PyArray_As1D(&op2, (char **)&wsave, &nsave, PyArray_DOUBLE) == -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; + for (i=0; i<nrepeats; i++) { + cfftf(npts, dptr, wsave); + dptr += npts*2; + } + 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 *self, PyObject *args) +{ + PyObject *op1, *op2; + PyArrayObject *data; + double *wsave, *dptr; + int npts, nsave, 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; + if (PyArray_As1D(&op2, (char **)&wsave, &nsave, PyArray_DOUBLE) == -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; + for (i=0; i<nrepeats; i++) { + cfftb(npts, dptr, wsave); + dptr += npts*2; + } + 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 *self, PyObject *args) +{ + PyArrayObject *op; + int dim, n; + + if (!PyArg_ParseTuple(args, "i", &n)) return NULL; + + dim = 4*n+15; /*Magic size needed by cffti*/ + /*Create a 1 dimensional array of dimensions of type double*/ + op = (PyArrayObject *)PyArray_FromDims(1, &dim, PyArray_DOUBLE); + if (op == NULL) return NULL; + + cffti(n, (double *)((PyArrayObject*)op)->data); + + return (PyObject *)op; +} + +static char fftpack_rfftf__doc__[] =""; + +PyObject * +fftpack_rfftf(PyObject *self, PyObject *args) +{ + PyObject *op1, *op2; + PyArrayObject *data, *ret; + double *wsave, *dptr, *rptr; + int npts, nsave, 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; + + if (PyArray_As1D(&op2, (char **)&wsave, &nsave, PyArray_DOUBLE) == -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; + + 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; + } + 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 *self, PyObject *args) +{ + PyObject *op1, *op2; + PyArrayObject *data, *ret; + double *wsave, *dptr, *rptr; + int npts, nsave, 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); + + if (PyArray_As1D(&op2, (char **)&wsave, &nsave, PyArray_DOUBLE) == -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; + + 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; + } + 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 *self, PyObject *args) +{ + PyArrayObject *op; + int dim, n; + + if (!PyArg_ParseTuple(args, "i", &n)) return NULL; + + dim = 2*n+15; /*Magic size needed by rffti*/ + /*Create a 1 dimensional array of dimensions of type double*/ + op = (PyArrayObject *)PyArray_FromDims(1, &dim, PyArray_DOUBLE); + if (op == NULL) return NULL; + + rffti(n, (double *)((PyArrayObject*)op)->data); + + 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..4e1da2abd --- /dev/null +++ b/numpy/fft/helper.py @@ -0,0 +1,66 @@ +""" +Discrete Fourier Transforms - helper.py +""" +# Created by Pearu Peterson, September 2002 + +__all__ = ['fftshift','ifftshift','fftfreq'] + +from numpy.core import asarray, concatenate, arange, take, \ + array, integer +import types + +def fftshift(x,axes=None): + """ fftshift(x, axes=None) -> y + + Shift zero-frequency component to center of spectrum. + + This function swaps half-spaces for all axes listed (defaults to all). + + Notes: + If len(x) is even then the Nyquist component is y[0]. + """ + 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): + """ ifftshift(x,axes=None) - > y + + Inverse of 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): + """ fftfreq(n, d=1.0) -> f + + DFT 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 + """ + assert isinstance(n,types.IntType) or isinstance(n, integer) + k = range(0,(n-1)/2+1)+range(-(n/2),0) + return array(k,'d')/(n*d) diff --git a/numpy/fft/info.py b/numpy/fft/info.py new file mode 100644 index 000000000..4aa9b7d2c --- /dev/null +++ b/numpy/fft/info.py @@ -0,0 +1,29 @@ +"""\ +Core FFT routines +================== + + Standard FFTs + + fft + ifft + fft2 + ifft2 + fftn + ifftn + + Real FFTs + + refft + irefft + refft2 + irefft2 + refftn + irefftn + + Hermite FFTs + + hfft + ihfft +""" + +depends = ['core'] diff --git a/numpy/fft/old.py b/numpy/fft/old.py new file mode 100644 index 000000000..86723af95 --- /dev/null +++ b/numpy/fft/old.py @@ -0,0 +1,19 @@ + +__all__ = ['fft', 'fft2d', 'fftnd', 'hermite_fft', 'inverse_fft', 'inverse_fft2d', + 'inverse_fftnd', 'inverse_hermite_fft', 'inverse_real_fft', 'inverse_real_fft2d', + 'inverse_real_fftnd', 'real_fft', 'real_fft2d', 'real_fftnd'] + +from fftpack import fft +from fftpack import fft2 as fft2d +from fftpack import fftn as fftnd +from fftpack import hfft as hermite_fft +from fftpack import ifft as inverse_fft +from fftpack import ifft2 as inverse_fft2d +from fftpack import ifftn as inverse_fftnd +from fftpack import ihfft as inverse_hermite_fft +from fftpack import irfft as inverse_real_fft +from fftpack import irfft2 as inverse_real_fft2d +from fftpack import irfftn as inverse_real_fftnd +from fftpack import rfft as real_fft +from fftpack import rfft2 as real_fft2d +from fftpack import rfftn as real_fftnd diff --git a/numpy/fft/setup.py b/numpy/fft/setup.py new file mode 100644 index 000000000..4d1de069d --- /dev/null +++ b/numpy/fft/setup.py @@ -0,0 +1,20 @@ + +from os.path import join + +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(top_path='').todict()) diff --git a/numpy/fft/tests/test_helper.py b/numpy/fft/tests/test_helper.py new file mode 100644 index 000000000..3d02c01df --- /dev/null +++ b/numpy/fft/tests/test_helper.py @@ -0,0 +1,45 @@ +#!/usr/bin/env python +# Copied from fftpack.helper by Pearu Peterson, October 2005 +""" Test functions for fftpack.helper module +""" + +import sys +from numpy.testing import * +set_package_path() +from numpy.fft import fftshift,ifftshift,fftfreq +del sys.path[0] + +from numpy import pi + +def random(size): + return rand(*size) + +class test_fftshift(NumpyTestCase): + + def check_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 check_inverse(self): + for n in [1,4,9,100,211]: + x = random((n,)) + assert_array_almost_equal(ifftshift(fftshift(x)),x) + +class test_fftfreq(NumpyTestCase): + + def check_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__": + NumpyTest().run() |