diff options
Diffstat (limited to 'numpy/dft/fftpack.py')
| -rw-r--r-- | numpy/dft/fftpack.py | 323 |
1 files changed, 0 insertions, 323 deletions
diff --git a/numpy/dft/fftpack.py b/numpy/dft/fftpack.py deleted file mode 100644 index 1cb24f2b7..000000000 --- a/numpy/dft/fftpack.py +++ /dev/null @@ -1,323 +0,0 @@ -""" -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') |