demo + utils venv
This commit is contained in:
d3m1g0d
@@ -0,0 +1,114 @@
|
||||
"""
|
||||
==================================================
|
||||
Discrete Fourier transforms (:mod:`scipy.fftpack`)
|
||||
==================================================
|
||||
|
||||
Fast Fourier Transforms (FFTs)
|
||||
==============================
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
fft - Fast (discrete) Fourier Transform (FFT)
|
||||
ifft - Inverse FFT
|
||||
fft2 - Two dimensional FFT
|
||||
ifft2 - Two dimensional inverse FFT
|
||||
fftn - n-dimensional FFT
|
||||
ifftn - n-dimensional inverse FFT
|
||||
rfft - FFT of strictly real-valued sequence
|
||||
irfft - Inverse of rfft
|
||||
dct - Discrete cosine transform
|
||||
idct - Inverse discrete cosine transform
|
||||
dctn - n-dimensional Discrete cosine transform
|
||||
idctn - n-dimensional Inverse discrete cosine transform
|
||||
dst - Discrete sine transform
|
||||
idst - Inverse discrete sine transform
|
||||
dstn - n-dimensional Discrete sine transform
|
||||
idstn - n-dimensional Inverse discrete sine transform
|
||||
|
||||
Differential and pseudo-differential operators
|
||||
==============================================
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
diff - Differentiation and integration of periodic sequences
|
||||
tilbert - Tilbert transform: cs_diff(x,h,h)
|
||||
itilbert - Inverse Tilbert transform: sc_diff(x,h,h)
|
||||
hilbert - Hilbert transform: cs_diff(x,inf,inf)
|
||||
ihilbert - Inverse Hilbert transform: sc_diff(x,inf,inf)
|
||||
cs_diff - cosh/sinh pseudo-derivative of periodic sequences
|
||||
sc_diff - sinh/cosh pseudo-derivative of periodic sequences
|
||||
ss_diff - sinh/sinh pseudo-derivative of periodic sequences
|
||||
cc_diff - cosh/cosh pseudo-derivative of periodic sequences
|
||||
shift - Shift periodic sequences
|
||||
|
||||
Helper functions
|
||||
================
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
fftshift - Shift the zero-frequency component to the center of the spectrum
|
||||
ifftshift - The inverse of `fftshift`
|
||||
fftfreq - Return the Discrete Fourier Transform sample frequencies
|
||||
rfftfreq - DFT sample frequencies (for usage with rfft, irfft)
|
||||
next_fast_len - Find the optimal length to zero-pad an FFT for speed
|
||||
|
||||
Note that ``fftshift``, ``ifftshift`` and ``fftfreq`` are numpy functions
|
||||
exposed by ``fftpack``; importing them from ``numpy`` should be preferred.
|
||||
|
||||
Convolutions (:mod:`scipy.fftpack.convolve`)
|
||||
============================================
|
||||
|
||||
.. module:: scipy.fftpack.convolve
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
convolve
|
||||
convolve_z
|
||||
init_convolution_kernel
|
||||
destroy_convolve_cache
|
||||
|
||||
"""
|
||||
|
||||
# List of possibly useful functions in scipy.fftpack._fftpack:
|
||||
# drfft
|
||||
# zfft
|
||||
# zrfft
|
||||
# zfftnd
|
||||
# destroy_drfft_cache
|
||||
# destroy_zfft_cache
|
||||
# destroy_zfftnd_cache
|
||||
|
||||
from __future__ import division, print_function, absolute_import
|
||||
|
||||
|
||||
__all__ = ['fft','ifft','fftn','ifftn','rfft','irfft',
|
||||
'fft2','ifft2',
|
||||
'diff',
|
||||
'tilbert','itilbert','hilbert','ihilbert',
|
||||
'sc_diff','cs_diff','cc_diff','ss_diff',
|
||||
'shift',
|
||||
'fftfreq', 'rfftfreq',
|
||||
'fftshift', 'ifftshift',
|
||||
'next_fast_len',
|
||||
]
|
||||
|
||||
from .basic import *
|
||||
from .pseudo_diffs import *
|
||||
from .helper import *
|
||||
|
||||
from numpy.dual import register_func
|
||||
for k in ['fft', 'ifft', 'fftn', 'ifftn', 'fft2', 'ifft2']:
|
||||
register_func(k, eval(k))
|
||||
del k, register_func
|
||||
|
||||
from .realtransforms import *
|
||||
__all__.extend(['dct', 'idct', 'dst', 'idst', 'dctn', 'idctn', 'dstn',
|
||||
'idstn'])
|
||||
|
||||
from scipy._lib._testutils import PytestTester
|
||||
test = PytestTester(__name__)
|
||||
del PytestTester
|
||||
BIN
Binary file not shown.
Binary file not shown.
BIN
Binary file not shown.
BIN
Binary file not shown.
BIN
Binary file not shown.
Binary file not shown.
Executable
BIN
Binary file not shown.
@@ -0,0 +1,702 @@
|
||||
"""
|
||||
Discrete Fourier Transforms - basic.py
|
||||
"""
|
||||
# Created by Pearu Peterson, August,September 2002
|
||||
from __future__ import division, print_function, absolute_import
|
||||
|
||||
__all__ = ['fft','ifft','fftn','ifftn','rfft','irfft',
|
||||
'fft2','ifft2']
|
||||
|
||||
from numpy import swapaxes, zeros
|
||||
import numpy
|
||||
from . import _fftpack
|
||||
from scipy.fftpack.helper import _init_nd_shape_and_axes_sorted
|
||||
|
||||
import atexit
|
||||
atexit.register(_fftpack.destroy_zfft_cache)
|
||||
atexit.register(_fftpack.destroy_zfftnd_cache)
|
||||
atexit.register(_fftpack.destroy_drfft_cache)
|
||||
atexit.register(_fftpack.destroy_cfft_cache)
|
||||
atexit.register(_fftpack.destroy_cfftnd_cache)
|
||||
atexit.register(_fftpack.destroy_rfft_cache)
|
||||
del atexit
|
||||
|
||||
|
||||
def istype(arr, typeclass):
|
||||
return issubclass(arr.dtype.type, typeclass)
|
||||
|
||||
|
||||
def _datacopied(arr, original):
|
||||
"""
|
||||
Strict check for `arr` not sharing any data with `original`,
|
||||
under the assumption that arr = asarray(original)
|
||||
|
||||
"""
|
||||
if arr is original:
|
||||
return False
|
||||
if not isinstance(original, numpy.ndarray) and hasattr(original, '__array__'):
|
||||
return False
|
||||
return arr.base is None
|
||||
|
||||
# XXX: single precision FFTs partially disabled due to accuracy issues
|
||||
# for large prime-sized inputs.
|
||||
#
|
||||
# See http://permalink.gmane.org/gmane.comp.python.scientific.devel/13834
|
||||
# ("fftpack test failures for 0.8.0b1", Ralf Gommers, 17 Jun 2010,
|
||||
# @ scipy-dev)
|
||||
#
|
||||
# These should be re-enabled once the problems are resolved
|
||||
|
||||
|
||||
def _is_safe_size(n):
|
||||
"""
|
||||
Is the size of FFT such that FFTPACK can handle it in single precision
|
||||
with sufficient accuracy?
|
||||
|
||||
Composite numbers of 2, 3, and 5 are accepted, as FFTPACK has those
|
||||
"""
|
||||
n = int(n)
|
||||
|
||||
if n == 0:
|
||||
return True
|
||||
|
||||
# Divide by 3 until you can't, then by 5 until you can't
|
||||
for c in (3, 5):
|
||||
while n % c == 0:
|
||||
n //= c
|
||||
|
||||
# Return True if the remainder is a power of 2
|
||||
return not n & (n-1)
|
||||
|
||||
|
||||
def _fake_crfft(x, n, *a, **kw):
|
||||
if _is_safe_size(n):
|
||||
return _fftpack.crfft(x, n, *a, **kw)
|
||||
else:
|
||||
return _fftpack.zrfft(x, n, *a, **kw).astype(numpy.complex64)
|
||||
|
||||
|
||||
def _fake_cfft(x, n, *a, **kw):
|
||||
if _is_safe_size(n):
|
||||
return _fftpack.cfft(x, n, *a, **kw)
|
||||
else:
|
||||
return _fftpack.zfft(x, n, *a, **kw).astype(numpy.complex64)
|
||||
|
||||
|
||||
def _fake_rfft(x, n, *a, **kw):
|
||||
if _is_safe_size(n):
|
||||
return _fftpack.rfft(x, n, *a, **kw)
|
||||
else:
|
||||
return _fftpack.drfft(x, n, *a, **kw).astype(numpy.float32)
|
||||
|
||||
|
||||
def _fake_cfftnd(x, shape, *a, **kw):
|
||||
if numpy.all(list(map(_is_safe_size, shape))):
|
||||
return _fftpack.cfftnd(x, shape, *a, **kw)
|
||||
else:
|
||||
return _fftpack.zfftnd(x, shape, *a, **kw).astype(numpy.complex64)
|
||||
|
||||
|
||||
_DTYPE_TO_FFT = {
|
||||
# numpy.dtype(numpy.float32): _fftpack.crfft,
|
||||
numpy.dtype(numpy.float32): _fake_crfft,
|
||||
numpy.dtype(numpy.float64): _fftpack.zrfft,
|
||||
# numpy.dtype(numpy.complex64): _fftpack.cfft,
|
||||
numpy.dtype(numpy.complex64): _fake_cfft,
|
||||
numpy.dtype(numpy.complex128): _fftpack.zfft,
|
||||
}
|
||||
|
||||
_DTYPE_TO_RFFT = {
|
||||
# numpy.dtype(numpy.float32): _fftpack.rfft,
|
||||
numpy.dtype(numpy.float32): _fake_rfft,
|
||||
numpy.dtype(numpy.float64): _fftpack.drfft,
|
||||
}
|
||||
|
||||
_DTYPE_TO_FFTN = {
|
||||
# numpy.dtype(numpy.complex64): _fftpack.cfftnd,
|
||||
numpy.dtype(numpy.complex64): _fake_cfftnd,
|
||||
numpy.dtype(numpy.complex128): _fftpack.zfftnd,
|
||||
# numpy.dtype(numpy.float32): _fftpack.cfftnd,
|
||||
numpy.dtype(numpy.float32): _fake_cfftnd,
|
||||
numpy.dtype(numpy.float64): _fftpack.zfftnd,
|
||||
}
|
||||
|
||||
|
||||
def _asfarray(x):
|
||||
"""Like numpy asfarray, except that it does not modify x dtype if x is
|
||||
already an array with a float dtype, and do not cast complex types to
|
||||
real."""
|
||||
if hasattr(x, "dtype") and x.dtype.char in numpy.typecodes["AllFloat"]:
|
||||
# 'dtype' attribute does not ensure that the
|
||||
# object is an ndarray (e.g. Series class
|
||||
# from the pandas library)
|
||||
if x.dtype == numpy.half:
|
||||
# no half-precision routines, so convert to single precision
|
||||
return numpy.asarray(x, dtype=numpy.float32)
|
||||
return numpy.asarray(x, dtype=x.dtype)
|
||||
else:
|
||||
# We cannot use asfarray directly because it converts sequences of
|
||||
# complex to sequence of real
|
||||
ret = numpy.asarray(x)
|
||||
if ret.dtype == numpy.half:
|
||||
return numpy.asarray(ret, dtype=numpy.float32)
|
||||
elif ret.dtype.char not in numpy.typecodes["AllFloat"]:
|
||||
return numpy.asfarray(x)
|
||||
return ret
|
||||
|
||||
|
||||
def _fix_shape(x, n, axis):
|
||||
""" Internal auxiliary function for _raw_fft, _raw_fftnd."""
|
||||
s = list(x.shape)
|
||||
if s[axis] > n:
|
||||
index = [slice(None)]*len(s)
|
||||
index[axis] = slice(0,n)
|
||||
x = x[tuple(index)]
|
||||
return x, False
|
||||
else:
|
||||
index = [slice(None)]*len(s)
|
||||
index[axis] = slice(0,s[axis])
|
||||
s[axis] = n
|
||||
z = zeros(s,x.dtype.char)
|
||||
z[tuple(index)] = x
|
||||
return z, True
|
||||
|
||||
|
||||
def _raw_fft(x, n, axis, direction, overwrite_x, work_function):
|
||||
""" Internal auxiliary function for fft, ifft, rfft, irfft."""
|
||||
if n is None:
|
||||
n = x.shape[axis]
|
||||
elif n != x.shape[axis]:
|
||||
x, copy_made = _fix_shape(x,n,axis)
|
||||
overwrite_x = overwrite_x or copy_made
|
||||
|
||||
if n < 1:
|
||||
raise ValueError("Invalid number of FFT data points "
|
||||
"(%d) specified." % n)
|
||||
|
||||
if axis == -1 or axis == len(x.shape)-1:
|
||||
r = work_function(x,n,direction,overwrite_x=overwrite_x)
|
||||
else:
|
||||
x = swapaxes(x, axis, -1)
|
||||
r = work_function(x,n,direction,overwrite_x=overwrite_x)
|
||||
r = swapaxes(r, axis, -1)
|
||||
return r
|
||||
|
||||
|
||||
def fft(x, n=None, axis=-1, overwrite_x=False):
|
||||
"""
|
||||
Return discrete Fourier transform of real or complex sequence.
|
||||
|
||||
The returned complex array contains ``y(0), y(1),..., y(n-1)`` where
|
||||
|
||||
``y(j) = (x * exp(-2*pi*sqrt(-1)*j*np.arange(n)/n)).sum()``.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
Array to Fourier transform.
|
||||
n : int, optional
|
||||
Length of the Fourier transform. If ``n < x.shape[axis]``, `x` is
|
||||
truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
|
||||
default results in ``n = x.shape[axis]``.
|
||||
axis : int, optional
|
||||
Axis along which the fft's are computed; the default is over the
|
||||
last axis (i.e., ``axis=-1``).
|
||||
overwrite_x : bool, optional
|
||||
If True, the contents of `x` can be destroyed; the default is False.
|
||||
|
||||
Returns
|
||||
-------
|
||||
z : complex ndarray
|
||||
with the elements::
|
||||
|
||||
[y(0),y(1),..,y(n/2),y(1-n/2),...,y(-1)] if n is even
|
||||
[y(0),y(1),..,y((n-1)/2),y(-(n-1)/2),...,y(-1)] if n is odd
|
||||
|
||||
where::
|
||||
|
||||
y(j) = sum[k=0..n-1] x[k] * exp(-sqrt(-1)*j*k* 2*pi/n), j = 0..n-1
|
||||
|
||||
See Also
|
||||
--------
|
||||
ifft : Inverse FFT
|
||||
rfft : FFT of a real sequence
|
||||
|
||||
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]`` contains the
|
||||
positive-frequency terms, and ``A[n/2:]`` 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].
|
||||
To rearrange the fft output so that the zero-frequency component is
|
||||
centered, like [-4, -3, -2, -1, 0, 1, 2, 3], use `fftshift`.
|
||||
|
||||
Both single and double precision routines are implemented. Half precision
|
||||
inputs will be converted to single precision. Non floating-point inputs
|
||||
will be converted to double precision. Long-double precision inputs are
|
||||
not supported.
|
||||
|
||||
This function is most efficient when `n` is a power of two, and least
|
||||
efficient when `n` is prime.
|
||||
|
||||
Note that if ``x`` is real-valued then ``A[j] == A[n-j].conjugate()``.
|
||||
If ``x`` is real-valued and ``n`` is even then ``A[n/2]`` is real.
|
||||
|
||||
If the data type of `x` is real, a "real FFT" algorithm is automatically
|
||||
used, which roughly halves the computation time. To increase efficiency
|
||||
a little further, use `rfft`, which does the same calculation, but only
|
||||
outputs half of the symmetrical spectrum. If the data is both real and
|
||||
symmetrical, the `dct` can again double the efficiency, by generating
|
||||
half of the spectrum from half of the signal.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from scipy.fftpack import fft, ifft
|
||||
>>> x = np.arange(5)
|
||||
>>> np.allclose(fft(ifft(x)), x, atol=1e-15) # within numerical accuracy.
|
||||
True
|
||||
|
||||
"""
|
||||
tmp = _asfarray(x)
|
||||
|
||||
try:
|
||||
work_function = _DTYPE_TO_FFT[tmp.dtype]
|
||||
except KeyError:
|
||||
raise ValueError("type %s is not supported" % tmp.dtype)
|
||||
|
||||
if not (istype(tmp, numpy.complex64) or istype(tmp, numpy.complex128)):
|
||||
overwrite_x = 1
|
||||
|
||||
overwrite_x = overwrite_x or _datacopied(tmp, x)
|
||||
|
||||
if n is None:
|
||||
n = tmp.shape[axis]
|
||||
elif n != tmp.shape[axis]:
|
||||
tmp, copy_made = _fix_shape(tmp,n,axis)
|
||||
overwrite_x = overwrite_x or copy_made
|
||||
|
||||
if n < 1:
|
||||
raise ValueError("Invalid number of FFT data points "
|
||||
"(%d) specified." % n)
|
||||
|
||||
if axis == -1 or axis == len(tmp.shape) - 1:
|
||||
return work_function(tmp,n,1,0,overwrite_x)
|
||||
|
||||
tmp = swapaxes(tmp, axis, -1)
|
||||
tmp = work_function(tmp,n,1,0,overwrite_x)
|
||||
return swapaxes(tmp, axis, -1)
|
||||
|
||||
|
||||
def ifft(x, n=None, axis=-1, overwrite_x=False):
|
||||
"""
|
||||
Return discrete inverse Fourier transform of real or complex sequence.
|
||||
|
||||
The returned complex array contains ``y(0), y(1),..., y(n-1)`` where
|
||||
|
||||
``y(j) = (x * exp(2*pi*sqrt(-1)*j*np.arange(n)/n)).mean()``.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
Transformed data to invert.
|
||||
n : int, optional
|
||||
Length of the inverse Fourier transform. If ``n < x.shape[axis]``,
|
||||
`x` is truncated. If ``n > x.shape[axis]``, `x` is zero-padded.
|
||||
The default results in ``n = x.shape[axis]``.
|
||||
axis : int, optional
|
||||
Axis along which the ifft's are computed; the default is over the
|
||||
last axis (i.e., ``axis=-1``).
|
||||
overwrite_x : bool, optional
|
||||
If True, the contents of `x` can be destroyed; the default is False.
|
||||
|
||||
Returns
|
||||
-------
|
||||
ifft : ndarray of floats
|
||||
The inverse discrete Fourier transform.
|
||||
|
||||
See Also
|
||||
--------
|
||||
fft : Forward FFT
|
||||
|
||||
Notes
|
||||
-----
|
||||
Both single and double precision routines are implemented. Half precision
|
||||
inputs will be converted to single precision. Non floating-point inputs
|
||||
will be converted to double precision. Long-double precision inputs are
|
||||
not supported.
|
||||
|
||||
This function is most efficient when `n` is a power of two, and least
|
||||
efficient when `n` is prime.
|
||||
|
||||
If the data type of `x` is real, a "real IFFT" algorithm is automatically
|
||||
used, which roughly halves the computation time.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from scipy.fftpack import fft, ifft
|
||||
>>> import numpy as np
|
||||
>>> x = np.arange(5)
|
||||
>>> np.allclose(ifft(fft(x)), x, atol=1e-15) # within numerical accuracy.
|
||||
True
|
||||
|
||||
"""
|
||||
tmp = _asfarray(x)
|
||||
|
||||
try:
|
||||
work_function = _DTYPE_TO_FFT[tmp.dtype]
|
||||
except KeyError:
|
||||
raise ValueError("type %s is not supported" % tmp.dtype)
|
||||
|
||||
if not (istype(tmp, numpy.complex64) or istype(tmp, numpy.complex128)):
|
||||
overwrite_x = 1
|
||||
|
||||
overwrite_x = overwrite_x or _datacopied(tmp, x)
|
||||
|
||||
if n is None:
|
||||
n = tmp.shape[axis]
|
||||
elif n != tmp.shape[axis]:
|
||||
tmp, copy_made = _fix_shape(tmp,n,axis)
|
||||
overwrite_x = overwrite_x or copy_made
|
||||
|
||||
if n < 1:
|
||||
raise ValueError("Invalid number of FFT data points "
|
||||
"(%d) specified." % n)
|
||||
|
||||
if axis == -1 or axis == len(tmp.shape) - 1:
|
||||
return work_function(tmp,n,-1,1,overwrite_x)
|
||||
|
||||
tmp = swapaxes(tmp, axis, -1)
|
||||
tmp = work_function(tmp,n,-1,1,overwrite_x)
|
||||
return swapaxes(tmp, axis, -1)
|
||||
|
||||
|
||||
def rfft(x, n=None, axis=-1, overwrite_x=False):
|
||||
"""
|
||||
Discrete Fourier transform of a real sequence.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like, real-valued
|
||||
The data to transform.
|
||||
n : int, optional
|
||||
Defines the length of the Fourier transform. If `n` is not specified
|
||||
(the default) then ``n = x.shape[axis]``. If ``n < x.shape[axis]``,
|
||||
`x` is truncated, if ``n > x.shape[axis]``, `x` is zero-padded.
|
||||
axis : int, optional
|
||||
The axis along which the transform is applied. The default is the
|
||||
last axis.
|
||||
overwrite_x : bool, optional
|
||||
If set to true, the contents of `x` can be overwritten. Default is
|
||||
False.
|
||||
|
||||
Returns
|
||||
-------
|
||||
z : real ndarray
|
||||
The returned real array contains::
|
||||
|
||||
[y(0),Re(y(1)),Im(y(1)),...,Re(y(n/2))] if n is even
|
||||
[y(0),Re(y(1)),Im(y(1)),...,Re(y(n/2)),Im(y(n/2))] if n is odd
|
||||
|
||||
where::
|
||||
|
||||
y(j) = sum[k=0..n-1] x[k] * exp(-sqrt(-1)*j*k*2*pi/n)
|
||||
j = 0..n-1
|
||||
|
||||
See Also
|
||||
--------
|
||||
fft, irfft, numpy.fft.rfft
|
||||
|
||||
Notes
|
||||
-----
|
||||
Within numerical accuracy, ``y == rfft(irfft(y))``.
|
||||
|
||||
Both single and double precision routines are implemented. Half precision
|
||||
inputs will be converted to single precision. Non floating-point inputs
|
||||
will be converted to double precision. Long-double precision inputs are
|
||||
not supported.
|
||||
|
||||
To get an output with a complex datatype, consider using the related
|
||||
function `numpy.fft.rfft`.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from scipy.fftpack import fft, rfft
|
||||
>>> a = [9, -9, 1, 3]
|
||||
>>> fft(a)
|
||||
array([ 4. +0.j, 8.+12.j, 16. +0.j, 8.-12.j])
|
||||
>>> rfft(a)
|
||||
array([ 4., 8., 12., 16.])
|
||||
|
||||
"""
|
||||
tmp = _asfarray(x)
|
||||
|
||||
if not numpy.isrealobj(tmp):
|
||||
raise TypeError("1st argument must be real sequence")
|
||||
|
||||
try:
|
||||
work_function = _DTYPE_TO_RFFT[tmp.dtype]
|
||||
except KeyError:
|
||||
raise ValueError("type %s is not supported" % tmp.dtype)
|
||||
|
||||
overwrite_x = overwrite_x or _datacopied(tmp, x)
|
||||
|
||||
return _raw_fft(tmp,n,axis,1,overwrite_x,work_function)
|
||||
|
||||
|
||||
def irfft(x, n=None, axis=-1, overwrite_x=False):
|
||||
"""
|
||||
Return inverse discrete Fourier transform of real sequence x.
|
||||
|
||||
The contents of `x` are interpreted as the output of the `rfft`
|
||||
function.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
Transformed data to invert.
|
||||
n : int, optional
|
||||
Length of the inverse Fourier transform.
|
||||
If n < x.shape[axis], x is truncated.
|
||||
If n > x.shape[axis], x is zero-padded.
|
||||
The default results in n = x.shape[axis].
|
||||
axis : int, optional
|
||||
Axis along which the ifft's are computed; the default is over
|
||||
the last axis (i.e., axis=-1).
|
||||
overwrite_x : bool, optional
|
||||
If True, the contents of `x` can be destroyed; the default is False.
|
||||
|
||||
Returns
|
||||
-------
|
||||
irfft : ndarray of floats
|
||||
The inverse discrete Fourier transform.
|
||||
|
||||
See Also
|
||||
--------
|
||||
rfft, ifft, numpy.fft.irfft
|
||||
|
||||
Notes
|
||||
-----
|
||||
The returned real array contains::
|
||||
|
||||
[y(0),y(1),...,y(n-1)]
|
||||
|
||||
where for n is even::
|
||||
|
||||
y(j) = 1/n (sum[k=1..n/2-1] (x[2*k-1]+sqrt(-1)*x[2*k])
|
||||
* exp(sqrt(-1)*j*k* 2*pi/n)
|
||||
+ c.c. + x[0] + (-1)**(j) x[n-1])
|
||||
|
||||
and for n is odd::
|
||||
|
||||
y(j) = 1/n (sum[k=1..(n-1)/2] (x[2*k-1]+sqrt(-1)*x[2*k])
|
||||
* exp(sqrt(-1)*j*k* 2*pi/n)
|
||||
+ c.c. + x[0])
|
||||
|
||||
c.c. denotes complex conjugate of preceding expression.
|
||||
|
||||
For details on input parameters, see `rfft`.
|
||||
|
||||
To process (conjugate-symmetric) frequency-domain data with a complex
|
||||
datatype, consider using the related function `numpy.fft.irfft`.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from scipy.fftpack import rfft, irfft
|
||||
>>> a = [1.0, 2.0, 3.0, 4.0, 5.0]
|
||||
>>> irfft(a)
|
||||
array([ 2.6 , -3.16405192, 1.24398433, -1.14955713, 1.46962473])
|
||||
>>> irfft(rfft(a))
|
||||
array([1., 2., 3., 4., 5.])
|
||||
|
||||
"""
|
||||
tmp = _asfarray(x)
|
||||
if not numpy.isrealobj(tmp):
|
||||
raise TypeError("1st argument must be real sequence")
|
||||
|
||||
try:
|
||||
work_function = _DTYPE_TO_RFFT[tmp.dtype]
|
||||
except KeyError:
|
||||
raise ValueError("type %s is not supported" % tmp.dtype)
|
||||
|
||||
overwrite_x = overwrite_x or _datacopied(tmp, x)
|
||||
|
||||
return _raw_fft(tmp,n,axis,-1,overwrite_x,work_function)
|
||||
|
||||
|
||||
def _raw_fftnd(x, s, axes, direction, overwrite_x, work_function):
|
||||
"""Internal auxiliary function for fftnd, ifftnd."""
|
||||
noaxes = axes is None
|
||||
s, axes = _init_nd_shape_and_axes_sorted(x, s, axes)
|
||||
|
||||
# No need to swap axes, array is in C order
|
||||
if noaxes:
|
||||
for ax in axes:
|
||||
x, copy_made = _fix_shape(x, s[ax], ax)
|
||||
overwrite_x = overwrite_x or copy_made
|
||||
return work_function(x, s, direction, overwrite_x=overwrite_x)
|
||||
|
||||
# Swap the request axes, last first (i.e. First swap the axis which ends up
|
||||
# at -1, then at -2, etc...), such as the request axes on which the
|
||||
# operation is carried become the last ones
|
||||
for i in range(1, axes.size+1):
|
||||
x = numpy.swapaxes(x, axes[-i], -i)
|
||||
|
||||
# We can now operate on the axes waxes, the p last axes (p = len(axes)), by
|
||||
# fixing the shape of the input array to 1 for any axis the fft is not
|
||||
# carried upon.
|
||||
waxes = list(range(x.ndim - axes.size, x.ndim))
|
||||
shape = numpy.ones(x.ndim)
|
||||
shape[waxes] = s
|
||||
|
||||
for i in range(len(waxes)):
|
||||
x, copy_made = _fix_shape(x, s[i], waxes[i])
|
||||
overwrite_x = overwrite_x or copy_made
|
||||
|
||||
r = work_function(x, shape, direction, overwrite_x=overwrite_x)
|
||||
|
||||
# reswap in the reverse order (first axis first, etc...) to get original
|
||||
# order
|
||||
for i in range(len(axes), 0, -1):
|
||||
r = numpy.swapaxes(r, -i, axes[-i])
|
||||
|
||||
return r
|
||||
|
||||
|
||||
def fftn(x, shape=None, axes=None, overwrite_x=False):
|
||||
"""
|
||||
Return multidimensional discrete Fourier transform.
|
||||
|
||||
The returned array contains::
|
||||
|
||||
y[j_1,..,j_d] = sum[k_1=0..n_1-1, ..., k_d=0..n_d-1]
|
||||
x[k_1,..,k_d] * prod[i=1..d] exp(-sqrt(-1)*2*pi/n_i * j_i * k_i)
|
||||
|
||||
where d = len(x.shape) and n = x.shape.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The (n-dimensional) array to transform.
|
||||
shape : int or array_like of ints or None, optional
|
||||
The shape of the result. If both `shape` and `axes` (see below) are
|
||||
None, `shape` is ``x.shape``; if `shape` is None but `axes` is
|
||||
not None, then `shape` is ``scipy.take(x.shape, axes, axis=0)``.
|
||||
If ``shape[i] > x.shape[i]``, the i-th dimension is padded with zeros.
|
||||
If ``shape[i] < x.shape[i]``, the i-th dimension is truncated to
|
||||
length ``shape[i]``.
|
||||
If any element of `shape` is -1, the size of the corresponding
|
||||
dimension of `x` is used.
|
||||
axes : int or array_like of ints or None, optional
|
||||
The axes of `x` (`y` if `shape` is not None) along which the
|
||||
transform is applied.
|
||||
The default is over all axes.
|
||||
overwrite_x : bool, optional
|
||||
If True, the contents of `x` can be destroyed. Default is False.
|
||||
|
||||
Returns
|
||||
-------
|
||||
y : complex-valued n-dimensional numpy array
|
||||
The (n-dimensional) DFT of the input array.
|
||||
|
||||
See Also
|
||||
--------
|
||||
ifftn
|
||||
|
||||
Notes
|
||||
-----
|
||||
If ``x`` is real-valued, then
|
||||
``y[..., j_i, ...] == y[..., n_i-j_i, ...].conjugate()``.
|
||||
|
||||
Both single and double precision routines are implemented. Half precision
|
||||
inputs will be converted to single precision. Non floating-point inputs
|
||||
will be converted to double precision. Long-double precision inputs are
|
||||
not supported.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from scipy.fftpack import fftn, ifftn
|
||||
>>> y = (-np.arange(16), 8 - np.arange(16), np.arange(16))
|
||||
>>> np.allclose(y, fftn(ifftn(y)))
|
||||
True
|
||||
|
||||
"""
|
||||
return _raw_fftn_dispatch(x, shape, axes, overwrite_x, 1)
|
||||
|
||||
|
||||
def _raw_fftn_dispatch(x, shape, axes, overwrite_x, direction):
|
||||
tmp = _asfarray(x)
|
||||
|
||||
try:
|
||||
work_function = _DTYPE_TO_FFTN[tmp.dtype]
|
||||
except KeyError:
|
||||
raise ValueError("type %s is not supported" % tmp.dtype)
|
||||
|
||||
if not (istype(tmp, numpy.complex64) or istype(tmp, numpy.complex128)):
|
||||
overwrite_x = 1
|
||||
|
||||
overwrite_x = overwrite_x or _datacopied(tmp, x)
|
||||
return _raw_fftnd(tmp, shape, axes, direction, overwrite_x, work_function)
|
||||
|
||||
|
||||
def ifftn(x, shape=None, axes=None, overwrite_x=False):
|
||||
"""
|
||||
Return inverse multi-dimensional discrete Fourier transform.
|
||||
|
||||
The sequence can be of an arbitrary type.
|
||||
|
||||
The returned array contains::
|
||||
|
||||
y[j_1,..,j_d] = 1/p * sum[k_1=0..n_1-1, ..., k_d=0..n_d-1]
|
||||
x[k_1,..,k_d] * prod[i=1..d] exp(sqrt(-1)*2*pi/n_i * j_i * k_i)
|
||||
|
||||
where ``d = len(x.shape)``, ``n = x.shape``, and ``p = prod[i=1..d] n_i``.
|
||||
|
||||
For description of parameters see `fftn`.
|
||||
|
||||
See Also
|
||||
--------
|
||||
fftn : for detailed information.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from scipy.fftpack import fftn, ifftn
|
||||
>>> import numpy as np
|
||||
>>> y = (-np.arange(16), 8 - np.arange(16), np.arange(16))
|
||||
>>> np.allclose(y, ifftn(fftn(y)))
|
||||
True
|
||||
|
||||
"""
|
||||
return _raw_fftn_dispatch(x, shape, axes, overwrite_x, -1)
|
||||
|
||||
|
||||
def fft2(x, shape=None, axes=(-2,-1), overwrite_x=False):
|
||||
"""
|
||||
2-D discrete Fourier transform.
|
||||
|
||||
Return the two-dimensional discrete Fourier transform of the 2-D argument
|
||||
`x`.
|
||||
|
||||
See Also
|
||||
--------
|
||||
fftn : for detailed information.
|
||||
|
||||
"""
|
||||
return fftn(x,shape,axes,overwrite_x)
|
||||
|
||||
|
||||
def ifft2(x, shape=None, axes=(-2,-1), overwrite_x=False):
|
||||
"""
|
||||
2-D discrete inverse Fourier transform of real or complex sequence.
|
||||
|
||||
Return inverse two-dimensional discrete Fourier transform of
|
||||
arbitrary type sequence x.
|
||||
|
||||
See `ifft` for more information.
|
||||
|
||||
See also
|
||||
--------
|
||||
fft2, ifft
|
||||
|
||||
"""
|
||||
return ifftn(x,shape,axes,overwrite_x)
|
||||
Executable
BIN
Binary file not shown.
@@ -0,0 +1,277 @@
|
||||
from __future__ import division, print_function, absolute_import
|
||||
|
||||
import operator
|
||||
from numpy import (arange, array, asarray, atleast_1d, intc, integer,
|
||||
isscalar, issubdtype, take, unique, where)
|
||||
from numpy.fft.helper import fftshift, ifftshift, fftfreq
|
||||
from bisect import bisect_left
|
||||
|
||||
__all__ = ['fftshift', 'ifftshift', 'fftfreq', 'rfftfreq', 'next_fast_len']
|
||||
|
||||
|
||||
def rfftfreq(n, d=1.0):
|
||||
"""DFT sample frequencies (for usage with rfft, irfft).
|
||||
|
||||
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,1,2,2,...,n/2-1,n/2-1,n/2]/(d*n) if n is even
|
||||
f = [0,1,1,2,2,...,n/2-1,n/2-1,n/2,n/2]/(d*n) if n is odd
|
||||
|
||||
Parameters
|
||||
----------
|
||||
n : int
|
||||
Window length.
|
||||
d : scalar, optional
|
||||
Sample spacing. Default is 1.
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : ndarray
|
||||
The array of length `n`, containing the sample frequencies.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from scipy import fftpack
|
||||
>>> sig = np.array([-2, 8, 6, 4, 1, 0, 3, 5], dtype=float)
|
||||
>>> sig_fft = fftpack.rfft(sig)
|
||||
>>> n = sig_fft.size
|
||||
>>> timestep = 0.1
|
||||
>>> freq = fftpack.rfftfreq(n, d=timestep)
|
||||
>>> freq
|
||||
array([ 0. , 1.25, 1.25, 2.5 , 2.5 , 3.75, 3.75, 5. ])
|
||||
|
||||
"""
|
||||
n = operator.index(n)
|
||||
if n < 0:
|
||||
raise ValueError("n = %s is not valid. "
|
||||
"n must be a nonnegative integer." % n)
|
||||
|
||||
return (arange(1, n + 1, dtype=int) // 2) / float(n * d)
|
||||
|
||||
|
||||
def next_fast_len(target):
|
||||
"""
|
||||
Find the next fast size of input data to `fft`, for zero-padding, etc.
|
||||
|
||||
SciPy's FFTPACK has efficient functions for radix {2, 3, 4, 5}, so this
|
||||
returns the next composite of the prime factors 2, 3, and 5 which is
|
||||
greater than or equal to `target`. (These are also known as 5-smooth
|
||||
numbers, regular numbers, or Hamming numbers.)
|
||||
|
||||
Parameters
|
||||
----------
|
||||
target : int
|
||||
Length to start searching from. Must be a positive integer.
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : int
|
||||
The first 5-smooth number greater than or equal to `target`.
|
||||
|
||||
Notes
|
||||
-----
|
||||
.. versionadded:: 0.18.0
|
||||
|
||||
Examples
|
||||
--------
|
||||
On a particular machine, an FFT of prime length takes 133 ms:
|
||||
|
||||
>>> from scipy import fftpack
|
||||
>>> min_len = 10007 # prime length is worst case for speed
|
||||
>>> a = np.random.randn(min_len)
|
||||
>>> b = fftpack.fft(a)
|
||||
|
||||
Zero-padding to the next 5-smooth length reduces computation time to
|
||||
211 us, a speedup of 630 times:
|
||||
|
||||
>>> fftpack.helper.next_fast_len(min_len)
|
||||
10125
|
||||
>>> b = fftpack.fft(a, 10125)
|
||||
|
||||
Rounding up to the next power of 2 is not optimal, taking 367 us to
|
||||
compute, 1.7 times as long as the 5-smooth size:
|
||||
|
||||
>>> b = fftpack.fft(a, 16384)
|
||||
|
||||
"""
|
||||
hams = (8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48,
|
||||
50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, 108, 120, 125, 128,
|
||||
135, 144, 150, 160, 162, 180, 192, 200, 216, 225, 240, 243, 250,
|
||||
256, 270, 288, 300, 320, 324, 360, 375, 384, 400, 405, 432, 450,
|
||||
480, 486, 500, 512, 540, 576, 600, 625, 640, 648, 675, 720, 729,
|
||||
750, 768, 800, 810, 864, 900, 960, 972, 1000, 1024, 1080, 1125,
|
||||
1152, 1200, 1215, 1250, 1280, 1296, 1350, 1440, 1458, 1500, 1536,
|
||||
1600, 1620, 1728, 1800, 1875, 1920, 1944, 2000, 2025, 2048, 2160,
|
||||
2187, 2250, 2304, 2400, 2430, 2500, 2560, 2592, 2700, 2880, 2916,
|
||||
3000, 3072, 3125, 3200, 3240, 3375, 3456, 3600, 3645, 3750, 3840,
|
||||
3888, 4000, 4050, 4096, 4320, 4374, 4500, 4608, 4800, 4860, 5000,
|
||||
5120, 5184, 5400, 5625, 5760, 5832, 6000, 6075, 6144, 6250, 6400,
|
||||
6480, 6561, 6750, 6912, 7200, 7290, 7500, 7680, 7776, 8000, 8100,
|
||||
8192, 8640, 8748, 9000, 9216, 9375, 9600, 9720, 10000)
|
||||
|
||||
target = int(target)
|
||||
|
||||
if target <= 6:
|
||||
return target
|
||||
|
||||
# Quickly check if it's already a power of 2
|
||||
if not (target & (target-1)):
|
||||
return target
|
||||
|
||||
# Get result quickly for small sizes, since FFT itself is similarly fast.
|
||||
if target <= hams[-1]:
|
||||
return hams[bisect_left(hams, target)]
|
||||
|
||||
match = float('inf') # Anything found will be smaller
|
||||
p5 = 1
|
||||
while p5 < target:
|
||||
p35 = p5
|
||||
while p35 < target:
|
||||
# Ceiling integer division, avoiding conversion to float
|
||||
# (quotient = ceil(target / p35))
|
||||
quotient = -(-target // p35)
|
||||
|
||||
# Quickly find next power of 2 >= quotient
|
||||
p2 = 2**((quotient - 1).bit_length())
|
||||
|
||||
N = p2 * p35
|
||||
if N == target:
|
||||
return N
|
||||
elif N < match:
|
||||
match = N
|
||||
p35 *= 3
|
||||
if p35 == target:
|
||||
return p35
|
||||
if p35 < match:
|
||||
match = p35
|
||||
p5 *= 5
|
||||
if p5 == target:
|
||||
return p5
|
||||
if p5 < match:
|
||||
match = p5
|
||||
return match
|
||||
|
||||
|
||||
def _init_nd_shape_and_axes(x, shape, axes):
|
||||
"""Handle shape and axes arguments for n-dimensional transforms.
|
||||
|
||||
Returns the shape and axes in a standard form, taking into account negative
|
||||
values and checking for various potential errors.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The input array.
|
||||
shape : int or array_like of ints or None
|
||||
The shape of the result. If both `shape` and `axes` (see below) are
|
||||
None, `shape` is ``x.shape``; if `shape` is None but `axes` is
|
||||
not None, then `shape` is ``scipy.take(x.shape, axes, axis=0)``.
|
||||
If `shape` is -1, the size of the corresponding dimension of `x` is
|
||||
used.
|
||||
axes : int or array_like of ints or None
|
||||
Axes along which the calculation is computed.
|
||||
The default is over all axes.
|
||||
Negative indices are automatically converted to their positive
|
||||
counterpart.
|
||||
|
||||
Returns
|
||||
-------
|
||||
shape : array
|
||||
The shape of the result. It is a 1D integer array.
|
||||
axes : array
|
||||
The shape of the result. It is a 1D integer array.
|
||||
|
||||
"""
|
||||
x = asarray(x)
|
||||
noshape = shape is None
|
||||
noaxes = axes is None
|
||||
|
||||
if noaxes:
|
||||
axes = arange(x.ndim, dtype=intc)
|
||||
else:
|
||||
axes = atleast_1d(axes)
|
||||
|
||||
if axes.size == 0:
|
||||
axes = axes.astype(intc)
|
||||
|
||||
if not axes.ndim == 1:
|
||||
raise ValueError("when given, axes values must be a scalar or vector")
|
||||
if not issubdtype(axes.dtype, integer):
|
||||
raise ValueError("when given, axes values must be integers")
|
||||
|
||||
axes = where(axes < 0, axes + x.ndim, axes)
|
||||
|
||||
if axes.size != 0 and (axes.max() >= x.ndim or axes.min() < 0):
|
||||
raise ValueError("axes exceeds dimensionality of input")
|
||||
if axes.size != 0 and unique(axes).shape != axes.shape:
|
||||
raise ValueError("all axes must be unique")
|
||||
|
||||
if not noshape:
|
||||
shape = atleast_1d(shape)
|
||||
elif isscalar(x):
|
||||
shape = array([], dtype=intc)
|
||||
elif noaxes:
|
||||
shape = array(x.shape, dtype=intc)
|
||||
else:
|
||||
shape = take(x.shape, axes)
|
||||
|
||||
if shape.size == 0:
|
||||
shape = shape.astype(intc)
|
||||
|
||||
if shape.ndim != 1:
|
||||
raise ValueError("when given, shape values must be a scalar or vector")
|
||||
if not issubdtype(shape.dtype, integer):
|
||||
raise ValueError("when given, shape values must be integers")
|
||||
if axes.shape != shape.shape:
|
||||
raise ValueError("when given, axes and shape arguments"
|
||||
" have to be of the same length")
|
||||
|
||||
shape = where(shape == -1, array(x.shape)[axes], shape)
|
||||
|
||||
if shape.size != 0 and (shape < 1).any():
|
||||
raise ValueError(
|
||||
"invalid number of data points ({0}) specified".format(shape))
|
||||
|
||||
return shape, axes
|
||||
|
||||
|
||||
def _init_nd_shape_and_axes_sorted(x, shape, axes):
|
||||
"""Handle and sort shape and axes arguments for n-dimensional transforms.
|
||||
|
||||
This is identical to `_init_nd_shape_and_axes`, except the axes are
|
||||
returned in sorted order and the shape is reordered to match.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The input array.
|
||||
shape : int or array_like of ints or None
|
||||
The shape of the result. If both `shape` and `axes` (see below) are
|
||||
None, `shape` is ``x.shape``; if `shape` is None but `axes` is
|
||||
not None, then `shape` is ``scipy.take(x.shape, axes, axis=0)``.
|
||||
If `shape` is -1, the size of the corresponding dimension of `x` is
|
||||
used.
|
||||
axes : int or array_like of ints or None
|
||||
Axes along which the calculation is computed.
|
||||
The default is over all axes.
|
||||
Negative indices are automatically converted to their positive
|
||||
counterpart.
|
||||
|
||||
Returns
|
||||
-------
|
||||
shape : array
|
||||
The shape of the result. It is a 1D integer array.
|
||||
axes : array
|
||||
The shape of the result. It is a 1D integer array.
|
||||
|
||||
"""
|
||||
noaxes = axes is None
|
||||
shape, axes = _init_nd_shape_and_axes(x, shape, axes)
|
||||
|
||||
if not noaxes:
|
||||
shape = shape[axes.argsort()]
|
||||
axes.sort()
|
||||
|
||||
return shape, axes
|
||||
@@ -0,0 +1,557 @@
|
||||
"""
|
||||
Differential and pseudo-differential operators.
|
||||
"""
|
||||
# Created by Pearu Peterson, September 2002
|
||||
from __future__ import division, print_function, absolute_import
|
||||
|
||||
|
||||
__all__ = ['diff',
|
||||
'tilbert','itilbert','hilbert','ihilbert',
|
||||
'cs_diff','cc_diff','sc_diff','ss_diff',
|
||||
'shift']
|
||||
|
||||
from numpy import pi, asarray, sin, cos, sinh, cosh, tanh, iscomplexobj
|
||||
from . import convolve
|
||||
|
||||
from scipy.fftpack.basic import _datacopied
|
||||
|
||||
import atexit
|
||||
atexit.register(convolve.destroy_convolve_cache)
|
||||
del atexit
|
||||
|
||||
|
||||
_cache = {}
|
||||
|
||||
|
||||
def diff(x,order=1,period=None, _cache=_cache):
|
||||
"""
|
||||
Return k-th derivative (or integral) of a periodic sequence x.
|
||||
|
||||
If x_j and y_j are Fourier coefficients of periodic functions x
|
||||
and y, respectively, then::
|
||||
|
||||
y_j = pow(sqrt(-1)*j*2*pi/period, order) * x_j
|
||||
y_0 = 0 if order is not 0.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
Input array.
|
||||
order : int, optional
|
||||
The order of differentiation. Default order is 1. If order is
|
||||
negative, then integration is carried out under the assumption
|
||||
that ``x_0 == 0``.
|
||||
period : float, optional
|
||||
The assumed period of the sequence. Default is ``2*pi``.
|
||||
|
||||
Notes
|
||||
-----
|
||||
If ``sum(x, axis=0) = 0`` then ``diff(diff(x, k), -k) == x`` (within
|
||||
numerical accuracy).
|
||||
|
||||
For odd order and even ``len(x)``, the Nyquist mode is taken zero.
|
||||
|
||||
"""
|
||||
tmp = asarray(x)
|
||||
if order == 0:
|
||||
return tmp
|
||||
if iscomplexobj(tmp):
|
||||
return diff(tmp.real,order,period)+1j*diff(tmp.imag,order,period)
|
||||
if period is not None:
|
||||
c = 2*pi/period
|
||||
else:
|
||||
c = 1.0
|
||||
n = len(x)
|
||||
omega = _cache.get((n,order,c))
|
||||
if omega is None:
|
||||
if len(_cache) > 20:
|
||||
while _cache:
|
||||
_cache.popitem()
|
||||
|
||||
def kernel(k,order=order,c=c):
|
||||
if k:
|
||||
return pow(c*k,order)
|
||||
return 0
|
||||
omega = convolve.init_convolution_kernel(n,kernel,d=order,
|
||||
zero_nyquist=1)
|
||||
_cache[(n,order,c)] = omega
|
||||
overwrite_x = _datacopied(tmp, x)
|
||||
return convolve.convolve(tmp,omega,swap_real_imag=order % 2,
|
||||
overwrite_x=overwrite_x)
|
||||
|
||||
|
||||
del _cache
|
||||
|
||||
|
||||
_cache = {}
|
||||
|
||||
|
||||
def tilbert(x, h, period=None, _cache=_cache):
|
||||
"""
|
||||
Return h-Tilbert transform of a periodic sequence x.
|
||||
|
||||
If x_j and y_j are Fourier coefficients of periodic functions x
|
||||
and y, respectively, then::
|
||||
|
||||
y_j = sqrt(-1)*coth(j*h*2*pi/period) * x_j
|
||||
y_0 = 0
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The input array to transform.
|
||||
h : float
|
||||
Defines the parameter of the Tilbert transform.
|
||||
period : float, optional
|
||||
The assumed period of the sequence. Default period is ``2*pi``.
|
||||
|
||||
Returns
|
||||
-------
|
||||
tilbert : ndarray
|
||||
The result of the transform.
|
||||
|
||||
Notes
|
||||
-----
|
||||
If ``sum(x, axis=0) == 0`` and ``n = len(x)`` is odd then
|
||||
``tilbert(itilbert(x)) == x``.
|
||||
|
||||
If ``2 * pi * h / period`` is approximately 10 or larger, then
|
||||
numerically ``tilbert == hilbert``
|
||||
(theoretically oo-Tilbert == Hilbert).
|
||||
|
||||
For even ``len(x)``, the Nyquist mode of ``x`` is taken zero.
|
||||
|
||||
"""
|
||||
tmp = asarray(x)
|
||||
if iscomplexobj(tmp):
|
||||
return tilbert(tmp.real, h, period) + \
|
||||
1j * tilbert(tmp.imag, h, period)
|
||||
|
||||
if period is not None:
|
||||
h = h * 2 * pi / period
|
||||
|
||||
n = len(x)
|
||||
omega = _cache.get((n, h))
|
||||
if omega is None:
|
||||
if len(_cache) > 20:
|
||||
while _cache:
|
||||
_cache.popitem()
|
||||
|
||||
def kernel(k, h=h):
|
||||
if k:
|
||||
return 1.0/tanh(h*k)
|
||||
|
||||
return 0
|
||||
|
||||
omega = convolve.init_convolution_kernel(n, kernel, d=1)
|
||||
_cache[(n,h)] = omega
|
||||
|
||||
overwrite_x = _datacopied(tmp, x)
|
||||
return convolve.convolve(tmp,omega,swap_real_imag=1,overwrite_x=overwrite_x)
|
||||
|
||||
|
||||
del _cache
|
||||
|
||||
|
||||
_cache = {}
|
||||
|
||||
|
||||
def itilbert(x,h,period=None, _cache=_cache):
|
||||
"""
|
||||
Return inverse h-Tilbert transform of a periodic sequence x.
|
||||
|
||||
If ``x_j`` and ``y_j`` are Fourier coefficients of periodic functions x
|
||||
and y, respectively, then::
|
||||
|
||||
y_j = -sqrt(-1)*tanh(j*h*2*pi/period) * x_j
|
||||
y_0 = 0
|
||||
|
||||
For more details, see `tilbert`.
|
||||
|
||||
"""
|
||||
tmp = asarray(x)
|
||||
if iscomplexobj(tmp):
|
||||
return itilbert(tmp.real,h,period) + \
|
||||
1j*itilbert(tmp.imag,h,period)
|
||||
if period is not None:
|
||||
h = h*2*pi/period
|
||||
n = len(x)
|
||||
omega = _cache.get((n,h))
|
||||
if omega is None:
|
||||
if len(_cache) > 20:
|
||||
while _cache:
|
||||
_cache.popitem()
|
||||
|
||||
def kernel(k,h=h):
|
||||
if k:
|
||||
return -tanh(h*k)
|
||||
return 0
|
||||
omega = convolve.init_convolution_kernel(n,kernel,d=1)
|
||||
_cache[(n,h)] = omega
|
||||
overwrite_x = _datacopied(tmp, x)
|
||||
return convolve.convolve(tmp,omega,swap_real_imag=1,overwrite_x=overwrite_x)
|
||||
|
||||
|
||||
del _cache
|
||||
|
||||
|
||||
_cache = {}
|
||||
|
||||
|
||||
def hilbert(x, _cache=_cache):
|
||||
"""
|
||||
Return Hilbert transform of a periodic sequence x.
|
||||
|
||||
If x_j and y_j are Fourier coefficients of periodic functions x
|
||||
and y, respectively, then::
|
||||
|
||||
y_j = sqrt(-1)*sign(j) * x_j
|
||||
y_0 = 0
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The input array, should be periodic.
|
||||
_cache : dict, optional
|
||||
Dictionary that contains the kernel used to do a convolution with.
|
||||
|
||||
Returns
|
||||
-------
|
||||
y : ndarray
|
||||
The transformed input.
|
||||
|
||||
See Also
|
||||
--------
|
||||
scipy.signal.hilbert : Compute the analytic signal, using the Hilbert
|
||||
transform.
|
||||
|
||||
Notes
|
||||
-----
|
||||
If ``sum(x, axis=0) == 0`` then ``hilbert(ihilbert(x)) == x``.
|
||||
|
||||
For even len(x), the Nyquist mode of x is taken zero.
|
||||
|
||||
The sign of the returned transform does not have a factor -1 that is more
|
||||
often than not found in the definition of the Hilbert transform. Note also
|
||||
that `scipy.signal.hilbert` does have an extra -1 factor compared to this
|
||||
function.
|
||||
|
||||
"""
|
||||
tmp = asarray(x)
|
||||
if iscomplexobj(tmp):
|
||||
return hilbert(tmp.real)+1j*hilbert(tmp.imag)
|
||||
n = len(x)
|
||||
omega = _cache.get(n)
|
||||
if omega is None:
|
||||
if len(_cache) > 20:
|
||||
while _cache:
|
||||
_cache.popitem()
|
||||
|
||||
def kernel(k):
|
||||
if k > 0:
|
||||
return 1.0
|
||||
elif k < 0:
|
||||
return -1.0
|
||||
return 0.0
|
||||
omega = convolve.init_convolution_kernel(n,kernel,d=1)
|
||||
_cache[n] = omega
|
||||
overwrite_x = _datacopied(tmp, x)
|
||||
return convolve.convolve(tmp,omega,swap_real_imag=1,overwrite_x=overwrite_x)
|
||||
|
||||
|
||||
del _cache
|
||||
|
||||
|
||||
def ihilbert(x):
|
||||
"""
|
||||
Return inverse Hilbert transform of a periodic sequence x.
|
||||
|
||||
If ``x_j`` and ``y_j`` are Fourier coefficients of periodic functions x
|
||||
and y, respectively, then::
|
||||
|
||||
y_j = -sqrt(-1)*sign(j) * x_j
|
||||
y_0 = 0
|
||||
|
||||
"""
|
||||
return -hilbert(x)
|
||||
|
||||
|
||||
_cache = {}
|
||||
|
||||
|
||||
def cs_diff(x, a, b, period=None, _cache=_cache):
|
||||
"""
|
||||
Return (a,b)-cosh/sinh pseudo-derivative of a periodic sequence.
|
||||
|
||||
If ``x_j`` and ``y_j`` are Fourier coefficients of periodic functions x
|
||||
and y, respectively, then::
|
||||
|
||||
y_j = -sqrt(-1)*cosh(j*a*2*pi/period)/sinh(j*b*2*pi/period) * x_j
|
||||
y_0 = 0
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The array to take the pseudo-derivative from.
|
||||
a, b : float
|
||||
Defines the parameters of the cosh/sinh pseudo-differential
|
||||
operator.
|
||||
period : float, optional
|
||||
The period of the sequence. Default period is ``2*pi``.
|
||||
|
||||
Returns
|
||||
-------
|
||||
cs_diff : ndarray
|
||||
Pseudo-derivative of periodic sequence `x`.
|
||||
|
||||
Notes
|
||||
-----
|
||||
For even len(`x`), the Nyquist mode of `x` is taken as zero.
|
||||
|
||||
"""
|
||||
tmp = asarray(x)
|
||||
if iscomplexobj(tmp):
|
||||
return cs_diff(tmp.real,a,b,period) + \
|
||||
1j*cs_diff(tmp.imag,a,b,period)
|
||||
if period is not None:
|
||||
a = a*2*pi/period
|
||||
b = b*2*pi/period
|
||||
n = len(x)
|
||||
omega = _cache.get((n,a,b))
|
||||
if omega is None:
|
||||
if len(_cache) > 20:
|
||||
while _cache:
|
||||
_cache.popitem()
|
||||
|
||||
def kernel(k,a=a,b=b):
|
||||
if k:
|
||||
return -cosh(a*k)/sinh(b*k)
|
||||
return 0
|
||||
omega = convolve.init_convolution_kernel(n,kernel,d=1)
|
||||
_cache[(n,a,b)] = omega
|
||||
overwrite_x = _datacopied(tmp, x)
|
||||
return convolve.convolve(tmp,omega,swap_real_imag=1,overwrite_x=overwrite_x)
|
||||
|
||||
|
||||
del _cache
|
||||
|
||||
|
||||
_cache = {}
|
||||
|
||||
|
||||
def sc_diff(x, a, b, period=None, _cache=_cache):
|
||||
"""
|
||||
Return (a,b)-sinh/cosh pseudo-derivative of a periodic sequence x.
|
||||
|
||||
If x_j and y_j are Fourier coefficients of periodic functions x
|
||||
and y, respectively, then::
|
||||
|
||||
y_j = sqrt(-1)*sinh(j*a*2*pi/period)/cosh(j*b*2*pi/period) * x_j
|
||||
y_0 = 0
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
Input array.
|
||||
a,b : float
|
||||
Defines the parameters of the sinh/cosh pseudo-differential
|
||||
operator.
|
||||
period : float, optional
|
||||
The period of the sequence x. Default is 2*pi.
|
||||
|
||||
Notes
|
||||
-----
|
||||
``sc_diff(cs_diff(x,a,b),b,a) == x``
|
||||
For even ``len(x)``, the Nyquist mode of x is taken as zero.
|
||||
|
||||
"""
|
||||
tmp = asarray(x)
|
||||
if iscomplexobj(tmp):
|
||||
return sc_diff(tmp.real,a,b,period) + \
|
||||
1j*sc_diff(tmp.imag,a,b,period)
|
||||
if period is not None:
|
||||
a = a*2*pi/period
|
||||
b = b*2*pi/period
|
||||
n = len(x)
|
||||
omega = _cache.get((n,a,b))
|
||||
if omega is None:
|
||||
if len(_cache) > 20:
|
||||
while _cache:
|
||||
_cache.popitem()
|
||||
|
||||
def kernel(k,a=a,b=b):
|
||||
if k:
|
||||
return sinh(a*k)/cosh(b*k)
|
||||
return 0
|
||||
omega = convolve.init_convolution_kernel(n,kernel,d=1)
|
||||
_cache[(n,a,b)] = omega
|
||||
overwrite_x = _datacopied(tmp, x)
|
||||
return convolve.convolve(tmp,omega,swap_real_imag=1,overwrite_x=overwrite_x)
|
||||
|
||||
|
||||
del _cache
|
||||
|
||||
|
||||
_cache = {}
|
||||
|
||||
|
||||
def ss_diff(x, a, b, period=None, _cache=_cache):
|
||||
"""
|
||||
Return (a,b)-sinh/sinh pseudo-derivative of a periodic sequence x.
|
||||
|
||||
If x_j and y_j are Fourier coefficients of periodic functions x
|
||||
and y, respectively, then::
|
||||
|
||||
y_j = sinh(j*a*2*pi/period)/sinh(j*b*2*pi/period) * x_j
|
||||
y_0 = a/b * x_0
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The array to take the pseudo-derivative from.
|
||||
a,b
|
||||
Defines the parameters of the sinh/sinh pseudo-differential
|
||||
operator.
|
||||
period : float, optional
|
||||
The period of the sequence x. Default is ``2*pi``.
|
||||
|
||||
Notes
|
||||
-----
|
||||
``ss_diff(ss_diff(x,a,b),b,a) == x``
|
||||
|
||||
"""
|
||||
tmp = asarray(x)
|
||||
if iscomplexobj(tmp):
|
||||
return ss_diff(tmp.real,a,b,period) + \
|
||||
1j*ss_diff(tmp.imag,a,b,period)
|
||||
if period is not None:
|
||||
a = a*2*pi/period
|
||||
b = b*2*pi/period
|
||||
n = len(x)
|
||||
omega = _cache.get((n,a,b))
|
||||
if omega is None:
|
||||
if len(_cache) > 20:
|
||||
while _cache:
|
||||
_cache.popitem()
|
||||
|
||||
def kernel(k,a=a,b=b):
|
||||
if k:
|
||||
return sinh(a*k)/sinh(b*k)
|
||||
return float(a)/b
|
||||
omega = convolve.init_convolution_kernel(n,kernel)
|
||||
_cache[(n,a,b)] = omega
|
||||
overwrite_x = _datacopied(tmp, x)
|
||||
return convolve.convolve(tmp,omega,overwrite_x=overwrite_x)
|
||||
|
||||
|
||||
del _cache
|
||||
|
||||
|
||||
_cache = {}
|
||||
|
||||
|
||||
def cc_diff(x, a, b, period=None, _cache=_cache):
|
||||
"""
|
||||
Return (a,b)-cosh/cosh pseudo-derivative of a periodic sequence.
|
||||
|
||||
If x_j and y_j are Fourier coefficients of periodic functions x
|
||||
and y, respectively, then::
|
||||
|
||||
y_j = cosh(j*a*2*pi/period)/cosh(j*b*2*pi/period) * x_j
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The array to take the pseudo-derivative from.
|
||||
a,b : float
|
||||
Defines the parameters of the sinh/sinh pseudo-differential
|
||||
operator.
|
||||
period : float, optional
|
||||
The period of the sequence x. Default is ``2*pi``.
|
||||
|
||||
Returns
|
||||
-------
|
||||
cc_diff : ndarray
|
||||
Pseudo-derivative of periodic sequence `x`.
|
||||
|
||||
Notes
|
||||
-----
|
||||
``cc_diff(cc_diff(x,a,b),b,a) == x``
|
||||
|
||||
"""
|
||||
tmp = asarray(x)
|
||||
if iscomplexobj(tmp):
|
||||
return cc_diff(tmp.real,a,b,period) + \
|
||||
1j*cc_diff(tmp.imag,a,b,period)
|
||||
if period is not None:
|
||||
a = a*2*pi/period
|
||||
b = b*2*pi/period
|
||||
n = len(x)
|
||||
omega = _cache.get((n,a,b))
|
||||
if omega is None:
|
||||
if len(_cache) > 20:
|
||||
while _cache:
|
||||
_cache.popitem()
|
||||
|
||||
def kernel(k,a=a,b=b):
|
||||
return cosh(a*k)/cosh(b*k)
|
||||
omega = convolve.init_convolution_kernel(n,kernel)
|
||||
_cache[(n,a,b)] = omega
|
||||
overwrite_x = _datacopied(tmp, x)
|
||||
return convolve.convolve(tmp,omega,overwrite_x=overwrite_x)
|
||||
|
||||
|
||||
del _cache
|
||||
|
||||
|
||||
_cache = {}
|
||||
|
||||
|
||||
def shift(x, a, period=None, _cache=_cache):
|
||||
"""
|
||||
Shift periodic sequence x by a: y(u) = x(u+a).
|
||||
|
||||
If x_j and y_j are Fourier coefficients of periodic functions x
|
||||
and y, respectively, then::
|
||||
|
||||
y_j = exp(j*a*2*pi/period*sqrt(-1)) * x_f
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The array to take the pseudo-derivative from.
|
||||
a : float
|
||||
Defines the parameters of the sinh/sinh pseudo-differential
|
||||
period : float, optional
|
||||
The period of the sequences x and y. Default period is ``2*pi``.
|
||||
"""
|
||||
tmp = asarray(x)
|
||||
if iscomplexobj(tmp):
|
||||
return shift(tmp.real,a,period)+1j*shift(tmp.imag,a,period)
|
||||
if period is not None:
|
||||
a = a*2*pi/period
|
||||
n = len(x)
|
||||
omega = _cache.get((n,a))
|
||||
if omega is None:
|
||||
if len(_cache) > 20:
|
||||
while _cache:
|
||||
_cache.popitem()
|
||||
|
||||
def kernel_real(k,a=a):
|
||||
return cos(a*k)
|
||||
|
||||
def kernel_imag(k,a=a):
|
||||
return sin(a*k)
|
||||
omega_real = convolve.init_convolution_kernel(n,kernel_real,d=0,
|
||||
zero_nyquist=0)
|
||||
omega_imag = convolve.init_convolution_kernel(n,kernel_imag,d=1,
|
||||
zero_nyquist=0)
|
||||
_cache[(n,a)] = omega_real,omega_imag
|
||||
else:
|
||||
omega_real,omega_imag = omega
|
||||
overwrite_x = _datacopied(tmp, x)
|
||||
return convolve.convolve_z(tmp,omega_real,omega_imag,
|
||||
overwrite_x=overwrite_x)
|
||||
|
||||
|
||||
del _cache
|
||||
@@ -0,0 +1,739 @@
|
||||
"""
|
||||
Real spectrum transforms (DCT, DST, MDCT)
|
||||
"""
|
||||
from __future__ import division, print_function, absolute_import
|
||||
|
||||
|
||||
__all__ = ['dct', 'idct', 'dst', 'idst', 'dctn', 'idctn', 'dstn', 'idstn']
|
||||
|
||||
import numpy as np
|
||||
from scipy.fftpack import _fftpack
|
||||
from scipy.fftpack.basic import _datacopied, _fix_shape, _asfarray
|
||||
from scipy.fftpack.helper import _init_nd_shape_and_axes
|
||||
|
||||
import atexit
|
||||
atexit.register(_fftpack.destroy_ddct1_cache)
|
||||
atexit.register(_fftpack.destroy_ddct2_cache)
|
||||
atexit.register(_fftpack.destroy_ddct4_cache)
|
||||
atexit.register(_fftpack.destroy_dct1_cache)
|
||||
atexit.register(_fftpack.destroy_dct2_cache)
|
||||
atexit.register(_fftpack.destroy_dct4_cache)
|
||||
|
||||
atexit.register(_fftpack.destroy_ddst1_cache)
|
||||
atexit.register(_fftpack.destroy_ddst2_cache)
|
||||
atexit.register(_fftpack.destroy_dst1_cache)
|
||||
atexit.register(_fftpack.destroy_dst2_cache)
|
||||
|
||||
|
||||
def dctn(x, type=2, shape=None, axes=None, norm=None, overwrite_x=False):
|
||||
"""
|
||||
Return multidimensional Discrete Cosine Transform along the specified axes.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The input array.
|
||||
type : {1, 2, 3, 4}, optional
|
||||
Type of the DCT (see Notes). Default type is 2.
|
||||
shape : int or array_like of ints or None, optional
|
||||
The shape of the result. If both `shape` and `axes` (see below) are
|
||||
None, `shape` is ``x.shape``; if `shape` is None but `axes` is
|
||||
not None, then `shape` is ``scipy.take(x.shape, axes, axis=0)``.
|
||||
If ``shape[i] > x.shape[i]``, the i-th dimension is padded with zeros.
|
||||
If ``shape[i] < x.shape[i]``, the i-th dimension is truncated to
|
||||
length ``shape[i]``.
|
||||
If any element of `shape` is -1, the size of the corresponding
|
||||
dimension of `x` is used.
|
||||
axes : int or array_like of ints or None, optional
|
||||
Axes along which the DCT is computed.
|
||||
The default is over all axes.
|
||||
norm : {None, 'ortho'}, optional
|
||||
Normalization mode (see Notes). Default is None.
|
||||
overwrite_x : bool, optional
|
||||
If True, the contents of `x` can be destroyed; the default is False.
|
||||
|
||||
Returns
|
||||
-------
|
||||
y : ndarray of real
|
||||
The transformed input array.
|
||||
|
||||
See Also
|
||||
--------
|
||||
idctn : Inverse multidimensional DCT
|
||||
|
||||
Notes
|
||||
-----
|
||||
For full details of the DCT types and normalization modes, as well as
|
||||
references, see `dct`.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from scipy.fftpack import dctn, idctn
|
||||
>>> y = np.random.randn(16, 16)
|
||||
>>> np.allclose(y, idctn(dctn(y, norm='ortho'), norm='ortho'))
|
||||
True
|
||||
|
||||
"""
|
||||
x = np.asanyarray(x)
|
||||
shape, axes = _init_nd_shape_and_axes(x, shape, axes)
|
||||
for n, ax in zip(shape, axes):
|
||||
x = dct(x, type=type, n=n, axis=ax, norm=norm, overwrite_x=overwrite_x)
|
||||
return x
|
||||
|
||||
|
||||
def idctn(x, type=2, shape=None, axes=None, norm=None, overwrite_x=False):
|
||||
"""
|
||||
Return multidimensional Discrete Cosine Transform along the specified axes.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The input array.
|
||||
type : {1, 2, 3, 4}, optional
|
||||
Type of the DCT (see Notes). Default type is 2.
|
||||
shape : int or array_like of ints or None, optional
|
||||
The shape of the result. If both `shape` and `axes` (see below) are
|
||||
None, `shape` is ``x.shape``; if `shape` is None but `axes` is
|
||||
not None, then `shape` is ``scipy.take(x.shape, axes, axis=0)``.
|
||||
If ``shape[i] > x.shape[i]``, the i-th dimension is padded with zeros.
|
||||
If ``shape[i] < x.shape[i]``, the i-th dimension is truncated to
|
||||
length ``shape[i]``.
|
||||
If any element of `shape` is -1, the size of the corresponding
|
||||
dimension of `x` is used.
|
||||
axes : int or array_like of ints or None, optional
|
||||
Axes along which the IDCT is computed.
|
||||
The default is over all axes.
|
||||
norm : {None, 'ortho'}, optional
|
||||
Normalization mode (see Notes). Default is None.
|
||||
overwrite_x : bool, optional
|
||||
If True, the contents of `x` can be destroyed; the default is False.
|
||||
|
||||
Returns
|
||||
-------
|
||||
y : ndarray of real
|
||||
The transformed input array.
|
||||
|
||||
See Also
|
||||
--------
|
||||
dctn : multidimensional DCT
|
||||
|
||||
Notes
|
||||
-----
|
||||
For full details of the IDCT types and normalization modes, as well as
|
||||
references, see `idct`.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from scipy.fftpack import dctn, idctn
|
||||
>>> y = np.random.randn(16, 16)
|
||||
>>> np.allclose(y, idctn(dctn(y, norm='ortho'), norm='ortho'))
|
||||
True
|
||||
|
||||
"""
|
||||
x = np.asanyarray(x)
|
||||
shape, axes = _init_nd_shape_and_axes(x, shape, axes)
|
||||
for n, ax in zip(shape, axes):
|
||||
x = idct(x, type=type, n=n, axis=ax, norm=norm,
|
||||
overwrite_x=overwrite_x)
|
||||
return x
|
||||
|
||||
|
||||
def dstn(x, type=2, shape=None, axes=None, norm=None, overwrite_x=False):
|
||||
"""
|
||||
Return multidimensional Discrete Sine Transform along the specified axes.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The input array.
|
||||
type : {1, 2, 3, 4}, optional
|
||||
Type of the DCT (see Notes). Default type is 2.
|
||||
shape : int or array_like of ints or None, optional
|
||||
The shape of the result. If both `shape` and `axes` (see below) are
|
||||
None, `shape` is ``x.shape``; if `shape` is None but `axes` is
|
||||
not None, then `shape` is ``scipy.take(x.shape, axes, axis=0)``.
|
||||
If ``shape[i] > x.shape[i]``, the i-th dimension is padded with zeros.
|
||||
If ``shape[i] < x.shape[i]``, the i-th dimension is truncated to
|
||||
length ``shape[i]``.
|
||||
If any element of `shape` is -1, the size of the corresponding
|
||||
dimension of `x` is used.
|
||||
axes : int or array_like of ints or None, optional
|
||||
Axes along which the DCT is computed.
|
||||
The default is over all axes.
|
||||
norm : {None, 'ortho'}, optional
|
||||
Normalization mode (see Notes). Default is None.
|
||||
overwrite_x : bool, optional
|
||||
If True, the contents of `x` can be destroyed; the default is False.
|
||||
|
||||
Returns
|
||||
-------
|
||||
y : ndarray of real
|
||||
The transformed input array.
|
||||
|
||||
See Also
|
||||
--------
|
||||
idstn : Inverse multidimensional DST
|
||||
|
||||
Notes
|
||||
-----
|
||||
For full details of the DST types and normalization modes, as well as
|
||||
references, see `dst`.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from scipy.fftpack import dstn, idstn
|
||||
>>> y = np.random.randn(16, 16)
|
||||
>>> np.allclose(y, idstn(dstn(y, norm='ortho'), norm='ortho'))
|
||||
True
|
||||
|
||||
"""
|
||||
x = np.asanyarray(x)
|
||||
shape, axes = _init_nd_shape_and_axes(x, shape, axes)
|
||||
for n, ax in zip(shape, axes):
|
||||
x = dst(x, type=type, n=n, axis=ax, norm=norm, overwrite_x=overwrite_x)
|
||||
return x
|
||||
|
||||
|
||||
def idstn(x, type=2, shape=None, axes=None, norm=None, overwrite_x=False):
|
||||
"""
|
||||
Return multidimensional Discrete Sine Transform along the specified axes.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The input array.
|
||||
type : {1, 2, 3, 4}, optional
|
||||
Type of the DCT (see Notes). Default type is 2.
|
||||
shape : int or array_like of ints or None, optional
|
||||
The shape of the result. If both `shape` and `axes` (see below) are
|
||||
None, `shape` is ``x.shape``; if `shape` is None but `axes` is
|
||||
not None, then `shape` is ``scipy.take(x.shape, axes, axis=0)``.
|
||||
If ``shape[i] > x.shape[i]``, the i-th dimension is padded with zeros.
|
||||
If ``shape[i] < x.shape[i]``, the i-th dimension is truncated to
|
||||
length ``shape[i]``.
|
||||
If any element of `shape` is -1, the size of the corresponding
|
||||
dimension of `x` is used.
|
||||
axes : int or array_like of ints or None, optional
|
||||
Axes along which the IDCT is computed.
|
||||
The default is over all axes.
|
||||
norm : {None, 'ortho'}, optional
|
||||
Normalization mode (see Notes). Default is None.
|
||||
overwrite_x : bool, optional
|
||||
If True, the contents of `x` can be destroyed; the default is False.
|
||||
|
||||
Returns
|
||||
-------
|
||||
y : ndarray of real
|
||||
The transformed input array.
|
||||
|
||||
See Also
|
||||
--------
|
||||
dctn : multidimensional DST
|
||||
|
||||
Notes
|
||||
-----
|
||||
For full details of the IDST types and normalization modes, as well as
|
||||
references, see `idst`.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from scipy.fftpack import dstn, idstn
|
||||
>>> y = np.random.randn(16, 16)
|
||||
>>> np.allclose(y, idstn(dstn(y, norm='ortho'), norm='ortho'))
|
||||
True
|
||||
|
||||
"""
|
||||
x = np.asanyarray(x)
|
||||
shape, axes = _init_nd_shape_and_axes(x, shape, axes)
|
||||
for n, ax in zip(shape, axes):
|
||||
x = idst(x, type=type, n=n, axis=ax, norm=norm,
|
||||
overwrite_x=overwrite_x)
|
||||
return x
|
||||
|
||||
|
||||
def dct(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False):
|
||||
"""
|
||||
Return the Discrete Cosine Transform of arbitrary type sequence x.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The input array.
|
||||
type : {1, 2, 3, 4}, optional
|
||||
Type of the DCT (see Notes). Default type is 2.
|
||||
n : int, optional
|
||||
Length of the transform. If ``n < x.shape[axis]``, `x` is
|
||||
truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
|
||||
default results in ``n = x.shape[axis]``.
|
||||
axis : int, optional
|
||||
Axis along which the dct is computed; the default is over the
|
||||
last axis (i.e., ``axis=-1``).
|
||||
norm : {None, 'ortho'}, optional
|
||||
Normalization mode (see Notes). Default is None.
|
||||
overwrite_x : bool, optional
|
||||
If True, the contents of `x` can be destroyed; the default is False.
|
||||
|
||||
Returns
|
||||
-------
|
||||
y : ndarray of real
|
||||
The transformed input array.
|
||||
|
||||
See Also
|
||||
--------
|
||||
idct : Inverse DCT
|
||||
|
||||
Notes
|
||||
-----
|
||||
For a single dimension array ``x``, ``dct(x, norm='ortho')`` is equal to
|
||||
MATLAB ``dct(x)``.
|
||||
|
||||
There are theoretically 8 types of the DCT, only the first 4 types are
|
||||
implemented in scipy. 'The' DCT generally refers to DCT type 2, and 'the'
|
||||
Inverse DCT generally refers to DCT type 3.
|
||||
|
||||
**Type I**
|
||||
|
||||
There are several definitions of the DCT-I; we use the following
|
||||
(for ``norm=None``)::
|
||||
|
||||
N-2
|
||||
y[k] = x[0] + (-1)**k x[N-1] + 2 * sum x[n]*cos(pi*k*n/(N-1))
|
||||
n=1
|
||||
|
||||
If ``norm='ortho'``, ``x[0]`` and ``x[N-1]`` are multiplied by
|
||||
a scaling factor of ``sqrt(2)``, and ``y[k]`` is multiplied by a
|
||||
scaling factor `f`::
|
||||
|
||||
f = 0.5*sqrt(1/(N-1)) if k = 0 or N-1,
|
||||
f = 0.5*sqrt(2/(N-1)) otherwise.
|
||||
|
||||
.. versionadded:: 1.2.0
|
||||
Orthonormalization in DCT-I.
|
||||
|
||||
.. note::
|
||||
The DCT-I is only supported for input size > 1.
|
||||
|
||||
**Type II**
|
||||
|
||||
There are several definitions of the DCT-II; we use the following
|
||||
(for ``norm=None``)::
|
||||
|
||||
|
||||
N-1
|
||||
y[k] = 2* sum x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N.
|
||||
n=0
|
||||
|
||||
If ``norm='ortho'``, ``y[k]`` is multiplied by a scaling factor `f`::
|
||||
|
||||
f = sqrt(1/(4*N)) if k = 0,
|
||||
f = sqrt(1/(2*N)) otherwise.
|
||||
|
||||
Which makes the corresponding matrix of coefficients orthonormal
|
||||
(``OO' = Id``).
|
||||
|
||||
**Type III**
|
||||
|
||||
There are several definitions, we use the following
|
||||
(for ``norm=None``)::
|
||||
|
||||
N-1
|
||||
y[k] = x[0] + 2 * sum x[n]*cos(pi*(k+0.5)*n/N), 0 <= k < N.
|
||||
n=1
|
||||
|
||||
or, for ``norm='ortho'`` and 0 <= k < N::
|
||||
|
||||
N-1
|
||||
y[k] = x[0] / sqrt(N) + sqrt(2/N) * sum x[n]*cos(pi*(k+0.5)*n/N)
|
||||
n=1
|
||||
|
||||
The (unnormalized) DCT-III is the inverse of the (unnormalized) DCT-II, up
|
||||
to a factor `2N`. The orthonormalized DCT-III is exactly the inverse of
|
||||
the orthonormalized DCT-II.
|
||||
|
||||
**Type IV**
|
||||
|
||||
There are several definitions of the DCT-IV; we use the following
|
||||
(for ``norm=None``)::
|
||||
|
||||
|
||||
N-1
|
||||
y[k] = 2* sum x[n]*cos(pi*(2k+1)*(2n+1)/(4*N)), 0 <= k < N.
|
||||
n=0
|
||||
|
||||
If ``norm='ortho'``, ``y[k]`` is multiplied by a scaling factor `f`::
|
||||
|
||||
f = 0.5*sqrt(2/N)
|
||||
|
||||
.. versionadded:: 1.2.0
|
||||
Support for DCT-IV.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] 'A Fast Cosine Transform in One and Two Dimensions', by J.
|
||||
Makhoul, `IEEE Transactions on acoustics, speech and signal
|
||||
processing` vol. 28(1), pp. 27-34,
|
||||
:doi:`10.1109/TASSP.1980.1163351` (1980).
|
||||
.. [2] Wikipedia, "Discrete cosine transform",
|
||||
https://en.wikipedia.org/wiki/Discrete_cosine_transform
|
||||
|
||||
Examples
|
||||
--------
|
||||
The Type 1 DCT is equivalent to the FFT (though faster) for real,
|
||||
even-symmetrical inputs. The output is also real and even-symmetrical.
|
||||
Half of the FFT input is used to generate half of the FFT output:
|
||||
|
||||
>>> from scipy.fftpack import fft, dct
|
||||
>>> fft(np.array([4., 3., 5., 10., 5., 3.])).real
|
||||
array([ 30., -8., 6., -2., 6., -8.])
|
||||
>>> dct(np.array([4., 3., 5., 10.]), 1)
|
||||
array([ 30., -8., 6., -2.])
|
||||
|
||||
"""
|
||||
return _dct(x, type, n, axis, normalize=norm, overwrite_x=overwrite_x)
|
||||
|
||||
|
||||
def idct(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False):
|
||||
"""
|
||||
Return the Inverse Discrete Cosine Transform of an arbitrary type sequence.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The input array.
|
||||
type : {1, 2, 3, 4}, optional
|
||||
Type of the DCT (see Notes). Default type is 2.
|
||||
n : int, optional
|
||||
Length of the transform. If ``n < x.shape[axis]``, `x` is
|
||||
truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
|
||||
default results in ``n = x.shape[axis]``.
|
||||
axis : int, optional
|
||||
Axis along which the idct is computed; the default is over the
|
||||
last axis (i.e., ``axis=-1``).
|
||||
norm : {None, 'ortho'}, optional
|
||||
Normalization mode (see Notes). Default is None.
|
||||
overwrite_x : bool, optional
|
||||
If True, the contents of `x` can be destroyed; the default is False.
|
||||
|
||||
Returns
|
||||
-------
|
||||
idct : ndarray of real
|
||||
The transformed input array.
|
||||
|
||||
See Also
|
||||
--------
|
||||
dct : Forward DCT
|
||||
|
||||
Notes
|
||||
-----
|
||||
For a single dimension array `x`, ``idct(x, norm='ortho')`` is equal to
|
||||
MATLAB ``idct(x)``.
|
||||
|
||||
'The' IDCT is the IDCT of type 2, which is the same as DCT of type 3.
|
||||
|
||||
IDCT of type 1 is the DCT of type 1, IDCT of type 2 is the DCT of type
|
||||
3, and IDCT of type 3 is the DCT of type 2. IDCT of type 4 is the DCT
|
||||
of type 4. For the definition of these types, see `dct`.
|
||||
|
||||
Examples
|
||||
--------
|
||||
The Type 1 DCT is equivalent to the DFT for real, even-symmetrical
|
||||
inputs. The output is also real and even-symmetrical. Half of the IFFT
|
||||
input is used to generate half of the IFFT output:
|
||||
|
||||
>>> from scipy.fftpack import ifft, idct
|
||||
>>> ifft(np.array([ 30., -8., 6., -2., 6., -8.])).real
|
||||
array([ 4., 3., 5., 10., 5., 3.])
|
||||
>>> idct(np.array([ 30., -8., 6., -2.]), 1) / 6
|
||||
array([ 4., 3., 5., 10.])
|
||||
|
||||
"""
|
||||
# Inverse/forward type table
|
||||
_TP = {1:1, 2:3, 3:2, 4:4}
|
||||
return _dct(x, _TP[type], n, axis, normalize=norm, overwrite_x=overwrite_x)
|
||||
|
||||
|
||||
def _get_dct_fun(type, dtype):
|
||||
try:
|
||||
name = {'float64':'ddct%d', 'float32':'dct%d'}[dtype.name]
|
||||
except KeyError:
|
||||
raise ValueError("dtype %s not supported" % dtype)
|
||||
try:
|
||||
f = getattr(_fftpack, name % type)
|
||||
except AttributeError as e:
|
||||
raise ValueError(str(e) + ". Type %d not understood" % type)
|
||||
return f
|
||||
|
||||
|
||||
def _get_norm_mode(normalize):
|
||||
try:
|
||||
nm = {None:0, 'ortho':1}[normalize]
|
||||
except KeyError:
|
||||
raise ValueError("Unknown normalize mode %s" % normalize)
|
||||
return nm
|
||||
|
||||
|
||||
def __fix_shape(x, n, axis, dct_or_dst):
|
||||
tmp = _asfarray(x)
|
||||
copy_made = _datacopied(tmp, x)
|
||||
if n is None:
|
||||
n = tmp.shape[axis]
|
||||
elif n != tmp.shape[axis]:
|
||||
tmp, copy_made2 = _fix_shape(tmp, n, axis)
|
||||
copy_made = copy_made or copy_made2
|
||||
if n < 1:
|
||||
raise ValueError("Invalid number of %s data points "
|
||||
"(%d) specified." % (dct_or_dst, n))
|
||||
return tmp, n, copy_made
|
||||
|
||||
|
||||
def _raw_dct(x0, type, n, axis, nm, overwrite_x):
|
||||
f = _get_dct_fun(type, x0.dtype)
|
||||
return _eval_fun(f, x0, n, axis, nm, overwrite_x)
|
||||
|
||||
|
||||
def _raw_dst(x0, type, n, axis, nm, overwrite_x):
|
||||
f = _get_dst_fun(type, x0.dtype)
|
||||
return _eval_fun(f, x0, n, axis, nm, overwrite_x)
|
||||
|
||||
|
||||
def _eval_fun(f, tmp, n, axis, nm, overwrite_x):
|
||||
if axis == -1 or axis == len(tmp.shape) - 1:
|
||||
return f(tmp, n, nm, overwrite_x)
|
||||
|
||||
tmp = np.swapaxes(tmp, axis, -1)
|
||||
tmp = f(tmp, n, nm, overwrite_x)
|
||||
return np.swapaxes(tmp, axis, -1)
|
||||
|
||||
|
||||
def _dct(x, type, n=None, axis=-1, overwrite_x=False, normalize=None):
|
||||
"""
|
||||
Return Discrete Cosine Transform of arbitrary type sequence x.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
input array.
|
||||
n : int, optional
|
||||
Length of the transform. If ``n < x.shape[axis]``, `x` is
|
||||
truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
|
||||
default results in ``n = x.shape[axis]``.
|
||||
axis : int, optional
|
||||
Axis along which the dct is computed; the default is over the
|
||||
last axis (i.e., ``axis=-1``).
|
||||
overwrite_x : bool, optional
|
||||
If True, the contents of `x` can be destroyed; the default is False.
|
||||
|
||||
Returns
|
||||
-------
|
||||
z : ndarray
|
||||
|
||||
"""
|
||||
x0, n, copy_made = __fix_shape(x, n, axis, 'DCT')
|
||||
if type == 1 and n < 2:
|
||||
raise ValueError("DCT-I is not defined for size < 2")
|
||||
overwrite_x = overwrite_x or copy_made
|
||||
nm = _get_norm_mode(normalize)
|
||||
if np.iscomplexobj(x0):
|
||||
return (_raw_dct(x0.real, type, n, axis, nm, overwrite_x) + 1j *
|
||||
_raw_dct(x0.imag, type, n, axis, nm, overwrite_x))
|
||||
else:
|
||||
return _raw_dct(x0, type, n, axis, nm, overwrite_x)
|
||||
|
||||
|
||||
def dst(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False):
|
||||
"""
|
||||
Return the Discrete Sine Transform of arbitrary type sequence x.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The input array.
|
||||
type : {1, 2, 3, 4}, optional
|
||||
Type of the DST (see Notes). Default type is 2.
|
||||
n : int, optional
|
||||
Length of the transform. If ``n < x.shape[axis]``, `x` is
|
||||
truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
|
||||
default results in ``n = x.shape[axis]``.
|
||||
axis : int, optional
|
||||
Axis along which the dst is computed; the default is over the
|
||||
last axis (i.e., ``axis=-1``).
|
||||
norm : {None, 'ortho'}, optional
|
||||
Normalization mode (see Notes). Default is None.
|
||||
overwrite_x : bool, optional
|
||||
If True, the contents of `x` can be destroyed; the default is False.
|
||||
|
||||
Returns
|
||||
-------
|
||||
dst : ndarray of reals
|
||||
The transformed input array.
|
||||
|
||||
See Also
|
||||
--------
|
||||
idst : Inverse DST
|
||||
|
||||
Notes
|
||||
-----
|
||||
For a single dimension array ``x``.
|
||||
|
||||
There are theoretically 8 types of the DST for different combinations of
|
||||
even/odd boundary conditions and boundary off sets [1]_, only the first
|
||||
3 types are implemented in scipy.
|
||||
|
||||
**Type I**
|
||||
|
||||
There are several definitions of the DST-I; we use the following
|
||||
for ``norm=None``. DST-I assumes the input is odd around n=-1 and n=N. ::
|
||||
|
||||
N-1
|
||||
y[k] = 2 * sum x[n]*sin(pi*(k+1)*(n+1)/(N+1))
|
||||
n=0
|
||||
|
||||
Note that the DST-I is only supported for input size > 1
|
||||
The (unnormalized) DST-I is its own inverse, up to a factor `2(N+1)`.
|
||||
The orthonormalized DST-I is exactly its own inverse.
|
||||
|
||||
**Type II**
|
||||
|
||||
There are several definitions of the DST-II; we use the following
|
||||
for ``norm=None``. DST-II assumes the input is odd around n=-1/2 and
|
||||
n=N-1/2; the output is odd around k=-1 and even around k=N-1 ::
|
||||
|
||||
N-1
|
||||
y[k] = 2* sum x[n]*sin(pi*(k+1)*(n+0.5)/N), 0 <= k < N.
|
||||
n=0
|
||||
|
||||
if ``norm='ortho'``, ``y[k]`` is multiplied by a scaling factor `f` ::
|
||||
|
||||
f = sqrt(1/(4*N)) if k == 0
|
||||
f = sqrt(1/(2*N)) otherwise.
|
||||
|
||||
**Type III**
|
||||
|
||||
There are several definitions of the DST-III, we use the following
|
||||
(for ``norm=None``). DST-III assumes the input is odd around n=-1
|
||||
and even around n=N-1 ::
|
||||
|
||||
N-2
|
||||
y[k] = x[N-1]*(-1)**k + 2* sum x[n]*sin(pi*(k+0.5)*(n+1)/N), 0 <= k < N.
|
||||
n=0
|
||||
|
||||
The (unnormalized) DST-III is the inverse of the (unnormalized) DST-II, up
|
||||
to a factor `2N`. The orthonormalized DST-III is exactly the inverse of
|
||||
the orthonormalized DST-II.
|
||||
|
||||
.. versionadded:: 0.11.0
|
||||
|
||||
**Type IV**
|
||||
|
||||
There are several definitions of the DST-IV, we use the following
|
||||
(for ``norm=None``). DST-IV assumes the input is odd around n=-0.5
|
||||
and even around n=N-0.5 ::
|
||||
|
||||
N-1
|
||||
y[k] = 2* sum x[n]*sin(pi*(k+0.5)*(n+0.5)/N), 0 <= k < N.
|
||||
n=0
|
||||
|
||||
The (unnormalized) DST-IV is its own inverse, up
|
||||
to a factor `2N`. The orthonormalized DST-IV is exactly its own inverse.
|
||||
|
||||
.. versionadded:: 1.2.0
|
||||
Support for DST-IV.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] Wikipedia, "Discrete sine transform",
|
||||
https://en.wikipedia.org/wiki/Discrete_sine_transform
|
||||
|
||||
"""
|
||||
return _dst(x, type, n, axis, normalize=norm, overwrite_x=overwrite_x)
|
||||
|
||||
|
||||
def idst(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False):
|
||||
"""
|
||||
Return the Inverse Discrete Sine Transform of an arbitrary type sequence.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The input array.
|
||||
type : {1, 2, 3, 4}, optional
|
||||
Type of the DST (see Notes). Default type is 2.
|
||||
n : int, optional
|
||||
Length of the transform. If ``n < x.shape[axis]``, `x` is
|
||||
truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
|
||||
default results in ``n = x.shape[axis]``.
|
||||
axis : int, optional
|
||||
Axis along which the idst is computed; the default is over the
|
||||
last axis (i.e., ``axis=-1``).
|
||||
norm : {None, 'ortho'}, optional
|
||||
Normalization mode (see Notes). Default is None.
|
||||
overwrite_x : bool, optional
|
||||
If True, the contents of `x` can be destroyed; the default is False.
|
||||
|
||||
Returns
|
||||
-------
|
||||
idst : ndarray of real
|
||||
The transformed input array.
|
||||
|
||||
See Also
|
||||
--------
|
||||
dst : Forward DST
|
||||
|
||||
Notes
|
||||
-----
|
||||
'The' IDST is the IDST of type 2, which is the same as DST of type 3.
|
||||
|
||||
IDST of type 1 is the DST of type 1, IDST of type 2 is the DST of type
|
||||
3, and IDST of type 3 is the DST of type 2. For the definition of these
|
||||
types, see `dst`.
|
||||
|
||||
.. versionadded:: 0.11.0
|
||||
|
||||
"""
|
||||
# Inverse/forward type table
|
||||
_TP = {1:1, 2:3, 3:2, 4:4}
|
||||
return _dst(x, _TP[type], n, axis, normalize=norm, overwrite_x=overwrite_x)
|
||||
|
||||
|
||||
def _get_dst_fun(type, dtype):
|
||||
try:
|
||||
name = {'float64':'ddst%d', 'float32':'dst%d'}[dtype.name]
|
||||
except KeyError:
|
||||
raise ValueError("dtype %s not supported" % dtype)
|
||||
try:
|
||||
f = getattr(_fftpack, name % type)
|
||||
except AttributeError as e:
|
||||
raise ValueError(str(e) + ". Type %d not understood" % type)
|
||||
return f
|
||||
|
||||
|
||||
def _dst(x, type, n=None, axis=-1, overwrite_x=False, normalize=None):
|
||||
"""
|
||||
Return Discrete Sine Transform of arbitrary type sequence x.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
input array.
|
||||
n : int, optional
|
||||
Length of the transform.
|
||||
axis : int, optional
|
||||
Axis along which the dst is computed. (default=-1)
|
||||
overwrite_x : bool, optional
|
||||
If True the contents of x can be destroyed. (default=False)
|
||||
|
||||
Returns
|
||||
-------
|
||||
z : real ndarray
|
||||
|
||||
"""
|
||||
x0, n, copy_made = __fix_shape(x, n, axis, 'DST')
|
||||
if type == 1 and n < 2:
|
||||
raise ValueError("DST-I is not defined for size < 2")
|
||||
overwrite_x = overwrite_x or copy_made
|
||||
nm = _get_norm_mode(normalize)
|
||||
if np.iscomplexobj(x0):
|
||||
return (_raw_dst(x0.real, type, n, axis, nm, overwrite_x) + 1j *
|
||||
_raw_dst(x0.imag, type, n, axis, nm, overwrite_x))
|
||||
else:
|
||||
return _raw_dst(x0, type, n, axis, nm, overwrite_x)
|
||||
@@ -0,0 +1,40 @@
|
||||
# Created by Pearu Peterson, August 2002
|
||||
from __future__ import division, print_function, absolute_import
|
||||
|
||||
|
||||
from os.path import join
|
||||
|
||||
|
||||
def configuration(parent_package='',top_path=None):
|
||||
from numpy.distutils.misc_util import Configuration
|
||||
|
||||
config = Configuration('fftpack',parent_package, top_path)
|
||||
|
||||
config.add_data_dir('tests')
|
||||
|
||||
dfftpack_src = [join('src/dfftpack','*.f')]
|
||||
config.add_library('dfftpack', sources=dfftpack_src)
|
||||
|
||||
fftpack_src = [join('src/fftpack','*.f')]
|
||||
config.add_library('fftpack', sources=fftpack_src)
|
||||
|
||||
sources = ['fftpack.pyf','src/zfft.c','src/drfft.c','src/zrfft.c',
|
||||
'src/zfftnd.c', 'src/dct.c.src', 'src/dst.c.src']
|
||||
|
||||
config.add_extension('_fftpack',
|
||||
sources=sources,
|
||||
libraries=['dfftpack', 'fftpack'],
|
||||
include_dirs=['src'],
|
||||
depends=(dfftpack_src + fftpack_src))
|
||||
|
||||
config.add_extension('convolve',
|
||||
sources=['convolve.pyf','src/convolve.c'],
|
||||
libraries=['dfftpack'],
|
||||
depends=dfftpack_src,
|
||||
)
|
||||
return config
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
from numpy.distutils.core import setup
|
||||
setup(**configuration(top_path='').todict())
|
||||
@@ -0,0 +1,13 @@
|
||||
CC = gcc
|
||||
LD = gcc
|
||||
|
||||
fftw_single: fftw_dct.c
|
||||
$(CC) -W -Wall -DDCT_TEST_USE_SINGLE $< -o $@ -lfftw3f
|
||||
|
||||
fftw_double: fftw_dct.c
|
||||
$(CC) -W -Wall $< -o $@ -lfftw3
|
||||
|
||||
clean:
|
||||
rm -f fftw_single
|
||||
rm -f fftw_double
|
||||
rm -f *.o
|
||||
BIN
Binary file not shown.
BIN
Binary file not shown.
BIN
Binary file not shown.
BIN
Binary file not shown.
BIN
Binary file not shown.
BIN
Binary file not shown.
BIN
Binary file not shown.
BIN
Binary file not shown.
@@ -0,0 +1,138 @@
|
||||
#include <stdlib.h>
|
||||
#include <stdio.h>
|
||||
|
||||
#include <fftw3.h>
|
||||
|
||||
#ifdef DCT_TEST_USE_SINGLE
|
||||
typedef float float_prec;
|
||||
#define PF "%.7f"
|
||||
#define FFTW_PLAN fftwf_plan
|
||||
#define FFTW_MALLOC fftwf_malloc
|
||||
#define FFTW_FREE fftwf_free
|
||||
#define FFTW_PLAN_CREATE fftwf_plan_r2r_1d
|
||||
#define FFTW_EXECUTE fftwf_execute
|
||||
#define FFTW_DESTROY_PLAN fftwf_destroy_plan
|
||||
#define FFTW_CLEANUP fftwf_cleanup
|
||||
#else
|
||||
typedef double float_prec;
|
||||
#define PF "%.18f"
|
||||
#define FFTW_PLAN fftw_plan
|
||||
#define FFTW_MALLOC fftw_malloc
|
||||
#define FFTW_FREE fftw_free
|
||||
#define FFTW_PLAN_CREATE fftw_plan_r2r_1d
|
||||
#define FFTW_EXECUTE fftw_execute
|
||||
#define FFTW_DESTROY_PLAN fftw_destroy_plan
|
||||
#define FFTW_CLEANUP fftw_cleanup
|
||||
#endif
|
||||
|
||||
|
||||
enum type {
|
||||
DCT_I = 1,
|
||||
DCT_II = 2,
|
||||
DCT_III = 3,
|
||||
DCT_IV = 4,
|
||||
DST_I = 5,
|
||||
DST_II = 6,
|
||||
DST_III = 7,
|
||||
DST_IV = 8,
|
||||
};
|
||||
|
||||
int gen(int type, int sz)
|
||||
{
|
||||
float_prec *a, *b;
|
||||
FFTW_PLAN p;
|
||||
int i, tp;
|
||||
|
||||
a = FFTW_MALLOC(sizeof(*a) * sz);
|
||||
if (a == NULL) {
|
||||
fprintf(stderr, "failure\n");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
b = FFTW_MALLOC(sizeof(*b) * sz);
|
||||
if (b == NULL) {
|
||||
fprintf(stderr, "failure\n");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
switch(type) {
|
||||
case DCT_I:
|
||||
tp = FFTW_REDFT00;
|
||||
break;
|
||||
case DCT_II:
|
||||
tp = FFTW_REDFT10;
|
||||
break;
|
||||
case DCT_III:
|
||||
tp = FFTW_REDFT01;
|
||||
break;
|
||||
case DCT_IV:
|
||||
tp = FFTW_REDFT11;
|
||||
break;
|
||||
case DST_I:
|
||||
tp = FFTW_RODFT00;
|
||||
break;
|
||||
case DST_II:
|
||||
tp = FFTW_RODFT10;
|
||||
break;
|
||||
case DST_III:
|
||||
tp = FFTW_RODFT01;
|
||||
break;
|
||||
case DST_IV:
|
||||
tp = FFTW_RODFT11;
|
||||
break;
|
||||
default:
|
||||
fprintf(stderr, "unknown type\n");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
switch(type) {
|
||||
case DCT_I:
|
||||
case DCT_II:
|
||||
case DCT_III:
|
||||
case DCT_IV:
|
||||
for(i=0; i < sz; ++i) {
|
||||
a[i] = i;
|
||||
}
|
||||
break;
|
||||
case DST_I:
|
||||
case DST_II:
|
||||
case DST_III:
|
||||
case DST_IV:
|
||||
/* TODO: what should we do for dst's?*/
|
||||
for(i=0; i < sz; ++i) {
|
||||
a[i] = i;
|
||||
}
|
||||
break;
|
||||
default:
|
||||
fprintf(stderr, "unknown type\n");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
p = FFTW_PLAN_CREATE(sz, a, b, tp, FFTW_ESTIMATE);
|
||||
FFTW_EXECUTE(p);
|
||||
FFTW_DESTROY_PLAN(p);
|
||||
|
||||
for(i=0; i < sz; ++i) {
|
||||
printf(PF"\n", b[i]);
|
||||
}
|
||||
FFTW_FREE(b);
|
||||
FFTW_FREE(a);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
int main(int argc, char* argv[])
|
||||
{
|
||||
int n, tp;
|
||||
|
||||
if (argc < 3) {
|
||||
fprintf(stderr, "missing argument: program type n\n");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
tp = atoi(argv[1]);
|
||||
n = atoi(argv[2]);
|
||||
|
||||
gen(tp, n);
|
||||
FFTW_CLEANUP();
|
||||
|
||||
return 0;
|
||||
}
|
||||
Binary file not shown.
Binary file not shown.
@@ -0,0 +1,59 @@
|
||||
from __future__ import division, print_function, absolute_import
|
||||
|
||||
from subprocess import Popen, PIPE, STDOUT
|
||||
|
||||
import numpy as np
|
||||
|
||||
SZ = [2, 3, 4, 8, 12, 15, 16, 17, 32, 64, 128, 256, 512, 1024]
|
||||
|
||||
|
||||
def gen_data(dt):
|
||||
arrays = {}
|
||||
|
||||
if dt == np.double:
|
||||
pg = './fftw_double'
|
||||
elif dt == np.float32:
|
||||
pg = './fftw_single'
|
||||
else:
|
||||
raise ValueError("unknown: %s" % dt)
|
||||
# Generate test data using FFTW for reference
|
||||
for type in [1, 2, 3, 4, 5, 6, 7, 8]:
|
||||
arrays[type] = {}
|
||||
for sz in SZ:
|
||||
a = Popen([pg, str(type), str(sz)], stdout=PIPE, stderr=STDOUT)
|
||||
st = [i.strip() for i in a.stdout.readlines()]
|
||||
arrays[type][sz] = np.fromstring(",".join(st), sep=',', dtype=dt)
|
||||
|
||||
return arrays
|
||||
|
||||
|
||||
# generate single precision data
|
||||
data = gen_data(np.float32)
|
||||
filename = 'fftw_single_ref'
|
||||
# Save ref data into npz format
|
||||
d = {'sizes': SZ}
|
||||
for type in [1, 2, 3, 4]:
|
||||
for sz in SZ:
|
||||
d['dct_%d_%d' % (type, sz)] = data[type][sz]
|
||||
|
||||
d['sizes'] = SZ
|
||||
for type in [5, 6, 7, 8]:
|
||||
for sz in SZ:
|
||||
d['dst_%d_%d' % (type-4, sz)] = data[type][sz]
|
||||
np.savez(filename, **d)
|
||||
|
||||
|
||||
# generate double precision data
|
||||
data = gen_data(np.float64)
|
||||
filename = 'fftw_double_ref'
|
||||
# Save ref data into npz format
|
||||
d = {'sizes': SZ}
|
||||
for type in [1, 2, 3, 4]:
|
||||
for sz in SZ:
|
||||
d['dct_%d_%d' % (type, sz)] = data[type][sz]
|
||||
|
||||
d['sizes'] = SZ
|
||||
for type in [5, 6, 7, 8]:
|
||||
for sz in SZ:
|
||||
d['dst_%d_%d' % (type-4, sz)] = data[type][sz]
|
||||
np.savez(filename, **d)
|
||||
@@ -0,0 +1,21 @@
|
||||
x0 = linspace(0, 10, 11);
|
||||
x1 = linspace(0, 10, 15);
|
||||
x2 = linspace(0, 10, 16);
|
||||
x3 = linspace(0, 10, 17);
|
||||
|
||||
x4 = randn(32, 1);
|
||||
x5 = randn(64, 1);
|
||||
x6 = randn(128, 1);
|
||||
x7 = randn(256, 1);
|
||||
|
||||
y0 = dct(x0);
|
||||
y1 = dct(x1);
|
||||
y2 = dct(x2);
|
||||
y3 = dct(x3);
|
||||
y4 = dct(x4);
|
||||
y5 = dct(x5);
|
||||
y6 = dct(x6);
|
||||
y7 = dct(x7);
|
||||
|
||||
save('test.mat', 'x0', 'x1', 'x2', 'x3', 'x4', 'x5', 'x6', 'x7', ...
|
||||
'y0', 'y1', 'y2', 'y3', 'y4', 'y5', 'y6', 'y7');
|
||||
@@ -0,0 +1,8 @@
|
||||
from __future__ import division, print_function, absolute_import
|
||||
|
||||
import numpy as np
|
||||
from scipy.io import loadmat
|
||||
|
||||
m = loadmat('test.mat', squeeze_me=True, struct_as_record=True,
|
||||
mat_dtype=True)
|
||||
np.savez('test.npz', **m)
|
||||
Binary file not shown.
@@ -0,0 +1,972 @@
|
||||
# Created by Pearu Peterson, September 2002
|
||||
|
||||
from __future__ import division, print_function, absolute_import
|
||||
|
||||
__usage__ = """
|
||||
Build fftpack:
|
||||
python setup_fftpack.py build
|
||||
Run tests if scipy is installed:
|
||||
python -c 'import scipy;scipy.fftpack.test()'
|
||||
Run tests if fftpack is not installed:
|
||||
python tests/test_basic.py
|
||||
"""
|
||||
|
||||
from numpy.testing import (assert_, assert_equal, assert_array_almost_equal,
|
||||
assert_array_almost_equal_nulp, assert_array_less)
|
||||
import pytest
|
||||
from pytest import raises as assert_raises
|
||||
from scipy.fftpack import ifft, fft, fftn, ifftn, rfft, irfft, fft2
|
||||
from scipy.fftpack import _fftpack as fftpack
|
||||
from scipy.fftpack.basic import _is_safe_size
|
||||
|
||||
from numpy import (arange, add, array, asarray, zeros, dot, exp, pi,
|
||||
swapaxes, double, cdouble)
|
||||
import numpy as np
|
||||
import numpy.fft
|
||||
from numpy.random import rand
|
||||
|
||||
# "large" composite numbers supported by FFTPACK
|
||||
LARGE_COMPOSITE_SIZES = [
|
||||
2**13,
|
||||
2**5 * 3**5,
|
||||
2**3 * 3**3 * 5**2,
|
||||
]
|
||||
SMALL_COMPOSITE_SIZES = [
|
||||
2,
|
||||
2*3*5,
|
||||
2*2*3*3,
|
||||
]
|
||||
# prime
|
||||
LARGE_PRIME_SIZES = [
|
||||
2011
|
||||
]
|
||||
SMALL_PRIME_SIZES = [
|
||||
29
|
||||
]
|
||||
|
||||
|
||||
def _assert_close_in_norm(x, y, rtol, size, rdt):
|
||||
# helper function for testing
|
||||
err_msg = "size: %s rdt: %s" % (size, rdt)
|
||||
assert_array_less(np.linalg.norm(x - y), rtol*np.linalg.norm(x), err_msg)
|
||||
|
||||
|
||||
def random(size):
|
||||
return rand(*size)
|
||||
|
||||
|
||||
def get_mat(n):
|
||||
data = arange(n)
|
||||
data = add.outer(data, data)
|
||||
return data
|
||||
|
||||
|
||||
def direct_dft(x):
|
||||
x = asarray(x)
|
||||
n = len(x)
|
||||
y = zeros(n, dtype=cdouble)
|
||||
w = -arange(n)*(2j*pi/n)
|
||||
for i in range(n):
|
||||
y[i] = dot(exp(i*w), x)
|
||||
return y
|
||||
|
||||
|
||||
def direct_idft(x):
|
||||
x = asarray(x)
|
||||
n = len(x)
|
||||
y = zeros(n, dtype=cdouble)
|
||||
w = arange(n)*(2j*pi/n)
|
||||
for i in range(n):
|
||||
y[i] = dot(exp(i*w), x)/n
|
||||
return y
|
||||
|
||||
|
||||
def direct_dftn(x):
|
||||
x = asarray(x)
|
||||
for axis in range(len(x.shape)):
|
||||
x = fft(x, axis=axis)
|
||||
return x
|
||||
|
||||
|
||||
def direct_idftn(x):
|
||||
x = asarray(x)
|
||||
for axis in range(len(x.shape)):
|
||||
x = ifft(x, axis=axis)
|
||||
return x
|
||||
|
||||
|
||||
def direct_rdft(x):
|
||||
x = asarray(x)
|
||||
n = len(x)
|
||||
w = -arange(n)*(2j*pi/n)
|
||||
r = zeros(n, dtype=double)
|
||||
for i in range(n//2+1):
|
||||
y = dot(exp(i*w), x)
|
||||
if i:
|
||||
r[2*i-1] = y.real
|
||||
if 2*i < n:
|
||||
r[2*i] = y.imag
|
||||
else:
|
||||
r[0] = y.real
|
||||
return r
|
||||
|
||||
|
||||
def direct_irdft(x):
|
||||
x = asarray(x)
|
||||
n = len(x)
|
||||
x1 = zeros(n, dtype=cdouble)
|
||||
for i in range(n//2+1):
|
||||
if i:
|
||||
if 2*i < n:
|
||||
x1[i] = x[2*i-1] + 1j*x[2*i]
|
||||
x1[n-i] = x[2*i-1] - 1j*x[2*i]
|
||||
else:
|
||||
x1[i] = x[2*i-1]
|
||||
else:
|
||||
x1[0] = x[0]
|
||||
return direct_idft(x1).real
|
||||
|
||||
|
||||
class _TestFFTBase(object):
|
||||
def setup_method(self):
|
||||
self.cdt = None
|
||||
self.rdt = None
|
||||
np.random.seed(1234)
|
||||
|
||||
def test_definition(self):
|
||||
x = np.array([1,2,3,4+1j,1,2,3,4+2j], dtype=self.cdt)
|
||||
y = fft(x)
|
||||
assert_equal(y.dtype, self.cdt)
|
||||
y1 = direct_dft(x)
|
||||
assert_array_almost_equal(y,y1)
|
||||
x = np.array([1,2,3,4+0j,5], dtype=self.cdt)
|
||||
assert_array_almost_equal(fft(x),direct_dft(x))
|
||||
|
||||
def test_n_argument_real(self):
|
||||
x1 = np.array([1,2,3,4], dtype=self.rdt)
|
||||
x2 = np.array([1,2,3,4], dtype=self.rdt)
|
||||
y = fft([x1,x2],n=4)
|
||||
assert_equal(y.dtype, self.cdt)
|
||||
assert_equal(y.shape,(2,4))
|
||||
assert_array_almost_equal(y[0],direct_dft(x1))
|
||||
assert_array_almost_equal(y[1],direct_dft(x2))
|
||||
|
||||
def _test_n_argument_complex(self):
|
||||
x1 = np.array([1,2,3,4+1j], dtype=self.cdt)
|
||||
x2 = np.array([1,2,3,4+1j], dtype=self.cdt)
|
||||
y = fft([x1,x2],n=4)
|
||||
assert_equal(y.dtype, self.cdt)
|
||||
assert_equal(y.shape,(2,4))
|
||||
assert_array_almost_equal(y[0],direct_dft(x1))
|
||||
assert_array_almost_equal(y[1],direct_dft(x2))
|
||||
|
||||
def test_djbfft(self):
|
||||
for i in range(2,14):
|
||||
n = 2**i
|
||||
x = list(range(n))
|
||||
y = fftpack.zfft(x)
|
||||
y2 = numpy.fft.fft(x)
|
||||
assert_array_almost_equal(y,y2)
|
||||
y = fftpack.zrfft(x)
|
||||
assert_array_almost_equal(y,y2)
|
||||
|
||||
def test_invalid_sizes(self):
|
||||
assert_raises(ValueError, fft, [])
|
||||
assert_raises(ValueError, fft, [[1,1],[2,2]], -5)
|
||||
|
||||
def test__is_safe_size(self):
|
||||
vals = [(0, True), (1, True), (2, True), (3, True), (4, True), (5, True), (6, True), (7, False),
|
||||
(15, True), (16, True), (17, False), (18, True), (21, False), (25, True), (50, True),
|
||||
(120, True), (210, False)]
|
||||
for n, is_safe in vals:
|
||||
assert_equal(_is_safe_size(n), is_safe)
|
||||
|
||||
|
||||
class TestDoubleFFT(_TestFFTBase):
|
||||
def setup_method(self):
|
||||
self.cdt = np.cdouble
|
||||
self.rdt = np.double
|
||||
|
||||
|
||||
class TestSingleFFT(_TestFFTBase):
|
||||
def setup_method(self):
|
||||
self.cdt = np.complex64
|
||||
self.rdt = np.float32
|
||||
|
||||
@pytest.mark.xfail(run=False, reason="single-precision FFT implementation is partially disabled, until accuracy issues with large prime powers are resolved")
|
||||
def test_notice(self):
|
||||
pass
|
||||
|
||||
|
||||
class TestFloat16FFT(object):
|
||||
|
||||
def test_1_argument_real(self):
|
||||
x1 = np.array([1, 2, 3, 4], dtype=np.float16)
|
||||
y = fft(x1, n=4)
|
||||
assert_equal(y.dtype, np.complex64)
|
||||
assert_equal(y.shape, (4, ))
|
||||
assert_array_almost_equal(y, direct_dft(x1.astype(np.float32)))
|
||||
|
||||
def test_n_argument_real(self):
|
||||
x1 = np.array([1, 2, 3, 4], dtype=np.float16)
|
||||
x2 = np.array([1, 2, 3, 4], dtype=np.float16)
|
||||
y = fft([x1, x2], n=4)
|
||||
assert_equal(y.dtype, np.complex64)
|
||||
assert_equal(y.shape, (2, 4))
|
||||
assert_array_almost_equal(y[0], direct_dft(x1.astype(np.float32)))
|
||||
assert_array_almost_equal(y[1], direct_dft(x2.astype(np.float32)))
|
||||
|
||||
|
||||
class _TestIFFTBase(object):
|
||||
def setup_method(self):
|
||||
np.random.seed(1234)
|
||||
|
||||
def test_definition(self):
|
||||
x = np.array([1,2,3,4+1j,1,2,3,4+2j], self.cdt)
|
||||
y = ifft(x)
|
||||
y1 = direct_idft(x)
|
||||
assert_equal(y.dtype, self.cdt)
|
||||
assert_array_almost_equal(y,y1)
|
||||
|
||||
x = np.array([1,2,3,4+0j,5], self.cdt)
|
||||
assert_array_almost_equal(ifft(x),direct_idft(x))
|
||||
|
||||
def test_definition_real(self):
|
||||
x = np.array([1,2,3,4,1,2,3,4], self.rdt)
|
||||
y = ifft(x)
|
||||
assert_equal(y.dtype, self.cdt)
|
||||
y1 = direct_idft(x)
|
||||
assert_array_almost_equal(y,y1)
|
||||
|
||||
x = np.array([1,2,3,4,5], dtype=self.rdt)
|
||||
assert_equal(y.dtype, self.cdt)
|
||||
assert_array_almost_equal(ifft(x),direct_idft(x))
|
||||
|
||||
def test_djbfft(self):
|
||||
for i in range(2,14):
|
||||
n = 2**i
|
||||
x = list(range(n))
|
||||
y = fftpack.zfft(x,direction=-1)
|
||||
y2 = numpy.fft.ifft(x)
|
||||
assert_array_almost_equal(y,y2)
|
||||
y = fftpack.zrfft(x,direction=-1)
|
||||
assert_array_almost_equal(y,y2)
|
||||
|
||||
def test_random_complex(self):
|
||||
for size in [1,51,111,100,200,64,128,256,1024]:
|
||||
x = random([size]).astype(self.cdt)
|
||||
x = random([size]).astype(self.cdt) + 1j*x
|
||||
y1 = ifft(fft(x))
|
||||
y2 = fft(ifft(x))
|
||||
assert_equal(y1.dtype, self.cdt)
|
||||
assert_equal(y2.dtype, self.cdt)
|
||||
assert_array_almost_equal(y1, x)
|
||||
assert_array_almost_equal(y2, x)
|
||||
|
||||
def test_random_real(self):
|
||||
for size in [1,51,111,100,200,64,128,256,1024]:
|
||||
x = random([size]).astype(self.rdt)
|
||||
y1 = ifft(fft(x))
|
||||
y2 = fft(ifft(x))
|
||||
assert_equal(y1.dtype, self.cdt)
|
||||
assert_equal(y2.dtype, self.cdt)
|
||||
assert_array_almost_equal(y1, x)
|
||||
assert_array_almost_equal(y2, x)
|
||||
|
||||
def test_size_accuracy(self):
|
||||
# Sanity check for the accuracy for prime and non-prime sized inputs
|
||||
if self.rdt == np.float32:
|
||||
rtol = 1e-5
|
||||
elif self.rdt == np.float64:
|
||||
rtol = 1e-10
|
||||
|
||||
for size in LARGE_COMPOSITE_SIZES + LARGE_PRIME_SIZES:
|
||||
np.random.seed(1234)
|
||||
x = np.random.rand(size).astype(self.rdt)
|
||||
y = ifft(fft(x))
|
||||
_assert_close_in_norm(x, y, rtol, size, self.rdt)
|
||||
y = fft(ifft(x))
|
||||
_assert_close_in_norm(x, y, rtol, size, self.rdt)
|
||||
|
||||
x = (x + 1j*np.random.rand(size)).astype(self.cdt)
|
||||
y = ifft(fft(x))
|
||||
_assert_close_in_norm(x, y, rtol, size, self.rdt)
|
||||
y = fft(ifft(x))
|
||||
_assert_close_in_norm(x, y, rtol, size, self.rdt)
|
||||
|
||||
def test_invalid_sizes(self):
|
||||
assert_raises(ValueError, ifft, [])
|
||||
assert_raises(ValueError, ifft, [[1,1],[2,2]], -5)
|
||||
|
||||
|
||||
class TestDoubleIFFT(_TestIFFTBase):
|
||||
def setup_method(self):
|
||||
self.cdt = np.cdouble
|
||||
self.rdt = np.double
|
||||
|
||||
|
||||
class TestSingleIFFT(_TestIFFTBase):
|
||||
def setup_method(self):
|
||||
self.cdt = np.complex64
|
||||
self.rdt = np.float32
|
||||
|
||||
|
||||
class _TestRFFTBase(object):
|
||||
def setup_method(self):
|
||||
np.random.seed(1234)
|
||||
|
||||
def test_definition(self):
|
||||
for t in [[1, 2, 3, 4, 1, 2, 3, 4], [1, 2, 3, 4, 1, 2, 3, 4, 5]]:
|
||||
x = np.array(t, dtype=self.rdt)
|
||||
y = rfft(x)
|
||||
y1 = direct_rdft(x)
|
||||
assert_array_almost_equal(y,y1)
|
||||
assert_equal(y.dtype, self.rdt)
|
||||
|
||||
def test_djbfft(self):
|
||||
from numpy.fft import fft as numpy_fft
|
||||
for i in range(2,14):
|
||||
n = 2**i
|
||||
x = list(range(n))
|
||||
y2 = numpy_fft(x)
|
||||
y1 = zeros((n,),dtype=double)
|
||||
y1[0] = y2[0].real
|
||||
y1[-1] = y2[n//2].real
|
||||
for k in range(1, n//2):
|
||||
y1[2*k-1] = y2[k].real
|
||||
y1[2*k] = y2[k].imag
|
||||
y = fftpack.drfft(x)
|
||||
assert_array_almost_equal(y,y1)
|
||||
|
||||
def test_invalid_sizes(self):
|
||||
assert_raises(ValueError, rfft, [])
|
||||
assert_raises(ValueError, rfft, [[1,1],[2,2]], -5)
|
||||
|
||||
# See gh-5790
|
||||
class MockSeries(object):
|
||||
def __init__(self, data):
|
||||
self.data = np.asarray(data)
|
||||
|
||||
def __getattr__(self, item):
|
||||
try:
|
||||
return getattr(self.data, item)
|
||||
except AttributeError:
|
||||
raise AttributeError(("'MockSeries' object "
|
||||
"has no attribute '{attr}'".
|
||||
format(attr=item)))
|
||||
|
||||
def test_non_ndarray_with_dtype(self):
|
||||
x = np.array([1., 2., 3., 4., 5.])
|
||||
xs = _TestRFFTBase.MockSeries(x)
|
||||
|
||||
expected = [1, 2, 3, 4, 5]
|
||||
out = rfft(xs)
|
||||
|
||||
# Data should not have been overwritten
|
||||
assert_equal(x, expected)
|
||||
assert_equal(xs.data, expected)
|
||||
|
||||
class TestRFFTDouble(_TestRFFTBase):
|
||||
def setup_method(self):
|
||||
self.cdt = np.cdouble
|
||||
self.rdt = np.double
|
||||
|
||||
|
||||
class TestRFFTSingle(_TestRFFTBase):
|
||||
def setup_method(self):
|
||||
self.cdt = np.complex64
|
||||
self.rdt = np.float32
|
||||
|
||||
|
||||
class _TestIRFFTBase(object):
|
||||
def setup_method(self):
|
||||
np.random.seed(1234)
|
||||
|
||||
def test_definition(self):
|
||||
x1 = [1,2,3,4,1,2,3,4]
|
||||
x1_1 = [1,2+3j,4+1j,2+3j,4,2-3j,4-1j,2-3j]
|
||||
x2 = [1,2,3,4,1,2,3,4,5]
|
||||
x2_1 = [1,2+3j,4+1j,2+3j,4+5j,4-5j,2-3j,4-1j,2-3j]
|
||||
|
||||
def _test(x, xr):
|
||||
y = irfft(np.array(x, dtype=self.rdt))
|
||||
y1 = direct_irdft(x)
|
||||
assert_equal(y.dtype, self.rdt)
|
||||
assert_array_almost_equal(y,y1, decimal=self.ndec)
|
||||
assert_array_almost_equal(y,ifft(xr), decimal=self.ndec)
|
||||
|
||||
_test(x1, x1_1)
|
||||
_test(x2, x2_1)
|
||||
|
||||
def test_djbfft(self):
|
||||
from numpy.fft import ifft as numpy_ifft
|
||||
for i in range(2,14):
|
||||
n = 2**i
|
||||
x = list(range(n))
|
||||
x1 = zeros((n,),dtype=cdouble)
|
||||
x1[0] = x[0]
|
||||
for k in range(1, n//2):
|
||||
x1[k] = x[2*k-1]+1j*x[2*k]
|
||||
x1[n-k] = x[2*k-1]-1j*x[2*k]
|
||||
x1[n//2] = x[-1]
|
||||
y1 = numpy_ifft(x1)
|
||||
y = fftpack.drfft(x,direction=-1)
|
||||
assert_array_almost_equal(y,y1)
|
||||
|
||||
def test_random_real(self):
|
||||
for size in [1,51,111,100,200,64,128,256,1024]:
|
||||
x = random([size]).astype(self.rdt)
|
||||
y1 = irfft(rfft(x))
|
||||
y2 = rfft(irfft(x))
|
||||
assert_equal(y1.dtype, self.rdt)
|
||||
assert_equal(y2.dtype, self.rdt)
|
||||
assert_array_almost_equal(y1, x, decimal=self.ndec,
|
||||
err_msg="size=%d" % size)
|
||||
assert_array_almost_equal(y2, x, decimal=self.ndec,
|
||||
err_msg="size=%d" % size)
|
||||
|
||||
def test_size_accuracy(self):
|
||||
# Sanity check for the accuracy for prime and non-prime sized inputs
|
||||
if self.rdt == np.float32:
|
||||
rtol = 1e-5
|
||||
elif self.rdt == np.float64:
|
||||
rtol = 1e-10
|
||||
|
||||
for size in LARGE_COMPOSITE_SIZES + LARGE_PRIME_SIZES:
|
||||
np.random.seed(1234)
|
||||
x = np.random.rand(size).astype(self.rdt)
|
||||
y = irfft(rfft(x))
|
||||
_assert_close_in_norm(x, y, rtol, size, self.rdt)
|
||||
y = rfft(irfft(x))
|
||||
_assert_close_in_norm(x, y, rtol, size, self.rdt)
|
||||
|
||||
def test_invalid_sizes(self):
|
||||
assert_raises(ValueError, irfft, [])
|
||||
assert_raises(ValueError, irfft, [[1,1],[2,2]], -5)
|
||||
|
||||
|
||||
# self.ndec is bogus; we should have a assert_array_approx_equal for number of
|
||||
# significant digits
|
||||
|
||||
class TestIRFFTDouble(_TestIRFFTBase):
|
||||
def setup_method(self):
|
||||
self.cdt = np.cdouble
|
||||
self.rdt = np.double
|
||||
self.ndec = 14
|
||||
|
||||
|
||||
class TestIRFFTSingle(_TestIRFFTBase):
|
||||
def setup_method(self):
|
||||
self.cdt = np.complex64
|
||||
self.rdt = np.float32
|
||||
self.ndec = 5
|
||||
|
||||
|
||||
class Testfft2(object):
|
||||
def setup_method(self):
|
||||
np.random.seed(1234)
|
||||
|
||||
def test_regression_244(self):
|
||||
"""FFT returns wrong result with axes parameter."""
|
||||
# fftn (and hence fft2) used to break when both axes and shape were
|
||||
# used
|
||||
x = numpy.ones((4, 4, 2))
|
||||
y = fft2(x, shape=(8, 8), axes=(-3, -2))
|
||||
y_r = numpy.fft.fftn(x, s=(8, 8), axes=(-3, -2))
|
||||
assert_array_almost_equal(y, y_r)
|
||||
|
||||
def test_invalid_sizes(self):
|
||||
assert_raises(ValueError, fft2, [[]])
|
||||
assert_raises(ValueError, fft2, [[1, 1], [2, 2]], (4, -3))
|
||||
|
||||
|
||||
class TestFftnSingle(object):
|
||||
def setup_method(self):
|
||||
np.random.seed(1234)
|
||||
|
||||
def test_definition(self):
|
||||
x = [[1, 2, 3],
|
||||
[4, 5, 6],
|
||||
[7, 8, 9]]
|
||||
y = fftn(np.array(x, np.float32))
|
||||
assert_(y.dtype == np.complex64,
|
||||
msg="double precision output with single precision")
|
||||
|
||||
y_r = np.array(fftn(x), np.complex64)
|
||||
assert_array_almost_equal_nulp(y, y_r)
|
||||
|
||||
@pytest.mark.parametrize('size', SMALL_COMPOSITE_SIZES + SMALL_PRIME_SIZES)
|
||||
def test_size_accuracy_small(self, size):
|
||||
x = np.random.rand(size, size) + 1j*np.random.rand(size, size)
|
||||
y1 = fftn(x.real.astype(np.float32))
|
||||
y2 = fftn(x.real.astype(np.float64)).astype(np.complex64)
|
||||
|
||||
assert_equal(y1.dtype, np.complex64)
|
||||
assert_array_almost_equal_nulp(y1, y2, 2000)
|
||||
|
||||
@pytest.mark.parametrize('size', LARGE_COMPOSITE_SIZES + LARGE_PRIME_SIZES)
|
||||
def test_size_accuracy_large(self, size):
|
||||
x = np.random.rand(size, 3) + 1j*np.random.rand(size, 3)
|
||||
y1 = fftn(x.real.astype(np.float32))
|
||||
y2 = fftn(x.real.astype(np.float64)).astype(np.complex64)
|
||||
|
||||
assert_equal(y1.dtype, np.complex64)
|
||||
assert_array_almost_equal_nulp(y1, y2, 2000)
|
||||
|
||||
def test_definition_float16(self):
|
||||
x = [[1, 2, 3],
|
||||
[4, 5, 6],
|
||||
[7, 8, 9]]
|
||||
y = fftn(np.array(x, np.float16))
|
||||
assert_equal(y.dtype, np.complex64)
|
||||
y_r = np.array(fftn(x), np.complex64)
|
||||
assert_array_almost_equal_nulp(y, y_r)
|
||||
|
||||
@pytest.mark.parametrize('size', SMALL_COMPOSITE_SIZES + SMALL_PRIME_SIZES)
|
||||
def test_float16_input_small(self, size):
|
||||
x = np.random.rand(size, size) + 1j*np.random.rand(size, size)
|
||||
y1 = fftn(x.real.astype(np.float16))
|
||||
y2 = fftn(x.real.astype(np.float64)).astype(np.complex64)
|
||||
|
||||
assert_equal(y1.dtype, np.complex64)
|
||||
assert_array_almost_equal_nulp(y1, y2, 5e5)
|
||||
|
||||
@pytest.mark.parametrize('size', LARGE_COMPOSITE_SIZES + LARGE_PRIME_SIZES)
|
||||
def test_float16_input_large(self, size):
|
||||
x = np.random.rand(size, 3) + 1j*np.random.rand(size, 3)
|
||||
y1 = fftn(x.real.astype(np.float16))
|
||||
y2 = fftn(x.real.astype(np.float64)).astype(np.complex64)
|
||||
|
||||
assert_equal(y1.dtype, np.complex64)
|
||||
assert_array_almost_equal_nulp(y1, y2, 2e6)
|
||||
|
||||
|
||||
class TestFftn(object):
|
||||
def setup_method(self):
|
||||
np.random.seed(1234)
|
||||
|
||||
def test_definition(self):
|
||||
x = [[1, 2, 3],
|
||||
[4, 5, 6],
|
||||
[7, 8, 9]]
|
||||
y = fftn(x)
|
||||
assert_array_almost_equal(y, direct_dftn(x))
|
||||
|
||||
x = random((20, 26))
|
||||
assert_array_almost_equal(fftn(x), direct_dftn(x))
|
||||
|
||||
x = random((5, 4, 3, 20))
|
||||
assert_array_almost_equal(fftn(x), direct_dftn(x))
|
||||
|
||||
def test_axes_argument(self):
|
||||
# plane == ji_plane, x== kji_space
|
||||
plane1 = [[1, 2, 3],
|
||||
[4, 5, 6],
|
||||
[7, 8, 9]]
|
||||
plane2 = [[10, 11, 12],
|
||||
[13, 14, 15],
|
||||
[16, 17, 18]]
|
||||
plane3 = [[19, 20, 21],
|
||||
[22, 23, 24],
|
||||
[25, 26, 27]]
|
||||
ki_plane1 = [[1, 2, 3],
|
||||
[10, 11, 12],
|
||||
[19, 20, 21]]
|
||||
ki_plane2 = [[4, 5, 6],
|
||||
[13, 14, 15],
|
||||
[22, 23, 24]]
|
||||
ki_plane3 = [[7, 8, 9],
|
||||
[16, 17, 18],
|
||||
[25, 26, 27]]
|
||||
jk_plane1 = [[1, 10, 19],
|
||||
[4, 13, 22],
|
||||
[7, 16, 25]]
|
||||
jk_plane2 = [[2, 11, 20],
|
||||
[5, 14, 23],
|
||||
[8, 17, 26]]
|
||||
jk_plane3 = [[3, 12, 21],
|
||||
[6, 15, 24],
|
||||
[9, 18, 27]]
|
||||
kj_plane1 = [[1, 4, 7],
|
||||
[10, 13, 16], [19, 22, 25]]
|
||||
kj_plane2 = [[2, 5, 8],
|
||||
[11, 14, 17], [20, 23, 26]]
|
||||
kj_plane3 = [[3, 6, 9],
|
||||
[12, 15, 18], [21, 24, 27]]
|
||||
ij_plane1 = [[1, 4, 7],
|
||||
[2, 5, 8],
|
||||
[3, 6, 9]]
|
||||
ij_plane2 = [[10, 13, 16],
|
||||
[11, 14, 17],
|
||||
[12, 15, 18]]
|
||||
ij_plane3 = [[19, 22, 25],
|
||||
[20, 23, 26],
|
||||
[21, 24, 27]]
|
||||
ik_plane1 = [[1, 10, 19],
|
||||
[2, 11, 20],
|
||||
[3, 12, 21]]
|
||||
ik_plane2 = [[4, 13, 22],
|
||||
[5, 14, 23],
|
||||
[6, 15, 24]]
|
||||
ik_plane3 = [[7, 16, 25],
|
||||
[8, 17, 26],
|
||||
[9, 18, 27]]
|
||||
ijk_space = [jk_plane1, jk_plane2, jk_plane3]
|
||||
ikj_space = [kj_plane1, kj_plane2, kj_plane3]
|
||||
jik_space = [ik_plane1, ik_plane2, ik_plane3]
|
||||
jki_space = [ki_plane1, ki_plane2, ki_plane3]
|
||||
kij_space = [ij_plane1, ij_plane2, ij_plane3]
|
||||
x = array([plane1, plane2, plane3])
|
||||
|
||||
assert_array_almost_equal(fftn(x),
|
||||
fftn(x, axes=(-3, -2, -1))) # kji_space
|
||||
assert_array_almost_equal(fftn(x), fftn(x, axes=(0, 1, 2)))
|
||||
assert_array_almost_equal(fftn(x, axes=(0, 2)), fftn(x, axes=(0, -1)))
|
||||
y = fftn(x, axes=(2, 1, 0)) # ijk_space
|
||||
assert_array_almost_equal(swapaxes(y, -1, -3), fftn(ijk_space))
|
||||
y = fftn(x, axes=(2, 0, 1)) # ikj_space
|
||||
assert_array_almost_equal(swapaxes(swapaxes(y, -1, -3), -1, -2),
|
||||
fftn(ikj_space))
|
||||
y = fftn(x, axes=(1, 2, 0)) # jik_space
|
||||
assert_array_almost_equal(swapaxes(swapaxes(y, -1, -3), -3, -2),
|
||||
fftn(jik_space))
|
||||
y = fftn(x, axes=(1, 0, 2)) # jki_space
|
||||
assert_array_almost_equal(swapaxes(y, -2, -3), fftn(jki_space))
|
||||
y = fftn(x, axes=(0, 2, 1)) # kij_space
|
||||
assert_array_almost_equal(swapaxes(y, -2, -1), fftn(kij_space))
|
||||
|
||||
y = fftn(x, axes=(-2, -1)) # ji_plane
|
||||
assert_array_almost_equal(fftn(plane1), y[0])
|
||||
assert_array_almost_equal(fftn(plane2), y[1])
|
||||
assert_array_almost_equal(fftn(plane3), y[2])
|
||||
|
||||
y = fftn(x, axes=(1, 2)) # ji_plane
|
||||
assert_array_almost_equal(fftn(plane1), y[0])
|
||||
assert_array_almost_equal(fftn(plane2), y[1])
|
||||
assert_array_almost_equal(fftn(plane3), y[2])
|
||||
|
||||
y = fftn(x, axes=(-3, -2)) # kj_plane
|
||||
assert_array_almost_equal(fftn(x[:, :, 0]), y[:, :, 0])
|
||||
assert_array_almost_equal(fftn(x[:, :, 1]), y[:, :, 1])
|
||||
assert_array_almost_equal(fftn(x[:, :, 2]), y[:, :, 2])
|
||||
|
||||
y = fftn(x, axes=(-3, -1)) # ki_plane
|
||||
assert_array_almost_equal(fftn(x[:, 0, :]), y[:, 0, :])
|
||||
assert_array_almost_equal(fftn(x[:, 1, :]), y[:, 1, :])
|
||||
assert_array_almost_equal(fftn(x[:, 2, :]), y[:, 2, :])
|
||||
|
||||
y = fftn(x, axes=(-1, -2)) # ij_plane
|
||||
assert_array_almost_equal(fftn(ij_plane1), swapaxes(y[0], -2, -1))
|
||||
assert_array_almost_equal(fftn(ij_plane2), swapaxes(y[1], -2, -1))
|
||||
assert_array_almost_equal(fftn(ij_plane3), swapaxes(y[2], -2, -1))
|
||||
|
||||
y = fftn(x, axes=(-1, -3)) # ik_plane
|
||||
assert_array_almost_equal(fftn(ik_plane1),
|
||||
swapaxes(y[:, 0, :], -1, -2))
|
||||
assert_array_almost_equal(fftn(ik_plane2),
|
||||
swapaxes(y[:, 1, :], -1, -2))
|
||||
assert_array_almost_equal(fftn(ik_plane3),
|
||||
swapaxes(y[:, 2, :], -1, -2))
|
||||
|
||||
y = fftn(x, axes=(-2, -3)) # jk_plane
|
||||
assert_array_almost_equal(fftn(jk_plane1),
|
||||
swapaxes(y[:, :, 0], -1, -2))
|
||||
assert_array_almost_equal(fftn(jk_plane2),
|
||||
swapaxes(y[:, :, 1], -1, -2))
|
||||
assert_array_almost_equal(fftn(jk_plane3),
|
||||
swapaxes(y[:, :, 2], -1, -2))
|
||||
|
||||
y = fftn(x, axes=(-1,)) # i_line
|
||||
for i in range(3):
|
||||
for j in range(3):
|
||||
assert_array_almost_equal(fft(x[i, j, :]), y[i, j, :])
|
||||
y = fftn(x, axes=(-2,)) # j_line
|
||||
for i in range(3):
|
||||
for j in range(3):
|
||||
assert_array_almost_equal(fft(x[i, :, j]), y[i, :, j])
|
||||
y = fftn(x, axes=(0,)) # k_line
|
||||
for i in range(3):
|
||||
for j in range(3):
|
||||
assert_array_almost_equal(fft(x[:, i, j]), y[:, i, j])
|
||||
|
||||
y = fftn(x, axes=()) # point
|
||||
assert_array_almost_equal(y, x)
|
||||
|
||||
def test_shape_argument(self):
|
||||
small_x = [[1, 2, 3],
|
||||
[4, 5, 6]]
|
||||
large_x1 = [[1, 2, 3, 0],
|
||||
[4, 5, 6, 0],
|
||||
[0, 0, 0, 0],
|
||||
[0, 0, 0, 0]]
|
||||
|
||||
y = fftn(small_x, shape=(4, 4))
|
||||
assert_array_almost_equal(y, fftn(large_x1))
|
||||
|
||||
y = fftn(small_x, shape=(3, 4))
|
||||
assert_array_almost_equal(y, fftn(large_x1[:-1]))
|
||||
|
||||
def test_shape_axes_argument(self):
|
||||
small_x = [[1, 2, 3],
|
||||
[4, 5, 6],
|
||||
[7, 8, 9]]
|
||||
large_x1 = array([[1, 2, 3, 0],
|
||||
[4, 5, 6, 0],
|
||||
[7, 8, 9, 0],
|
||||
[0, 0, 0, 0]])
|
||||
y = fftn(small_x, shape=(4, 4), axes=(-2, -1))
|
||||
assert_array_almost_equal(y, fftn(large_x1))
|
||||
y = fftn(small_x, shape=(4, 4), axes=(-1, -2))
|
||||
|
||||
assert_array_almost_equal(y, swapaxes(
|
||||
fftn(swapaxes(large_x1, -1, -2)), -1, -2))
|
||||
|
||||
def test_shape_axes_argument2(self):
|
||||
# Change shape of the last axis
|
||||
x = numpy.random.random((10, 5, 3, 7))
|
||||
y = fftn(x, axes=(-1,), shape=(8,))
|
||||
assert_array_almost_equal(y, fft(x, axis=-1, n=8))
|
||||
|
||||
# Change shape of an arbitrary axis which is not the last one
|
||||
x = numpy.random.random((10, 5, 3, 7))
|
||||
y = fftn(x, axes=(-2,), shape=(8,))
|
||||
assert_array_almost_equal(y, fft(x, axis=-2, n=8))
|
||||
|
||||
# Change shape of axes: cf #244, where shape and axes were mixed up
|
||||
x = numpy.random.random((4, 4, 2))
|
||||
y = fftn(x, axes=(-3, -2), shape=(8, 8))
|
||||
assert_array_almost_equal(y,
|
||||
numpy.fft.fftn(x, axes=(-3, -2), s=(8, 8)))
|
||||
|
||||
def test_shape_argument_more(self):
|
||||
x = zeros((4, 4, 2))
|
||||
with assert_raises(ValueError,
|
||||
match="when given, axes and shape arguments"
|
||||
" have to be of the same length"):
|
||||
fftn(x, shape=(8, 8, 2, 1))
|
||||
|
||||
def test_invalid_sizes(self):
|
||||
with assert_raises(ValueError,
|
||||
match="invalid number of data points"
|
||||
r" \(\[1 0\]\) specified"):
|
||||
fftn([[]])
|
||||
|
||||
with assert_raises(ValueError,
|
||||
match="invalid number of data points"
|
||||
r" \(\[ 4 -3\]\) specified"):
|
||||
fftn([[1, 1], [2, 2]], (4, -3))
|
||||
|
||||
|
||||
class TestIfftn(object):
|
||||
dtype = None
|
||||
cdtype = None
|
||||
|
||||
def setup_method(self):
|
||||
np.random.seed(1234)
|
||||
|
||||
@pytest.mark.parametrize('dtype,cdtype,maxnlp',
|
||||
[(np.float64, np.complex128, 2000),
|
||||
(np.float32, np.complex64, 3500)])
|
||||
def test_definition(self, dtype, cdtype, maxnlp):
|
||||
x = np.array([[1, 2, 3],
|
||||
[4, 5, 6],
|
||||
[7, 8, 9]], dtype=dtype)
|
||||
y = ifftn(x)
|
||||
assert_equal(y.dtype, cdtype)
|
||||
assert_array_almost_equal_nulp(y, direct_idftn(x), maxnlp)
|
||||
|
||||
x = random((20, 26))
|
||||
assert_array_almost_equal_nulp(ifftn(x), direct_idftn(x), maxnlp)
|
||||
|
||||
x = random((5, 4, 3, 20))
|
||||
assert_array_almost_equal_nulp(ifftn(x), direct_idftn(x), maxnlp)
|
||||
|
||||
@pytest.mark.parametrize('maxnlp', [2000, 3500])
|
||||
@pytest.mark.parametrize('size', [1, 2, 51, 32, 64, 92])
|
||||
def test_random_complex(self, maxnlp, size):
|
||||
x = random([size, size]) + 1j*random([size, size])
|
||||
assert_array_almost_equal_nulp(ifftn(fftn(x)), x, maxnlp)
|
||||
assert_array_almost_equal_nulp(fftn(ifftn(x)), x, maxnlp)
|
||||
|
||||
def test_invalid_sizes(self):
|
||||
with assert_raises(ValueError,
|
||||
match="invalid number of data points"
|
||||
r" \(\[1 0\]\) specified"):
|
||||
ifftn([[]])
|
||||
|
||||
with assert_raises(ValueError,
|
||||
match="invalid number of data points"
|
||||
r" \(\[ 4 -3\]\) specified"):
|
||||
ifftn([[1, 1], [2, 2]], (4, -3))
|
||||
|
||||
|
||||
class TestLongDoubleFailure(object):
|
||||
def setup_method(self):
|
||||
np.random.seed(1234)
|
||||
|
||||
def test_complex(self):
|
||||
if np.dtype(np.longcomplex).itemsize == np.dtype(complex).itemsize:
|
||||
# longdouble == double; so fft is supported
|
||||
return
|
||||
|
||||
x = np.random.randn(10).astype(np.longdouble) + \
|
||||
1j * np.random.randn(10).astype(np.longdouble)
|
||||
|
||||
for f in [fft, ifft]:
|
||||
try:
|
||||
f(x)
|
||||
raise AssertionError("Type {0} not supported but does not fail" %
|
||||
np.longcomplex)
|
||||
except ValueError:
|
||||
pass
|
||||
|
||||
def test_real(self):
|
||||
if np.dtype(np.longdouble).itemsize == np.dtype(np.double).itemsize:
|
||||
# longdouble == double; so fft is supported
|
||||
return
|
||||
|
||||
x = np.random.randn(10).astype(np.longcomplex)
|
||||
|
||||
for f in [fft, ifft]:
|
||||
try:
|
||||
f(x)
|
||||
raise AssertionError("Type %r not supported but does not fail" %
|
||||
np.longcomplex)
|
||||
except ValueError:
|
||||
pass
|
||||
|
||||
|
||||
class FakeArray(object):
|
||||
def __init__(self, data):
|
||||
self._data = data
|
||||
self.__array_interface__ = data.__array_interface__
|
||||
|
||||
|
||||
class FakeArray2(object):
|
||||
def __init__(self, data):
|
||||
self._data = data
|
||||
|
||||
def __array__(self):
|
||||
return self._data
|
||||
|
||||
|
||||
class TestOverwrite(object):
|
||||
"""Check input overwrite behavior of the FFT functions."""
|
||||
|
||||
real_dtypes = [np.float32, np.float64]
|
||||
dtypes = real_dtypes + [np.complex64, np.complex128]
|
||||
fftsizes = [8, 16, 32]
|
||||
|
||||
def _check(self, x, routine, fftsize, axis, overwrite_x, should_overwrite):
|
||||
x2 = x.copy()
|
||||
for fake in [lambda x: x, FakeArray, FakeArray2]:
|
||||
routine(fake(x2), fftsize, axis, overwrite_x=overwrite_x)
|
||||
|
||||
sig = "%s(%s%r, %r, axis=%r, overwrite_x=%r)" % (
|
||||
routine.__name__, x.dtype, x.shape, fftsize, axis, overwrite_x)
|
||||
if not should_overwrite:
|
||||
assert_equal(x2, x, err_msg="spurious overwrite in %s" % sig)
|
||||
|
||||
def _check_1d(self, routine, dtype, shape, axis, overwritable_dtypes,
|
||||
fftsize, overwrite_x):
|
||||
np.random.seed(1234)
|
||||
if np.issubdtype(dtype, np.complexfloating):
|
||||
data = np.random.randn(*shape) + 1j*np.random.randn(*shape)
|
||||
else:
|
||||
data = np.random.randn(*shape)
|
||||
data = data.astype(dtype)
|
||||
|
||||
should_overwrite = (overwrite_x
|
||||
and dtype in overwritable_dtypes
|
||||
and fftsize <= shape[axis]
|
||||
and (len(shape) == 1 or
|
||||
(axis % len(shape) == len(shape)-1
|
||||
and fftsize == shape[axis])))
|
||||
self._check(data, routine, fftsize, axis,
|
||||
overwrite_x=overwrite_x,
|
||||
should_overwrite=should_overwrite)
|
||||
|
||||
@pytest.mark.parametrize('dtype', dtypes)
|
||||
@pytest.mark.parametrize('fftsize', fftsizes)
|
||||
@pytest.mark.parametrize('overwrite_x', [True, False])
|
||||
@pytest.mark.parametrize('shape,axes', [((16,), -1),
|
||||
((16, 2), 0),
|
||||
((2, 16), 1)])
|
||||
def test_fft_ifft(self, dtype, fftsize, overwrite_x, shape, axes):
|
||||
overwritable = (np.complex128, np.complex64)
|
||||
self._check_1d(fft, dtype, shape, axes, overwritable,
|
||||
fftsize, overwrite_x)
|
||||
self._check_1d(ifft, dtype, shape, axes, overwritable,
|
||||
fftsize, overwrite_x)
|
||||
|
||||
@pytest.mark.parametrize('dtype', real_dtypes)
|
||||
@pytest.mark.parametrize('fftsize', fftsizes)
|
||||
@pytest.mark.parametrize('overwrite_x', [True, False])
|
||||
@pytest.mark.parametrize('shape,axes', [((16,), -1),
|
||||
((16, 2), 0),
|
||||
((2, 16), 1)])
|
||||
def test_rfft_irfft(self, dtype, fftsize, overwrite_x, shape, axes):
|
||||
overwritable = self.real_dtypes
|
||||
self._check_1d(irfft, dtype, shape, axes, overwritable,
|
||||
fftsize, overwrite_x)
|
||||
self._check_1d(rfft, dtype, shape, axes, overwritable,
|
||||
fftsize, overwrite_x)
|
||||
|
||||
def _check_nd_one(self, routine, dtype, shape, axes, overwritable_dtypes,
|
||||
overwrite_x):
|
||||
np.random.seed(1234)
|
||||
if np.issubdtype(dtype, np.complexfloating):
|
||||
data = np.random.randn(*shape) + 1j*np.random.randn(*shape)
|
||||
else:
|
||||
data = np.random.randn(*shape)
|
||||
data = data.astype(dtype)
|
||||
|
||||
def fftshape_iter(shp):
|
||||
if len(shp) <= 0:
|
||||
yield ()
|
||||
else:
|
||||
for j in (shp[0]//2, shp[0], shp[0]*2):
|
||||
for rest in fftshape_iter(shp[1:]):
|
||||
yield (j,) + rest
|
||||
|
||||
if axes is None:
|
||||
part_shape = shape
|
||||
else:
|
||||
part_shape = tuple(np.take(shape, axes))
|
||||
|
||||
for fftshape in fftshape_iter(part_shape):
|
||||
should_overwrite = (overwrite_x
|
||||
and data.ndim == 1
|
||||
and np.all([x < y for x, y in zip(fftshape,
|
||||
part_shape)])
|
||||
and dtype in overwritable_dtypes)
|
||||
self._check(data, routine, fftshape, axes,
|
||||
overwrite_x=overwrite_x,
|
||||
should_overwrite=should_overwrite)
|
||||
if data.ndim > 1:
|
||||
# check fortran order: it never overwrites
|
||||
self._check(data.T, routine, fftshape, axes,
|
||||
overwrite_x=overwrite_x,
|
||||
should_overwrite=False)
|
||||
|
||||
@pytest.mark.parametrize('dtype', dtypes)
|
||||
@pytest.mark.parametrize('overwrite_x', [True, False])
|
||||
@pytest.mark.parametrize('shape,axes', [((16,), None),
|
||||
((16,), (0,)),
|
||||
((16, 2), (0,)),
|
||||
((2, 16), (1,)),
|
||||
((8, 16), None),
|
||||
((8, 16), (0, 1)),
|
||||
((8, 16, 2), (0, 1)),
|
||||
((8, 16, 2), (1, 2)),
|
||||
((8, 16, 2), (0,)),
|
||||
((8, 16, 2), (1,)),
|
||||
((8, 16, 2), (2,)),
|
||||
((8, 16, 2), None),
|
||||
((8, 16, 2), (0, 1, 2))])
|
||||
def test_fftn_ifftn(self, dtype, overwrite_x, shape, axes):
|
||||
overwritable = (np.complex128, np.complex64)
|
||||
self._check_nd_one(fftn, dtype, shape, axes, overwritable,
|
||||
overwrite_x)
|
||||
self._check_nd_one(ifftn, dtype, shape, axes, overwritable,
|
||||
overwrite_x)
|
||||
@@ -0,0 +1,416 @@
|
||||
# Created by Pearu Peterson, September 2002
|
||||
|
||||
from __future__ import division, print_function, absolute_import
|
||||
|
||||
__usage__ = """
|
||||
Build fftpack:
|
||||
python setup_fftpack.py build
|
||||
Run tests if scipy is installed:
|
||||
python -c 'import scipy;scipy.fftpack.test(<level>)'
|
||||
Run tests if fftpack is not installed:
|
||||
python tests/test_helper.py [<level>]
|
||||
"""
|
||||
|
||||
from pytest import raises as assert_raises
|
||||
from numpy.testing import assert_array_almost_equal, assert_equal, assert_
|
||||
from scipy.fftpack import fftshift,ifftshift,fftfreq,rfftfreq
|
||||
from scipy.fftpack.helper import (next_fast_len,
|
||||
_init_nd_shape_and_axes,
|
||||
_init_nd_shape_and_axes_sorted)
|
||||
|
||||
from numpy import pi, random
|
||||
import numpy as np
|
||||
|
||||
class TestFFTShift(object):
|
||||
|
||||
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.random((n,))
|
||||
assert_array_almost_equal(ifftshift(fftshift(x)),x)
|
||||
|
||||
|
||||
class TestFFTFreq(object):
|
||||
|
||||
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)
|
||||
|
||||
|
||||
class TestRFFTFreq(object):
|
||||
|
||||
def test_definition(self):
|
||||
x = [0,1,1,2,2,3,3,4,4]
|
||||
assert_array_almost_equal(9*rfftfreq(9),x)
|
||||
assert_array_almost_equal(9*pi*rfftfreq(9,pi),x)
|
||||
x = [0,1,1,2,2,3,3,4,4,5]
|
||||
assert_array_almost_equal(10*rfftfreq(10),x)
|
||||
assert_array_almost_equal(10*pi*rfftfreq(10,pi),x)
|
||||
|
||||
|
||||
class TestNextOptLen(object):
|
||||
|
||||
def test_next_opt_len(self):
|
||||
random.seed(1234)
|
||||
|
||||
def nums():
|
||||
for j in range(1, 1000):
|
||||
yield j
|
||||
yield 2**5 * 3**5 * 4**5 + 1
|
||||
|
||||
for n in nums():
|
||||
m = next_fast_len(n)
|
||||
msg = "n=%d, m=%d" % (n, m)
|
||||
|
||||
assert_(m >= n, msg)
|
||||
|
||||
# check regularity
|
||||
k = m
|
||||
for d in [2, 3, 5]:
|
||||
while True:
|
||||
a, b = divmod(k, d)
|
||||
if b == 0:
|
||||
k = a
|
||||
else:
|
||||
break
|
||||
assert_equal(k, 1, err_msg=msg)
|
||||
|
||||
def test_np_integers(self):
|
||||
ITYPES = [np.int16, np.int32, np.int64, np.uint16, np.uint32, np.uint64]
|
||||
for ityp in ITYPES:
|
||||
x = ityp(12345)
|
||||
testN = next_fast_len(x)
|
||||
assert_equal(testN, next_fast_len(int(x)))
|
||||
|
||||
def test_next_opt_len_strict(self):
|
||||
hams = {
|
||||
1: 1, 2: 2, 3: 3, 4: 4, 5: 5, 6: 6, 7: 8, 8: 8, 14: 15, 15: 15,
|
||||
16: 16, 17: 18, 1021: 1024, 1536: 1536, 51200000: 51200000,
|
||||
510183360: 510183360, 510183360 + 1: 512000000,
|
||||
511000000: 512000000,
|
||||
854296875: 854296875, 854296875 + 1: 859963392,
|
||||
196608000000: 196608000000, 196608000000 + 1: 196830000000,
|
||||
8789062500000: 8789062500000, 8789062500000 + 1: 8796093022208,
|
||||
206391214080000: 206391214080000,
|
||||
206391214080000 + 1: 206624260800000,
|
||||
470184984576000: 470184984576000,
|
||||
470184984576000 + 1: 470715894135000,
|
||||
7222041363087360: 7222041363087360,
|
||||
7222041363087360 + 1: 7230196133913600,
|
||||
# power of 5 5**23
|
||||
11920928955078125: 11920928955078125,
|
||||
11920928955078125 - 1: 11920928955078125,
|
||||
# power of 3 3**34
|
||||
16677181699666569: 16677181699666569,
|
||||
16677181699666569 - 1: 16677181699666569,
|
||||
# power of 2 2**54
|
||||
18014398509481984: 18014398509481984,
|
||||
18014398509481984 - 1: 18014398509481984,
|
||||
# above this, int(ceil(n)) == int(ceil(n+1))
|
||||
19200000000000000: 19200000000000000,
|
||||
19200000000000000 + 1: 19221679687500000,
|
||||
288230376151711744: 288230376151711744,
|
||||
288230376151711744 + 1: 288325195312500000,
|
||||
288325195312500000 - 1: 288325195312500000,
|
||||
288325195312500000: 288325195312500000,
|
||||
288325195312500000 + 1: 288555831593533440,
|
||||
# power of 3 3**83
|
||||
3990838394187339929534246675572349035227 - 1:
|
||||
3990838394187339929534246675572349035227,
|
||||
3990838394187339929534246675572349035227:
|
||||
3990838394187339929534246675572349035227,
|
||||
# power of 2 2**135
|
||||
43556142965880123323311949751266331066368 - 1:
|
||||
43556142965880123323311949751266331066368,
|
||||
43556142965880123323311949751266331066368:
|
||||
43556142965880123323311949751266331066368,
|
||||
# power of 5 5**57
|
||||
6938893903907228377647697925567626953125 - 1:
|
||||
6938893903907228377647697925567626953125,
|
||||
6938893903907228377647697925567626953125:
|
||||
6938893903907228377647697925567626953125,
|
||||
# http://www.drdobbs.com/228700538
|
||||
# 2**96 * 3**1 * 5**13
|
||||
290142196707511001929482240000000000000 - 1:
|
||||
290142196707511001929482240000000000000,
|
||||
290142196707511001929482240000000000000:
|
||||
290142196707511001929482240000000000000,
|
||||
290142196707511001929482240000000000000 + 1:
|
||||
290237644800000000000000000000000000000,
|
||||
# 2**36 * 3**69 * 5**7
|
||||
4479571262811807241115438439905203543080960000000 - 1:
|
||||
4479571262811807241115438439905203543080960000000,
|
||||
4479571262811807241115438439905203543080960000000:
|
||||
4479571262811807241115438439905203543080960000000,
|
||||
4479571262811807241115438439905203543080960000000 + 1:
|
||||
4480327901140333639941336854183943340032000000000,
|
||||
# 2**37 * 3**44 * 5**42
|
||||
30774090693237851027531250000000000000000000000000000000000000 - 1:
|
||||
30774090693237851027531250000000000000000000000000000000000000,
|
||||
30774090693237851027531250000000000000000000000000000000000000:
|
||||
30774090693237851027531250000000000000000000000000000000000000,
|
||||
30774090693237851027531250000000000000000000000000000000000000 + 1:
|
||||
30778180617309082445871527002041377406962596539492679680000000,
|
||||
}
|
||||
for x, y in hams.items():
|
||||
assert_equal(next_fast_len(x), y)
|
||||
|
||||
|
||||
class Test_init_nd_shape_and_axes(object):
|
||||
|
||||
def test_py_0d_defaults(self):
|
||||
x = 4
|
||||
shape = None
|
||||
axes = None
|
||||
|
||||
shape_expected = np.array([])
|
||||
axes_expected = np.array([])
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes_sorted(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
def test_np_0d_defaults(self):
|
||||
x = np.array(7.)
|
||||
shape = None
|
||||
axes = None
|
||||
|
||||
shape_expected = np.array([])
|
||||
axes_expected = np.array([])
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes_sorted(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
def test_py_1d_defaults(self):
|
||||
x = [1, 2, 3]
|
||||
shape = None
|
||||
axes = None
|
||||
|
||||
shape_expected = np.array([3])
|
||||
axes_expected = np.array([0])
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes_sorted(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
def test_np_1d_defaults(self):
|
||||
x = np.arange(0, 1, .1)
|
||||
shape = None
|
||||
axes = None
|
||||
|
||||
shape_expected = np.array([10])
|
||||
axes_expected = np.array([0])
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes_sorted(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
def test_py_2d_defaults(self):
|
||||
x = [[1, 2, 3, 4],
|
||||
[5, 6, 7, 8]]
|
||||
shape = None
|
||||
axes = None
|
||||
|
||||
shape_expected = np.array([2, 4])
|
||||
axes_expected = np.array([0, 1])
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes_sorted(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
def test_np_2d_defaults(self):
|
||||
x = np.arange(0, 1, .1).reshape(5, 2)
|
||||
shape = None
|
||||
axes = None
|
||||
|
||||
shape_expected = np.array([5, 2])
|
||||
axes_expected = np.array([0, 1])
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes_sorted(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
def test_np_5d_defaults(self):
|
||||
x = np.zeros([6, 2, 5, 3, 4])
|
||||
shape = None
|
||||
axes = None
|
||||
|
||||
shape_expected = np.array([6, 2, 5, 3, 4])
|
||||
axes_expected = np.array([0, 1, 2, 3, 4])
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes_sorted(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
def test_np_5d_set_shape(self):
|
||||
x = np.zeros([6, 2, 5, 3, 4])
|
||||
shape = [10, -1, -1, 1, 4]
|
||||
axes = None
|
||||
|
||||
shape_expected = np.array([10, 2, 5, 1, 4])
|
||||
axes_expected = np.array([0, 1, 2, 3, 4])
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes_sorted(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
def test_np_5d_set_axes(self):
|
||||
x = np.zeros([6, 2, 5, 3, 4])
|
||||
shape = None
|
||||
axes = [4, 1, 2]
|
||||
|
||||
shape_expected = np.array([4, 2, 5])
|
||||
axes_expected = np.array([4, 1, 2])
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
def test_np_5d_set_axes_sorted(self):
|
||||
x = np.zeros([6, 2, 5, 3, 4])
|
||||
shape = None
|
||||
axes = [4, 1, 2]
|
||||
|
||||
shape_expected = np.array([2, 5, 4])
|
||||
axes_expected = np.array([1, 2, 4])
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes_sorted(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
def test_np_5d_set_shape_axes(self):
|
||||
x = np.zeros([6, 2, 5, 3, 4])
|
||||
shape = [10, -1, 2]
|
||||
axes = [1, 0, 3]
|
||||
|
||||
shape_expected = np.array([10, 6, 2])
|
||||
axes_expected = np.array([1, 0, 3])
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
def test_np_5d_set_shape_axes_sorted(self):
|
||||
x = np.zeros([6, 2, 5, 3, 4])
|
||||
shape = [10, -1, 2]
|
||||
axes = [1, 0, 3]
|
||||
|
||||
shape_expected = np.array([6, 10, 2])
|
||||
axes_expected = np.array([0, 1, 3])
|
||||
|
||||
shape_res, axes_res = _init_nd_shape_and_axes_sorted(x, shape, axes)
|
||||
|
||||
assert_equal(shape_res, shape_expected)
|
||||
assert_equal(axes_res, axes_expected)
|
||||
|
||||
def test_errors(self):
|
||||
with assert_raises(ValueError,
|
||||
match="when given, axes values must be a scalar"
|
||||
" or vector"):
|
||||
_init_nd_shape_and_axes([0], shape=None, axes=[[1, 2], [3, 4]])
|
||||
|
||||
with assert_raises(ValueError,
|
||||
match="when given, axes values must be integers"):
|
||||
_init_nd_shape_and_axes([0], shape=None, axes=[1., 2., 3., 4.])
|
||||
|
||||
with assert_raises(ValueError,
|
||||
match="axes exceeds dimensionality of input"):
|
||||
_init_nd_shape_and_axes([0], shape=None, axes=[1])
|
||||
|
||||
with assert_raises(ValueError,
|
||||
match="axes exceeds dimensionality of input"):
|
||||
_init_nd_shape_and_axes([0], shape=None, axes=[-2])
|
||||
|
||||
with assert_raises(ValueError,
|
||||
match="all axes must be unique"):
|
||||
_init_nd_shape_and_axes([0], shape=None, axes=[0, 0])
|
||||
|
||||
with assert_raises(ValueError,
|
||||
match="when given, shape values must be a scalar "
|
||||
"or vector"):
|
||||
_init_nd_shape_and_axes([0], shape=[[1, 2], [3, 4]], axes=None)
|
||||
|
||||
with assert_raises(ValueError,
|
||||
match="when given, shape values must be integers"):
|
||||
_init_nd_shape_and_axes([0], shape=[1., 2., 3., 4.], axes=None)
|
||||
|
||||
with assert_raises(ValueError,
|
||||
match="when given, axes and shape arguments"
|
||||
" have to be of the same length"):
|
||||
_init_nd_shape_and_axes(np.zeros([1, 1, 1, 1]),
|
||||
shape=[1, 2, 3], axes=[1])
|
||||
|
||||
with assert_raises(ValueError,
|
||||
match="invalid number of data points"
|
||||
r" \(\[0\]\) specified"):
|
||||
_init_nd_shape_and_axes([0], shape=[0], axes=None)
|
||||
|
||||
with assert_raises(ValueError,
|
||||
match="invalid number of data points"
|
||||
r" \(\[-2\]\) specified"):
|
||||
_init_nd_shape_and_axes([0], shape=-2, axes=None)
|
||||
@@ -0,0 +1,33 @@
|
||||
"""Test possibility of patching fftpack with pyfftw.
|
||||
|
||||
No module source outside of scipy.fftpack should contain an import of
|
||||
the form `from scipy.fftpack import ...`, so that a simple replacement
|
||||
of scipy.fftpack by the corresponding fftw interface completely swaps
|
||||
the two FFT implementations.
|
||||
|
||||
Because this simply inspects source files, we only need to run the test
|
||||
on one version of Python.
|
||||
"""
|
||||
|
||||
|
||||
import sys
|
||||
if sys.version_info >= (3, 4):
|
||||
from pathlib import Path
|
||||
import re
|
||||
import tokenize
|
||||
from numpy.testing import assert_
|
||||
import scipy
|
||||
|
||||
class TestFFTPackImport(object):
|
||||
def test_fftpack_import(self):
|
||||
base = Path(scipy.__file__).parent
|
||||
regexp = r"\s*from.+\.fftpack import .*\n"
|
||||
for path in base.rglob("*.py"):
|
||||
if base / "fftpack" in path.parents:
|
||||
continue
|
||||
# use tokenize to auto-detect encoding on systems where no
|
||||
# default encoding is defined (e.g. LANG='C')
|
||||
with tokenize.open(str(path)) as file:
|
||||
assert_(all(not re.fullmatch(regexp, line)
|
||||
for line in file),
|
||||
"{0} contains an import from fftpack".format(path))
|
||||
@@ -0,0 +1,382 @@
|
||||
# Created by Pearu Peterson, September 2002
|
||||
|
||||
from __future__ import division, print_function, absolute_import
|
||||
|
||||
__usage__ = """
|
||||
Build fftpack:
|
||||
python setup_fftpack.py build
|
||||
Run tests if scipy is installed:
|
||||
python -c 'import scipy;scipy.fftpack.test(<level>)'
|
||||
Run tests if fftpack is not installed:
|
||||
python tests/test_pseudo_diffs.py [<level>]
|
||||
"""
|
||||
|
||||
from numpy.testing import (assert_equal, assert_almost_equal,
|
||||
assert_array_almost_equal)
|
||||
from scipy.fftpack import (diff, fft, ifft, tilbert, itilbert, hilbert,
|
||||
ihilbert, shift, fftfreq, cs_diff, sc_diff,
|
||||
ss_diff, cc_diff)
|
||||
|
||||
import numpy as np
|
||||
from numpy import arange, sin, cos, pi, exp, tanh, sum, sign
|
||||
from numpy.random import random
|
||||
|
||||
|
||||
def direct_diff(x,k=1,period=None):
|
||||
fx = fft(x)
|
||||
n = len(fx)
|
||||
if period is None:
|
||||
period = 2*pi
|
||||
w = fftfreq(n)*2j*pi/period*n
|
||||
if k < 0:
|
||||
w = 1 / w**k
|
||||
w[0] = 0.0
|
||||
else:
|
||||
w = w**k
|
||||
if n > 2000:
|
||||
w[250:n-250] = 0.0
|
||||
return ifft(w*fx).real
|
||||
|
||||
|
||||
def direct_tilbert(x,h=1,period=None):
|
||||
fx = fft(x)
|
||||
n = len(fx)
|
||||
if period is None:
|
||||
period = 2*pi
|
||||
w = fftfreq(n)*h*2*pi/period*n
|
||||
w[0] = 1
|
||||
w = 1j/tanh(w)
|
||||
w[0] = 0j
|
||||
return ifft(w*fx)
|
||||
|
||||
|
||||
def direct_itilbert(x,h=1,period=None):
|
||||
fx = fft(x)
|
||||
n = len(fx)
|
||||
if period is None:
|
||||
period = 2*pi
|
||||
w = fftfreq(n)*h*2*pi/period*n
|
||||
w = -1j*tanh(w)
|
||||
return ifft(w*fx)
|
||||
|
||||
|
||||
def direct_hilbert(x):
|
||||
fx = fft(x)
|
||||
n = len(fx)
|
||||
w = fftfreq(n)*n
|
||||
w = 1j*sign(w)
|
||||
return ifft(w*fx)
|
||||
|
||||
|
||||
def direct_ihilbert(x):
|
||||
return -direct_hilbert(x)
|
||||
|
||||
|
||||
def direct_shift(x,a,period=None):
|
||||
n = len(x)
|
||||
if period is None:
|
||||
k = fftfreq(n)*1j*n
|
||||
else:
|
||||
k = fftfreq(n)*2j*pi/period*n
|
||||
return ifft(fft(x)*exp(k*a)).real
|
||||
|
||||
|
||||
class TestDiff(object):
|
||||
|
||||
def test_definition(self):
|
||||
for n in [16,17,64,127,32]:
|
||||
x = arange(n)*2*pi/n
|
||||
assert_array_almost_equal(diff(sin(x)),direct_diff(sin(x)))
|
||||
assert_array_almost_equal(diff(sin(x),2),direct_diff(sin(x),2))
|
||||
assert_array_almost_equal(diff(sin(x),3),direct_diff(sin(x),3))
|
||||
assert_array_almost_equal(diff(sin(x),4),direct_diff(sin(x),4))
|
||||
assert_array_almost_equal(diff(sin(x),5),direct_diff(sin(x),5))
|
||||
assert_array_almost_equal(diff(sin(2*x),3),direct_diff(sin(2*x),3))
|
||||
assert_array_almost_equal(diff(sin(2*x),4),direct_diff(sin(2*x),4))
|
||||
assert_array_almost_equal(diff(cos(x)),direct_diff(cos(x)))
|
||||
assert_array_almost_equal(diff(cos(x),2),direct_diff(cos(x),2))
|
||||
assert_array_almost_equal(diff(cos(x),3),direct_diff(cos(x),3))
|
||||
assert_array_almost_equal(diff(cos(x),4),direct_diff(cos(x),4))
|
||||
assert_array_almost_equal(diff(cos(2*x)),direct_diff(cos(2*x)))
|
||||
assert_array_almost_equal(diff(sin(x*n/8)),direct_diff(sin(x*n/8)))
|
||||
assert_array_almost_equal(diff(cos(x*n/8)),direct_diff(cos(x*n/8)))
|
||||
for k in range(5):
|
||||
assert_array_almost_equal(diff(sin(4*x),k),direct_diff(sin(4*x),k))
|
||||
assert_array_almost_equal(diff(cos(4*x),k),direct_diff(cos(4*x),k))
|
||||
|
||||
def test_period(self):
|
||||
for n in [17,64]:
|
||||
x = arange(n)/float(n)
|
||||
assert_array_almost_equal(diff(sin(2*pi*x),period=1),
|
||||
2*pi*cos(2*pi*x))
|
||||
assert_array_almost_equal(diff(sin(2*pi*x),3,period=1),
|
||||
-(2*pi)**3*cos(2*pi*x))
|
||||
|
||||
def test_sin(self):
|
||||
for n in [32,64,77]:
|
||||
x = arange(n)*2*pi/n
|
||||
assert_array_almost_equal(diff(sin(x)),cos(x))
|
||||
assert_array_almost_equal(diff(cos(x)),-sin(x))
|
||||
assert_array_almost_equal(diff(sin(x),2),-sin(x))
|
||||
assert_array_almost_equal(diff(sin(x),4),sin(x))
|
||||
assert_array_almost_equal(diff(sin(4*x)),4*cos(4*x))
|
||||
assert_array_almost_equal(diff(sin(sin(x))),cos(x)*cos(sin(x)))
|
||||
|
||||
def test_expr(self):
|
||||
for n in [64,77,100,128,256,512,1024,2048,4096,8192][:5]:
|
||||
x = arange(n)*2*pi/n
|
||||
f = sin(x)*cos(4*x)+exp(sin(3*x))
|
||||
df = cos(x)*cos(4*x)-4*sin(x)*sin(4*x)+3*cos(3*x)*exp(sin(3*x))
|
||||
ddf = -17*sin(x)*cos(4*x)-8*cos(x)*sin(4*x)\
|
||||
- 9*sin(3*x)*exp(sin(3*x))+9*cos(3*x)**2*exp(sin(3*x))
|
||||
d1 = diff(f)
|
||||
assert_array_almost_equal(d1,df)
|
||||
assert_array_almost_equal(diff(df),ddf)
|
||||
assert_array_almost_equal(diff(f,2),ddf)
|
||||
assert_array_almost_equal(diff(ddf,-1),df)
|
||||
|
||||
def test_expr_large(self):
|
||||
for n in [2048,4096]:
|
||||
x = arange(n)*2*pi/n
|
||||
f = sin(x)*cos(4*x)+exp(sin(3*x))
|
||||
df = cos(x)*cos(4*x)-4*sin(x)*sin(4*x)+3*cos(3*x)*exp(sin(3*x))
|
||||
ddf = -17*sin(x)*cos(4*x)-8*cos(x)*sin(4*x)\
|
||||
- 9*sin(3*x)*exp(sin(3*x))+9*cos(3*x)**2*exp(sin(3*x))
|
||||
assert_array_almost_equal(diff(f),df)
|
||||
assert_array_almost_equal(diff(df),ddf)
|
||||
assert_array_almost_equal(diff(ddf,-1),df)
|
||||
assert_array_almost_equal(diff(f,2),ddf)
|
||||
|
||||
def test_int(self):
|
||||
n = 64
|
||||
x = arange(n)*2*pi/n
|
||||
assert_array_almost_equal(diff(sin(x),-1),-cos(x))
|
||||
assert_array_almost_equal(diff(sin(x),-2),-sin(x))
|
||||
assert_array_almost_equal(diff(sin(x),-4),sin(x))
|
||||
assert_array_almost_equal(diff(2*cos(2*x),-1),sin(2*x))
|
||||
|
||||
def test_random_even(self):
|
||||
for k in [0,2,4,6]:
|
||||
for n in [60,32,64,56,55]:
|
||||
f = random((n,))
|
||||
af = sum(f,axis=0)/n
|
||||
f = f-af
|
||||
# zeroing Nyquist mode:
|
||||
f = diff(diff(f,1),-1)
|
||||
assert_almost_equal(sum(f,axis=0),0.0)
|
||||
assert_array_almost_equal(diff(diff(f,k),-k),f)
|
||||
assert_array_almost_equal(diff(diff(f,-k),k),f)
|
||||
|
||||
def test_random_odd(self):
|
||||
for k in [0,1,2,3,4,5,6]:
|
||||
for n in [33,65,55]:
|
||||
f = random((n,))
|
||||
af = sum(f,axis=0)/n
|
||||
f = f-af
|
||||
assert_almost_equal(sum(f,axis=0),0.0)
|
||||
assert_array_almost_equal(diff(diff(f,k),-k),f)
|
||||
assert_array_almost_equal(diff(diff(f,-k),k),f)
|
||||
|
||||
def test_zero_nyquist(self):
|
||||
for k in [0,1,2,3,4,5,6]:
|
||||
for n in [32,33,64,56,55]:
|
||||
f = random((n,))
|
||||
af = sum(f,axis=0)/n
|
||||
f = f-af
|
||||
# zeroing Nyquist mode:
|
||||
f = diff(diff(f,1),-1)
|
||||
assert_almost_equal(sum(f,axis=0),0.0)
|
||||
assert_array_almost_equal(diff(diff(f,k),-k),f)
|
||||
assert_array_almost_equal(diff(diff(f,-k),k),f)
|
||||
|
||||
|
||||
class TestTilbert(object):
|
||||
|
||||
def test_definition(self):
|
||||
for h in [0.1,0.5,1,5.5,10]:
|
||||
for n in [16,17,64,127]:
|
||||
x = arange(n)*2*pi/n
|
||||
y = tilbert(sin(x),h)
|
||||
y1 = direct_tilbert(sin(x),h)
|
||||
assert_array_almost_equal(y,y1)
|
||||
assert_array_almost_equal(tilbert(sin(x),h),
|
||||
direct_tilbert(sin(x),h))
|
||||
assert_array_almost_equal(tilbert(sin(2*x),h),
|
||||
direct_tilbert(sin(2*x),h))
|
||||
|
||||
def test_random_even(self):
|
||||
for h in [0.1,0.5,1,5.5,10]:
|
||||
for n in [32,64,56]:
|
||||
f = random((n,))
|
||||
af = sum(f,axis=0)/n
|
||||
f = f-af
|
||||
assert_almost_equal(sum(f,axis=0),0.0)
|
||||
assert_array_almost_equal(direct_tilbert(direct_itilbert(f,h),h),f)
|
||||
|
||||
def test_random_odd(self):
|
||||
for h in [0.1,0.5,1,5.5,10]:
|
||||
for n in [33,65,55]:
|
||||
f = random((n,))
|
||||
af = sum(f,axis=0)/n
|
||||
f = f-af
|
||||
assert_almost_equal(sum(f,axis=0),0.0)
|
||||
assert_array_almost_equal(itilbert(tilbert(f,h),h),f)
|
||||
assert_array_almost_equal(tilbert(itilbert(f,h),h),f)
|
||||
|
||||
|
||||
class TestITilbert(object):
|
||||
|
||||
def test_definition(self):
|
||||
for h in [0.1,0.5,1,5.5,10]:
|
||||
for n in [16,17,64,127]:
|
||||
x = arange(n)*2*pi/n
|
||||
y = itilbert(sin(x),h)
|
||||
y1 = direct_itilbert(sin(x),h)
|
||||
assert_array_almost_equal(y,y1)
|
||||
assert_array_almost_equal(itilbert(sin(x),h),
|
||||
direct_itilbert(sin(x),h))
|
||||
assert_array_almost_equal(itilbert(sin(2*x),h),
|
||||
direct_itilbert(sin(2*x),h))
|
||||
|
||||
|
||||
class TestHilbert(object):
|
||||
|
||||
def test_definition(self):
|
||||
for n in [16,17,64,127]:
|
||||
x = arange(n)*2*pi/n
|
||||
y = hilbert(sin(x))
|
||||
y1 = direct_hilbert(sin(x))
|
||||
assert_array_almost_equal(y,y1)
|
||||
assert_array_almost_equal(hilbert(sin(2*x)),
|
||||
direct_hilbert(sin(2*x)))
|
||||
|
||||
def test_tilbert_relation(self):
|
||||
for n in [16,17,64,127]:
|
||||
x = arange(n)*2*pi/n
|
||||
f = sin(x)+cos(2*x)*sin(x)
|
||||
y = hilbert(f)
|
||||
y1 = direct_hilbert(f)
|
||||
assert_array_almost_equal(y,y1)
|
||||
y2 = tilbert(f,h=10)
|
||||
assert_array_almost_equal(y,y2)
|
||||
|
||||
def test_random_odd(self):
|
||||
for n in [33,65,55]:
|
||||
f = random((n,))
|
||||
af = sum(f,axis=0)/n
|
||||
f = f-af
|
||||
assert_almost_equal(sum(f,axis=0),0.0)
|
||||
assert_array_almost_equal(ihilbert(hilbert(f)),f)
|
||||
assert_array_almost_equal(hilbert(ihilbert(f)),f)
|
||||
|
||||
def test_random_even(self):
|
||||
for n in [32,64,56]:
|
||||
f = random((n,))
|
||||
af = sum(f,axis=0)/n
|
||||
f = f-af
|
||||
# zeroing Nyquist mode:
|
||||
f = diff(diff(f,1),-1)
|
||||
assert_almost_equal(sum(f,axis=0),0.0)
|
||||
assert_array_almost_equal(direct_hilbert(direct_ihilbert(f)),f)
|
||||
assert_array_almost_equal(hilbert(ihilbert(f)),f)
|
||||
|
||||
|
||||
class TestIHilbert(object):
|
||||
|
||||
def test_definition(self):
|
||||
for n in [16,17,64,127]:
|
||||
x = arange(n)*2*pi/n
|
||||
y = ihilbert(sin(x))
|
||||
y1 = direct_ihilbert(sin(x))
|
||||
assert_array_almost_equal(y,y1)
|
||||
assert_array_almost_equal(ihilbert(sin(2*x)),
|
||||
direct_ihilbert(sin(2*x)))
|
||||
|
||||
def test_itilbert_relation(self):
|
||||
for n in [16,17,64,127]:
|
||||
x = arange(n)*2*pi/n
|
||||
f = sin(x)+cos(2*x)*sin(x)
|
||||
y = ihilbert(f)
|
||||
y1 = direct_ihilbert(f)
|
||||
assert_array_almost_equal(y,y1)
|
||||
y2 = itilbert(f,h=10)
|
||||
assert_array_almost_equal(y,y2)
|
||||
|
||||
|
||||
class TestShift(object):
|
||||
|
||||
def test_definition(self):
|
||||
for n in [18,17,64,127,32,2048,256]:
|
||||
x = arange(n)*2*pi/n
|
||||
for a in [0.1,3]:
|
||||
assert_array_almost_equal(shift(sin(x),a),direct_shift(sin(x),a))
|
||||
assert_array_almost_equal(shift(sin(x),a),sin(x+a))
|
||||
assert_array_almost_equal(shift(cos(x),a),cos(x+a))
|
||||
assert_array_almost_equal(shift(cos(2*x)+sin(x),a),
|
||||
cos(2*(x+a))+sin(x+a))
|
||||
assert_array_almost_equal(shift(exp(sin(x)),a),exp(sin(x+a)))
|
||||
assert_array_almost_equal(shift(sin(x),2*pi),sin(x))
|
||||
assert_array_almost_equal(shift(sin(x),pi),-sin(x))
|
||||
assert_array_almost_equal(shift(sin(x),pi/2),cos(x))
|
||||
|
||||
|
||||
class TestOverwrite(object):
|
||||
"""Check input overwrite behavior """
|
||||
|
||||
real_dtypes = [np.float32, np.float64]
|
||||
dtypes = real_dtypes + [np.complex64, np.complex128]
|
||||
|
||||
def _check(self, x, routine, *args, **kwargs):
|
||||
x2 = x.copy()
|
||||
routine(x2, *args, **kwargs)
|
||||
sig = routine.__name__
|
||||
if args:
|
||||
sig += repr(args)
|
||||
if kwargs:
|
||||
sig += repr(kwargs)
|
||||
assert_equal(x2, x, err_msg="spurious overwrite in %s" % sig)
|
||||
|
||||
def _check_1d(self, routine, dtype, shape, *args, **kwargs):
|
||||
np.random.seed(1234)
|
||||
if np.issubdtype(dtype, np.complexfloating):
|
||||
data = np.random.randn(*shape) + 1j*np.random.randn(*shape)
|
||||
else:
|
||||
data = np.random.randn(*shape)
|
||||
data = data.astype(dtype)
|
||||
self._check(data, routine, *args, **kwargs)
|
||||
|
||||
def test_diff(self):
|
||||
for dtype in self.dtypes:
|
||||
self._check_1d(diff, dtype, (16,))
|
||||
|
||||
def test_tilbert(self):
|
||||
for dtype in self.dtypes:
|
||||
self._check_1d(tilbert, dtype, (16,), 1.6)
|
||||
|
||||
def test_itilbert(self):
|
||||
for dtype in self.dtypes:
|
||||
self._check_1d(itilbert, dtype, (16,), 1.6)
|
||||
|
||||
def test_hilbert(self):
|
||||
for dtype in self.dtypes:
|
||||
self._check_1d(hilbert, dtype, (16,))
|
||||
|
||||
def test_cs_diff(self):
|
||||
for dtype in self.dtypes:
|
||||
self._check_1d(cs_diff, dtype, (16,), 1.0, 4.0)
|
||||
|
||||
def test_sc_diff(self):
|
||||
for dtype in self.dtypes:
|
||||
self._check_1d(sc_diff, dtype, (16,), 1.0, 4.0)
|
||||
|
||||
def test_ss_diff(self):
|
||||
for dtype in self.dtypes:
|
||||
self._check_1d(ss_diff, dtype, (16,), 1.0, 4.0)
|
||||
|
||||
def test_cc_diff(self):
|
||||
for dtype in self.dtypes:
|
||||
self._check_1d(cc_diff, dtype, (16,), 1.0, 4.0)
|
||||
|
||||
def test_shift(self):
|
||||
for dtype in self.dtypes:
|
||||
self._check_1d(shift, dtype, (16,), 1.0)
|
||||
@@ -0,0 +1,829 @@
|
||||
from __future__ import division, print_function, absolute_import
|
||||
|
||||
from os.path import join, dirname
|
||||
|
||||
import numpy as np
|
||||
from numpy.testing import assert_array_almost_equal, assert_equal
|
||||
import pytest
|
||||
from pytest import raises as assert_raises
|
||||
|
||||
from scipy.fftpack.realtransforms import (
|
||||
dct, idct, dst, idst, dctn, idctn, dstn, idstn)
|
||||
|
||||
# Matlab reference data
|
||||
MDATA = np.load(join(dirname(__file__), 'test.npz'))
|
||||
X = [MDATA['x%d' % i] for i in range(8)]
|
||||
Y = [MDATA['y%d' % i] for i in range(8)]
|
||||
|
||||
# FFTW reference data: the data are organized as follows:
|
||||
# * SIZES is an array containing all available sizes
|
||||
# * for every type (1, 2, 3, 4) and every size, the array dct_type_size
|
||||
# contains the output of the DCT applied to the input np.linspace(0, size-1,
|
||||
# size)
|
||||
FFTWDATA_DOUBLE = np.load(join(dirname(__file__), 'fftw_double_ref.npz'))
|
||||
FFTWDATA_SINGLE = np.load(join(dirname(__file__), 'fftw_single_ref.npz'))
|
||||
FFTWDATA_SIZES = FFTWDATA_DOUBLE['sizes']
|
||||
|
||||
|
||||
def fftw_dct_ref(type, size, dt):
|
||||
x = np.linspace(0, size-1, size).astype(dt)
|
||||
dt = np.result_type(np.float32, dt)
|
||||
if dt == np.double:
|
||||
data = FFTWDATA_DOUBLE
|
||||
elif dt == np.float32:
|
||||
data = FFTWDATA_SINGLE
|
||||
else:
|
||||
raise ValueError()
|
||||
y = (data['dct_%d_%d' % (type, size)]).astype(dt)
|
||||
return x, y, dt
|
||||
|
||||
|
||||
def fftw_dst_ref(type, size, dt):
|
||||
x = np.linspace(0, size-1, size).astype(dt)
|
||||
dt = np.result_type(np.float32, dt)
|
||||
if dt == np.double:
|
||||
data = FFTWDATA_DOUBLE
|
||||
elif dt == np.float32:
|
||||
data = FFTWDATA_SINGLE
|
||||
else:
|
||||
raise ValueError()
|
||||
y = (data['dst_%d_%d' % (type, size)]).astype(dt)
|
||||
return x, y, dt
|
||||
|
||||
|
||||
def dct_2d_ref(x, **kwargs):
|
||||
"""Calculate reference values for testing dct2."""
|
||||
x = np.array(x, copy=True)
|
||||
for row in range(x.shape[0]):
|
||||
x[row, :] = dct(x[row, :], **kwargs)
|
||||
for col in range(x.shape[1]):
|
||||
x[:, col] = dct(x[:, col], **kwargs)
|
||||
return x
|
||||
|
||||
|
||||
def idct_2d_ref(x, **kwargs):
|
||||
"""Calculate reference values for testing idct2."""
|
||||
x = np.array(x, copy=True)
|
||||
for row in range(x.shape[0]):
|
||||
x[row, :] = idct(x[row, :], **kwargs)
|
||||
for col in range(x.shape[1]):
|
||||
x[:, col] = idct(x[:, col], **kwargs)
|
||||
return x
|
||||
|
||||
|
||||
def dst_2d_ref(x, **kwargs):
|
||||
"""Calculate reference values for testing dst2."""
|
||||
x = np.array(x, copy=True)
|
||||
for row in range(x.shape[0]):
|
||||
x[row, :] = dst(x[row, :], **kwargs)
|
||||
for col in range(x.shape[1]):
|
||||
x[:, col] = dst(x[:, col], **kwargs)
|
||||
return x
|
||||
|
||||
|
||||
def idst_2d_ref(x, **kwargs):
|
||||
"""Calculate reference values for testing idst2."""
|
||||
x = np.array(x, copy=True)
|
||||
for row in range(x.shape[0]):
|
||||
x[row, :] = idst(x[row, :], **kwargs)
|
||||
for col in range(x.shape[1]):
|
||||
x[:, col] = idst(x[:, col], **kwargs)
|
||||
return x
|
||||
|
||||
|
||||
def naive_dct1(x, norm=None):
|
||||
"""Calculate textbook definition version of DCT-I."""
|
||||
x = np.array(x, copy=True)
|
||||
N = len(x)
|
||||
M = N-1
|
||||
y = np.zeros(N)
|
||||
m0, m = 1, 2
|
||||
if norm == 'ortho':
|
||||
m0 = np.sqrt(1.0/M)
|
||||
m = np.sqrt(2.0/M)
|
||||
for k in range(N):
|
||||
for n in range(1, N-1):
|
||||
y[k] += m*x[n]*np.cos(np.pi*n*k/M)
|
||||
y[k] += m0 * x[0]
|
||||
y[k] += m0 * x[N-1] * (1 if k % 2 == 0 else -1)
|
||||
if norm == 'ortho':
|
||||
y[0] *= 1/np.sqrt(2)
|
||||
y[N-1] *= 1/np.sqrt(2)
|
||||
return y
|
||||
|
||||
|
||||
def naive_dst1(x, norm=None):
|
||||
"""Calculate textbook definition version of DST-I."""
|
||||
x = np.array(x, copy=True)
|
||||
N = len(x)
|
||||
M = N+1
|
||||
y = np.zeros(N)
|
||||
for k in range(N):
|
||||
for n in range(N):
|
||||
y[k] += 2*x[n]*np.sin(np.pi*(n+1.0)*(k+1.0)/M)
|
||||
if norm == 'ortho':
|
||||
y *= np.sqrt(0.5/M)
|
||||
return y
|
||||
|
||||
|
||||
def naive_dct4(x, norm=None):
|
||||
"""Calculate textbook definition version of DCT-IV."""
|
||||
x = np.array(x, copy=True)
|
||||
N = len(x)
|
||||
y = np.zeros(N)
|
||||
for k in range(N):
|
||||
for n in range(N):
|
||||
y[k] += x[n]*np.cos(np.pi*(n+0.5)*(k+0.5)/(N))
|
||||
if norm == 'ortho':
|
||||
y *= np.sqrt(2.0/N)
|
||||
else:
|
||||
y *= 2
|
||||
return y
|
||||
|
||||
|
||||
def naive_dst4(x, norm=None):
|
||||
"""Calculate textbook definition version of DST-IV."""
|
||||
x = np.array(x, copy=True)
|
||||
N = len(x)
|
||||
y = np.zeros(N)
|
||||
for k in range(N):
|
||||
for n in range(N):
|
||||
y[k] += x[n]*np.sin(np.pi*(n+0.5)*(k+0.5)/(N))
|
||||
if norm == 'ortho':
|
||||
y *= np.sqrt(2.0/N)
|
||||
else:
|
||||
y *= 2
|
||||
return y
|
||||
|
||||
|
||||
class TestComplex(object):
|
||||
def test_dct_complex64(self):
|
||||
y = dct(1j*np.arange(5, dtype=np.complex64))
|
||||
x = 1j*dct(np.arange(5))
|
||||
assert_array_almost_equal(x, y)
|
||||
|
||||
def test_dct_complex(self):
|
||||
y = dct(np.arange(5)*1j)
|
||||
x = 1j*dct(np.arange(5))
|
||||
assert_array_almost_equal(x, y)
|
||||
|
||||
def test_idct_complex(self):
|
||||
y = idct(np.arange(5)*1j)
|
||||
x = 1j*idct(np.arange(5))
|
||||
assert_array_almost_equal(x, y)
|
||||
|
||||
def test_dst_complex64(self):
|
||||
y = dst(np.arange(5, dtype=np.complex64)*1j)
|
||||
x = 1j*dst(np.arange(5))
|
||||
assert_array_almost_equal(x, y)
|
||||
|
||||
def test_dst_complex(self):
|
||||
y = dst(np.arange(5)*1j)
|
||||
x = 1j*dst(np.arange(5))
|
||||
assert_array_almost_equal(x, y)
|
||||
|
||||
def test_idst_complex(self):
|
||||
y = idst(np.arange(5)*1j)
|
||||
x = 1j*idst(np.arange(5))
|
||||
assert_array_almost_equal(x, y)
|
||||
|
||||
|
||||
class _TestDCTBase(object):
|
||||
def setup_method(self):
|
||||
self.rdt = None
|
||||
self.dec = 14
|
||||
self.type = None
|
||||
|
||||
def test_definition(self):
|
||||
for i in FFTWDATA_SIZES:
|
||||
x, yr, dt = fftw_dct_ref(self.type, i, self.rdt)
|
||||
y = dct(x, type=self.type)
|
||||
assert_equal(y.dtype, dt)
|
||||
# XXX: we divide by np.max(y) because the tests fail otherwise. We
|
||||
# should really use something like assert_array_approx_equal. The
|
||||
# difference is due to fftw using a better algorithm w.r.t error
|
||||
# propagation compared to the ones from fftpack.
|
||||
assert_array_almost_equal(y / np.max(y), yr / np.max(y), decimal=self.dec,
|
||||
err_msg="Size %d failed" % i)
|
||||
|
||||
def test_axis(self):
|
||||
nt = 2
|
||||
for i in [7, 8, 9, 16, 32, 64]:
|
||||
x = np.random.randn(nt, i)
|
||||
y = dct(x, type=self.type)
|
||||
for j in range(nt):
|
||||
assert_array_almost_equal(y[j], dct(x[j], type=self.type),
|
||||
decimal=self.dec)
|
||||
|
||||
x = x.T
|
||||
y = dct(x, axis=0, type=self.type)
|
||||
for j in range(nt):
|
||||
assert_array_almost_equal(y[:,j], dct(x[:,j], type=self.type),
|
||||
decimal=self.dec)
|
||||
|
||||
|
||||
class _TestDCTIBase(_TestDCTBase):
|
||||
def test_definition_ortho(self):
|
||||
# Test orthornomal mode.
|
||||
for i in range(len(X)):
|
||||
x = np.array(X[i], dtype=self.rdt)
|
||||
dt = np.result_type(np.float32, self.rdt)
|
||||
y = dct(x, norm='ortho', type=1)
|
||||
y2 = naive_dct1(x, norm='ortho')
|
||||
assert_equal(y.dtype, dt)
|
||||
assert_array_almost_equal(y / np.max(y), y2 / np.max(y), decimal=self.dec)
|
||||
|
||||
class _TestDCTIIBase(_TestDCTBase):
|
||||
def test_definition_matlab(self):
|
||||
# Test correspondence with matlab (orthornomal mode).
|
||||
for i in range(len(X)):
|
||||
dt = np.result_type(np.float32, self.rdt)
|
||||
x = np.array(X[i], dtype=dt)
|
||||
|
||||
yr = Y[i]
|
||||
y = dct(x, norm="ortho", type=2)
|
||||
assert_equal(y.dtype, dt)
|
||||
assert_array_almost_equal(y, yr, decimal=self.dec)
|
||||
|
||||
|
||||
class _TestDCTIIIBase(_TestDCTBase):
|
||||
def test_definition_ortho(self):
|
||||
# Test orthornomal mode.
|
||||
for i in range(len(X)):
|
||||
x = np.array(X[i], dtype=self.rdt)
|
||||
dt = np.result_type(np.float32, self.rdt)
|
||||
y = dct(x, norm='ortho', type=2)
|
||||
xi = dct(y, norm="ortho", type=3)
|
||||
assert_equal(xi.dtype, dt)
|
||||
assert_array_almost_equal(xi, x, decimal=self.dec)
|
||||
|
||||
class _TestDCTIVBase(_TestDCTBase):
|
||||
def test_definition_ortho(self):
|
||||
# Test orthornomal mode.
|
||||
for i in range(len(X)):
|
||||
x = np.array(X[i], dtype=self.rdt)
|
||||
dt = np.result_type(np.float32, self.rdt)
|
||||
y = dct(x, norm='ortho', type=4)
|
||||
y2 = naive_dct4(x, norm='ortho')
|
||||
assert_equal(y.dtype, dt)
|
||||
assert_array_almost_equal(y / np.max(y), y2 / np.max(y), decimal=self.dec)
|
||||
|
||||
|
||||
class TestDCTIDouble(_TestDCTIBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.double
|
||||
self.dec = 10
|
||||
self.type = 1
|
||||
|
||||
|
||||
class TestDCTIFloat(_TestDCTIBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.float32
|
||||
self.dec = 4
|
||||
self.type = 1
|
||||
|
||||
|
||||
class TestDCTIInt(_TestDCTIBase):
|
||||
def setup_method(self):
|
||||
self.rdt = int
|
||||
self.dec = 5
|
||||
self.type = 1
|
||||
|
||||
|
||||
class TestDCTIIDouble(_TestDCTIIBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.double
|
||||
self.dec = 10
|
||||
self.type = 2
|
||||
|
||||
|
||||
class TestDCTIIFloat(_TestDCTIIBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.float32
|
||||
self.dec = 5
|
||||
self.type = 2
|
||||
|
||||
|
||||
class TestDCTIIInt(_TestDCTIIBase):
|
||||
def setup_method(self):
|
||||
self.rdt = int
|
||||
self.dec = 5
|
||||
self.type = 2
|
||||
|
||||
|
||||
class TestDCTIIIDouble(_TestDCTIIIBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.double
|
||||
self.dec = 14
|
||||
self.type = 3
|
||||
|
||||
|
||||
class TestDCTIIIFloat(_TestDCTIIIBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.float32
|
||||
self.dec = 5
|
||||
self.type = 3
|
||||
|
||||
|
||||
class TestDCTIIIInt(_TestDCTIIIBase):
|
||||
def setup_method(self):
|
||||
self.rdt = int
|
||||
self.dec = 5
|
||||
self.type = 3
|
||||
|
||||
|
||||
class TestDCTIVDouble(_TestDCTIVBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.double
|
||||
self.dec = 12
|
||||
self.type = 3
|
||||
|
||||
|
||||
class TestDCTIVFloat(_TestDCTIVBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.float32
|
||||
self.dec = 5
|
||||
self.type = 3
|
||||
|
||||
|
||||
class TestDCTIVInt(_TestDCTIVBase):
|
||||
def setup_method(self):
|
||||
self.rdt = int
|
||||
self.dec = 5
|
||||
self.type = 3
|
||||
|
||||
|
||||
class _TestIDCTBase(object):
|
||||
def setup_method(self):
|
||||
self.rdt = None
|
||||
self.dec = 14
|
||||
self.type = None
|
||||
|
||||
def test_definition(self):
|
||||
for i in FFTWDATA_SIZES:
|
||||
xr, yr, dt = fftw_dct_ref(self.type, i, self.rdt)
|
||||
x = idct(yr, type=self.type)
|
||||
if self.type == 1:
|
||||
x /= 2 * (i-1)
|
||||
else:
|
||||
x /= 2 * i
|
||||
assert_equal(x.dtype, dt)
|
||||
# XXX: we divide by np.max(y) because the tests fail otherwise. We
|
||||
# should really use something like assert_array_approx_equal. The
|
||||
# difference is due to fftw using a better algorithm w.r.t error
|
||||
# propagation compared to the ones from fftpack.
|
||||
assert_array_almost_equal(x / np.max(x), xr / np.max(x), decimal=self.dec,
|
||||
err_msg="Size %d failed" % i)
|
||||
|
||||
|
||||
class TestIDCTIDouble(_TestIDCTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.double
|
||||
self.dec = 10
|
||||
self.type = 1
|
||||
|
||||
|
||||
class TestIDCTIFloat(_TestIDCTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.float32
|
||||
self.dec = 4
|
||||
self.type = 1
|
||||
|
||||
|
||||
class TestIDCTIInt(_TestIDCTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = int
|
||||
self.dec = 4
|
||||
self.type = 1
|
||||
|
||||
|
||||
class TestIDCTIIDouble(_TestIDCTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.double
|
||||
self.dec = 10
|
||||
self.type = 2
|
||||
|
||||
|
||||
class TestIDCTIIFloat(_TestIDCTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.float32
|
||||
self.dec = 5
|
||||
self.type = 2
|
||||
|
||||
|
||||
class TestIDCTIIInt(_TestIDCTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = int
|
||||
self.dec = 5
|
||||
self.type = 2
|
||||
|
||||
|
||||
class TestIDCTIIIDouble(_TestIDCTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.double
|
||||
self.dec = 14
|
||||
self.type = 3
|
||||
|
||||
|
||||
class TestIDCTIIIFloat(_TestIDCTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.float32
|
||||
self.dec = 5
|
||||
self.type = 3
|
||||
|
||||
|
||||
class TestIDCTIIIInt(_TestIDCTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = int
|
||||
self.dec = 5
|
||||
self.type = 3
|
||||
|
||||
class TestIDCTIVDouble(_TestIDCTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.double
|
||||
self.dec = 12
|
||||
self.type = 4
|
||||
|
||||
|
||||
class TestIDCTIVFloat(_TestIDCTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.float32
|
||||
self.dec = 5
|
||||
self.type = 4
|
||||
|
||||
|
||||
class TestIDCTIVInt(_TestIDCTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = int
|
||||
self.dec = 5
|
||||
self.type = 4
|
||||
|
||||
class _TestDSTBase(object):
|
||||
def setup_method(self):
|
||||
self.rdt = None # dtype
|
||||
self.dec = None # number of decimals to match
|
||||
self.type = None # dst type
|
||||
|
||||
def test_definition(self):
|
||||
for i in FFTWDATA_SIZES:
|
||||
xr, yr, dt = fftw_dst_ref(self.type, i, self.rdt)
|
||||
y = dst(xr, type=self.type)
|
||||
assert_equal(y.dtype, dt)
|
||||
# XXX: we divide by np.max(y) because the tests fail otherwise. We
|
||||
# should really use something like assert_array_approx_equal. The
|
||||
# difference is due to fftw using a better algorithm w.r.t error
|
||||
# propagation compared to the ones from fftpack.
|
||||
assert_array_almost_equal(y / np.max(y), yr / np.max(y), decimal=self.dec,
|
||||
err_msg="Size %d failed" % i)
|
||||
|
||||
|
||||
class _TestDSTIBase(_TestDSTBase):
|
||||
def test_definition_ortho(self):
|
||||
# Test orthornomal mode.
|
||||
for i in range(len(X)):
|
||||
x = np.array(X[i], dtype=self.rdt)
|
||||
dt = np.result_type(np.float32, self.rdt)
|
||||
y = dst(x, norm='ortho', type=1)
|
||||
y2 = naive_dst1(x, norm='ortho')
|
||||
assert_equal(y.dtype, dt)
|
||||
assert_array_almost_equal(y / np.max(y), y2 / np.max(y), decimal=self.dec)
|
||||
|
||||
class _TestDSTIVBase(_TestDSTBase):
|
||||
def test_definition_ortho(self):
|
||||
# Test orthornomal mode.
|
||||
for i in range(len(X)):
|
||||
x = np.array(X[i], dtype=self.rdt)
|
||||
dt = np.result_type(np.float32, self.rdt)
|
||||
y = dst(x, norm='ortho', type=4)
|
||||
y2 = naive_dst4(x, norm='ortho')
|
||||
assert_equal(y.dtype, dt)
|
||||
assert_array_almost_equal(y, y2, decimal=self.dec)
|
||||
|
||||
class TestDSTIDouble(_TestDSTIBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.double
|
||||
self.dec = 12
|
||||
self.type = 1
|
||||
|
||||
|
||||
class TestDSTIFloat(_TestDSTIBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.float32
|
||||
self.dec = 4
|
||||
self.type = 1
|
||||
|
||||
|
||||
class TestDSTIInt(_TestDSTIBase):
|
||||
def setup_method(self):
|
||||
self.rdt = int
|
||||
self.dec = 5
|
||||
self.type = 1
|
||||
|
||||
|
||||
class TestDSTIIDouble(_TestDSTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.double
|
||||
self.dec = 14
|
||||
self.type = 2
|
||||
|
||||
|
||||
class TestDSTIIFloat(_TestDSTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.float32
|
||||
self.dec = 6
|
||||
self.type = 2
|
||||
|
||||
|
||||
class TestDSTIIInt(_TestDSTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = int
|
||||
self.dec = 6
|
||||
self.type = 2
|
||||
|
||||
|
||||
class TestDSTIIIDouble(_TestDSTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.double
|
||||
self.dec = 14
|
||||
self.type = 3
|
||||
|
||||
|
||||
class TestDSTIIIFloat(_TestDSTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.float32
|
||||
self.dec = 7
|
||||
self.type = 3
|
||||
|
||||
|
||||
class TestDSTIIIInt(_TestDSTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = int
|
||||
self.dec = 7
|
||||
self.type = 3
|
||||
|
||||
|
||||
class TestDSTIVDouble(_TestDSTIVBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.double
|
||||
self.dec = 12
|
||||
self.type = 4
|
||||
|
||||
|
||||
class TestDSTIVFloat(_TestDSTIVBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.float32
|
||||
self.dec = 4
|
||||
self.type = 4
|
||||
|
||||
|
||||
class TestDSTIVInt(_TestDSTIVBase):
|
||||
def setup_method(self):
|
||||
self.rdt = int
|
||||
self.dec = 5
|
||||
self.type = 4
|
||||
|
||||
|
||||
class _TestIDSTBase(object):
|
||||
def setup_method(self):
|
||||
self.rdt = None
|
||||
self.dec = None
|
||||
self.type = None
|
||||
|
||||
def test_definition(self):
|
||||
for i in FFTWDATA_SIZES:
|
||||
xr, yr, dt = fftw_dst_ref(self.type, i, self.rdt)
|
||||
x = idst(yr, type=self.type)
|
||||
if self.type == 1:
|
||||
x /= 2 * (i+1)
|
||||
else:
|
||||
x /= 2 * i
|
||||
assert_equal(x.dtype, dt)
|
||||
# XXX: we divide by np.max(x) because the tests fail otherwise. We
|
||||
# should really use something like assert_array_approx_equal. The
|
||||
# difference is due to fftw using a better algorithm w.r.t error
|
||||
# propagation compared to the ones from fftpack.
|
||||
assert_array_almost_equal(x / np.max(x), xr / np.max(x), decimal=self.dec,
|
||||
err_msg="Size %d failed" % i)
|
||||
|
||||
|
||||
class TestIDSTIDouble(_TestIDSTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.double
|
||||
self.dec = 12
|
||||
self.type = 1
|
||||
|
||||
|
||||
class TestIDSTIFloat(_TestIDSTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.float32
|
||||
self.dec = 4
|
||||
self.type = 1
|
||||
|
||||
|
||||
class TestIDSTIInt(_TestIDSTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = int
|
||||
self.dec = 4
|
||||
self.type = 1
|
||||
|
||||
|
||||
class TestIDSTIIDouble(_TestIDSTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.double
|
||||
self.dec = 14
|
||||
self.type = 2
|
||||
|
||||
|
||||
class TestIDSTIIFloat(_TestIDSTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.float32
|
||||
self.dec = 6
|
||||
self.type = 2
|
||||
|
||||
|
||||
class TestIDSTIIInt(_TestIDSTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = int
|
||||
self.dec = 6
|
||||
self.type = 2
|
||||
|
||||
|
||||
class TestIDSTIIIDouble(_TestIDSTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.double
|
||||
self.dec = 14
|
||||
self.type = 3
|
||||
|
||||
|
||||
class TestIDSTIIIFloat(_TestIDSTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.float32
|
||||
self.dec = 6
|
||||
self.type = 3
|
||||
|
||||
|
||||
class TestIDSTIIIInt(_TestIDSTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = int
|
||||
self.dec = 6
|
||||
self.type = 3
|
||||
|
||||
|
||||
class TestIDSTIVDouble(_TestIDSTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.double
|
||||
self.dec = 12
|
||||
self.type = 4
|
||||
|
||||
|
||||
class TestIDSTIVFloat(_TestIDSTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = np.float32
|
||||
self.dec = 6
|
||||
self.type = 4
|
||||
|
||||
|
||||
class TestIDSTIVnt(_TestIDSTBase):
|
||||
def setup_method(self):
|
||||
self.rdt = int
|
||||
self.dec = 6
|
||||
self.type = 4
|
||||
|
||||
|
||||
class TestOverwrite(object):
|
||||
"""Check input overwrite behavior."""
|
||||
|
||||
real_dtypes = [np.float32, np.float64]
|
||||
|
||||
def _check(self, x, routine, type, fftsize, axis, norm, overwrite_x,
|
||||
should_overwrite, **kw):
|
||||
x2 = x.copy()
|
||||
routine(x2, type, fftsize, axis, norm, overwrite_x=overwrite_x)
|
||||
|
||||
sig = "%s(%s%r, %r, axis=%r, overwrite_x=%r)" % (
|
||||
routine.__name__, x.dtype, x.shape, fftsize, axis, overwrite_x)
|
||||
if not should_overwrite:
|
||||
assert_equal(x2, x, err_msg="spurious overwrite in %s" % sig)
|
||||
|
||||
def _check_1d(self, routine, dtype, shape, axis, overwritable_dtypes):
|
||||
np.random.seed(1234)
|
||||
if np.issubdtype(dtype, np.complexfloating):
|
||||
data = np.random.randn(*shape) + 1j*np.random.randn(*shape)
|
||||
else:
|
||||
data = np.random.randn(*shape)
|
||||
data = data.astype(dtype)
|
||||
|
||||
for type in [1, 2, 3, 4]:
|
||||
for overwrite_x in [True, False]:
|
||||
for norm in [None, 'ortho']:
|
||||
should_overwrite = (overwrite_x
|
||||
and dtype in overwritable_dtypes
|
||||
and (len(shape) == 1 or
|
||||
(axis % len(shape) == len(shape)-1
|
||||
)))
|
||||
self._check(data, routine, type, None, axis, norm,
|
||||
overwrite_x, should_overwrite)
|
||||
|
||||
def test_dct(self):
|
||||
overwritable = self.real_dtypes
|
||||
for dtype in self.real_dtypes:
|
||||
self._check_1d(dct, dtype, (16,), -1, overwritable)
|
||||
self._check_1d(dct, dtype, (16, 2), 0, overwritable)
|
||||
self._check_1d(dct, dtype, (2, 16), 1, overwritable)
|
||||
|
||||
def test_idct(self):
|
||||
overwritable = self.real_dtypes
|
||||
for dtype in self.real_dtypes:
|
||||
self._check_1d(idct, dtype, (16,), -1, overwritable)
|
||||
self._check_1d(idct, dtype, (16, 2), 0, overwritable)
|
||||
self._check_1d(idct, dtype, (2, 16), 1, overwritable)
|
||||
|
||||
def test_dst(self):
|
||||
overwritable = self.real_dtypes
|
||||
for dtype in self.real_dtypes:
|
||||
self._check_1d(dst, dtype, (16,), -1, overwritable)
|
||||
self._check_1d(dst, dtype, (16, 2), 0, overwritable)
|
||||
self._check_1d(dst, dtype, (2, 16), 1, overwritable)
|
||||
|
||||
def test_idst(self):
|
||||
overwritable = self.real_dtypes
|
||||
for dtype in self.real_dtypes:
|
||||
self._check_1d(idst, dtype, (16,), -1, overwritable)
|
||||
self._check_1d(idst, dtype, (16, 2), 0, overwritable)
|
||||
self._check_1d(idst, dtype, (2, 16), 1, overwritable)
|
||||
|
||||
|
||||
class Test_DCTN_IDCTN(object):
|
||||
dec = 14
|
||||
dct_type = [1, 2, 3, 4]
|
||||
norms = [None, 'ortho']
|
||||
rstate = np.random.RandomState(1234)
|
||||
shape = (32, 16)
|
||||
data = rstate.randn(*shape)
|
||||
|
||||
@pytest.mark.parametrize('fforward,finverse', [(dctn, idctn),
|
||||
(dstn, idstn)])
|
||||
@pytest.mark.parametrize('axes', [None,
|
||||
1, (1,), [1],
|
||||
0, (0,), [0],
|
||||
(0, 1), [0, 1],
|
||||
(-2, -1), [-2, -1]])
|
||||
@pytest.mark.parametrize('dct_type', dct_type)
|
||||
@pytest.mark.parametrize('norm', ['ortho'])
|
||||
def test_axes_round_trip(self, fforward, finverse, axes, dct_type, norm):
|
||||
tmp = fforward(self.data, type=dct_type, axes=axes, norm=norm)
|
||||
tmp = finverse(tmp, type=dct_type, axes=axes, norm=norm)
|
||||
assert_array_almost_equal(self.data, tmp, decimal=12)
|
||||
|
||||
@pytest.mark.parametrize('fforward,fforward_ref', [(dctn, dct_2d_ref),
|
||||
(dstn, dst_2d_ref)])
|
||||
@pytest.mark.parametrize('dct_type', dct_type)
|
||||
@pytest.mark.parametrize('norm', norms)
|
||||
def test_dctn_vs_2d_reference(self, fforward, fforward_ref,
|
||||
dct_type, norm):
|
||||
y1 = fforward(self.data, type=dct_type, axes=None, norm=norm)
|
||||
y2 = fforward_ref(self.data, type=dct_type, norm=norm)
|
||||
assert_array_almost_equal(y1, y2, decimal=11)
|
||||
|
||||
@pytest.mark.parametrize('finverse,finverse_ref', [(idctn, idct_2d_ref),
|
||||
(idstn, idst_2d_ref)])
|
||||
@pytest.mark.parametrize('dct_type', dct_type)
|
||||
@pytest.mark.parametrize('norm', [None, 'ortho'])
|
||||
def test_idctn_vs_2d_reference(self, finverse, finverse_ref,
|
||||
dct_type, norm):
|
||||
fdata = dctn(self.data, type=dct_type, norm=norm)
|
||||
y1 = finverse(fdata, type=dct_type, norm=norm)
|
||||
y2 = finverse_ref(fdata, type=dct_type, norm=norm)
|
||||
assert_array_almost_equal(y1, y2, decimal=11)
|
||||
|
||||
@pytest.mark.parametrize('fforward,finverse', [(dctn, idctn),
|
||||
(dstn, idstn)])
|
||||
def test_axes_and_shape(self, fforward, finverse):
|
||||
with assert_raises(ValueError,
|
||||
match="when given, axes and shape arguments"
|
||||
" have to be of the same length"):
|
||||
fforward(self.data, shape=self.data.shape[0], axes=(0, 1))
|
||||
|
||||
with assert_raises(ValueError,
|
||||
match="when given, axes and shape arguments"
|
||||
" have to be of the same length"):
|
||||
fforward(self.data, shape=self.data.shape[0], axes=None)
|
||||
|
||||
with assert_raises(ValueError,
|
||||
match="when given, axes and shape arguments"
|
||||
" have to be of the same length"):
|
||||
fforward(self.data, shape=self.data.shape, axes=0)
|
||||
|
||||
@pytest.mark.parametrize('fforward', [dctn, dstn])
|
||||
def test_shape(self, fforward):
|
||||
tmp = fforward(self.data, shape=(128, 128), axes=None)
|
||||
assert_equal(tmp.shape, (128, 128))
|
||||
|
||||
@pytest.mark.parametrize('fforward,finverse', [(dctn, idctn),
|
||||
(dstn, idstn)])
|
||||
@pytest.mark.parametrize('axes', [1, (1,), [1],
|
||||
0, (0,), [0]])
|
||||
def test_shape_is_none_with_axes(self, fforward, finverse, axes):
|
||||
tmp = fforward(self.data, shape=None, axes=axes, norm='ortho')
|
||||
tmp = finverse(tmp, shape=None, axes=axes, norm='ortho')
|
||||
assert_array_almost_equal(self.data, tmp, decimal=self.dec)
|
||||
Reference in New Issue
Block a user