diff --git a/pywafo/src/wafo/dctpack.py b/pywafo/src/wafo/dctpack.py
index de9aa72..79dbd70 100644
--- a/pywafo/src/wafo/dctpack.py
+++ b/pywafo/src/wafo/dctpack.py
@@ -1,112 +1,150 @@
import numpy as np
+from scipy.fftpack import dct as _dct
+from scipy.fftpack import idct as _idct
__all__ = ['dct', 'idct', 'dctn', 'idctn']
-def dct(x, n=None):
- """
- Discrete Cosine Transform
- N-1
- y[k] = 2* sum x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N.
- n=0
-
- Examples
+def dct(x, type=2, n=None, axis=-1, norm='ortho'):
+ '''
+ Return the Discrete Cosine Transform of arbitrary type sequence x.
+
+ Parameters
+ ----------
+ x : array_like
+ The input array.
+ type : {1, 2, 3}, optional
+ Type of the DCT (see Notes). Default type is 2.
+ n : int, optional
+ Length of the transform.
+ axis : int, optional
+ Axis over which to compute the transform.
+ norm : {None, 'ortho'}, optional
+ Normalization mode (see Notes). Default is 'ortho'.
+
+ Returns
+ -------
+ y : ndarray of real
+ The transformed input array.
+
+ See Also
--------
- >>> import numpy as np
- >>> x = np.arange(5)
- >>> np.abs(x-idct(dct(x)))<1e-14
- array([ True, True, True, True, True], dtype=bool)
- >>> np.abs(x-dct(idct(x)))<1e-14
- array([ True, True, True, True, True], dtype=bool)
-
- Reference
- ---------
- http://en.wikipedia.org/wiki/Discrete_cosine_transform
- http://users.ece.utexas.edu/~bevans/courses/ee381k/lectures/
- """
- fft = np.fft.fft
- x = np.atleast_1d(x)
-
- if n is None:
- n = x.shape[-1]
-
- if x.shape[-1] < n:
- n_shape = x.shape[:-1] + (n - x.shape[-1],)
- xx = np.hstack((x, np.zeros(n_shape)))
- else:
- xx = x[..., :n]
-
- real_x = np.all(np.isreal(xx))
- if (real_x and (np.remainder(n, 2) == 0)):
- xp = 2 * fft(np.hstack((xx[..., ::2], xx[..., ::-2])))
- else:
- xp = fft(np.hstack((xx, xx[..., ::-1])))
- xp = xp[..., :n]
-
- w = np.exp(-1j * np.arange(n) * np.pi / (2 * n))
-
- y = xp * w
-
- if real_x:
- return y.real
- else:
- return y
-
-def idct(x, n=None):
- """
- Inverse Discrete Cosine Transform
-
+ idct
+
+ 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 3 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
+
+ Only None is supported as normalization mode for DCT-I. Note also that 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
- x[k] = 1/N sum w[n]*y[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N.
- n=0
-
- w(0) = 1/2
- w(n) = 1 for n>0
-
- Examples
- --------
- >>> import numpy as np
- >>> x = np.arange(5)
- >>> np.abs(x-idct(dct(x)))<1e-14
- array([ True, True, True, True, True], dtype=bool)
- >>> np.abs(x-dct(idct(x)))<1e-14
- array([ True, True, True, True, True], dtype=bool)
-
- Reference
- ---------
+ 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(1/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.
+
+ References
+ ----------
+
http://en.wikipedia.org/wiki/Discrete_cosine_transform
- http://users.ece.utexas.edu/~bevans/courses/ee381k/lectures/
- """
-
- ifft = np.fft.ifft
- x = np.atleast_1d(x)
-
- if n is None:
- n = x.shape[-1]
-
- w = np.exp(1j * np.arange(n) * np.pi / (2 * n))
-
- if x.shape[-1] < n:
- n_shape = x.shape[:-1] + (n - x.shape[-1],)
- xx = np.hstack((x, np.zeros(n_shape))) * w
- else:
- xx = x[..., :n] * w
-
- real_x = np.all(np.isreal(x))
- if (real_x and (np.remainder(n, 2) == 0)):
- xx[..., 0] = xx[..., 0] * 0.5
- yp = ifft(xx)
- y = np.zeros(xx.shape, dtype=complex)
- y[..., ::2] = yp[..., :n / 2]
- y[..., ::-2] = yp[..., n / 2::]
+
+ '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, http://dx.doi.org/10.1109/TASSP.1980.1163351 (1980).
+ '''
+ farr = np.asfarray
+ if np.iscomplex(x).any():
+ return _dct(farr(x.real), type, n, axis, norm) + 1j*_dct(farr(x.imag), type, n, axis, norm)
else:
- yp = ifft(np.hstack((xx, np.zeros_like(xx[..., 0]), np.conj(xx[..., :0:-1]))))
- y = yp[..., :n]
+ return _dct(farr(x), type, n, axis, norm)
+def idct(x, type=2, n=None, axis=-1, norm='ortho'):
+ '''
+ Return the Inverse Discrete Cosine Transform of an arbitrary type sequence.
- if real_x:
- return y.real
+ Parameters
+ ----------
+ x : array_like
+ The input array.
+ type : {1, 2, 3}, optional
+ Type of the DCT (see Notes). Default type is 2.
+ n : int, optional
+ Length of the transform.
+ axis : int, optional
+ Axis over which to compute the transform.
+ norm : {None, 'ortho'}, optional
+ Normalization mode (see Notes). Default is 'ortho'.
+
+ Returns
+ -------
+ y : ndarray of real
+ The transformed input array.
+
+ See Also
+ --------
+ 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. For the definition of these types,
+ see `dct`.
+ '''
+ farr = np.asarray
+ if np.iscomplex(x).any():
+ return _idct(farr(x.real), type, n, axis, norm) + 1j*_idct(farr(x.imag), type, n, axis, norm)
else:
- return y
-
-def dctn(y, axis=None, w=None):
+ return _idct(farr(x), type, n, axis, norm)
+def dctn(x, type=2, axis=None, norm='ortho'):
'''
DCTN N-D discrete cosine transform.
@@ -145,172 +183,381 @@ def dctn(y, axis=None, w=None):
--------
idctn, dct, idct
'''
-
- y = np.atleast_1d(y)
- shape0 = y.shape
-
+ y = np.atleast_1d(x)
+ shape0 = y.shape
if axis is None:
y = y.squeeze() # Working across singleton dimensions is useless
- dimy = y.ndim
- if dimy==1:
- y = np.atleast_2d(y)
- y = y.T
- # Some modifications are required if Y is a vector
-# if isvector(y):
-# if y.shape[0]==1:
-# if axis==0:
-# return y, None
-# elif axis==1:
-# axis=0
-# y = y.T
-# elif axis==1:
-# return y, None
-
- if w is None:
- w = [0,] * dimy
- for dim in range(dimy):
- if axis is not None and dim!=axis:
- continue
- n = (dimy==1)*y.size + (dimy>1)*shape0[dim]
- #w{dim} = exp(1i*(0:n-1)'*pi/2/n);
- w[dim] = np.exp(1j * np.arange(n) * np.pi / (2 * n))
-
- # --- DCT algorithm ---
+ ndim = y.ndim
+ isvector = max(shape0)==y.size
+ if isvector:
+ if ndim==1:
+ y = np.atleast_2d(y)
+ y = y.T
+ elif y.shape[0]==1:
+ if axis==0:
+ return x
+ elif axis==1:
+ axis=0
+ y = y.T
+ elif axis==1:
+ return y
+
if np.iscomplex(y).any():
- y = np.complex(dctn(np.real(y),axis,w),dctn(np.imag(y),axis,w))
+ y = dctn(y.real, type, axis, norm) + 1j*dctn(y.imag, type, axis, norm)
else:
- for dim in range(dimy):
- y = shiftdim(y,1)
+ y = np.asfarray(y)
+ for dim in range(ndim):
+ y = y.transpose(np.roll(range(y.ndim), -1))
+ #y = shiftdim(y,1)
if axis is not None and dim!=axis:
- y = shiftdim(y, 1)
continue
- siz = y.shape
- n = siz[-1]
- y = y[...,np.r_[0:n:2, 2*int(n//2)-1:0:-2]]
- y = y.reshape((-1,n))
- y = y*np.sqrt(2*n);
- y = (np.fft.ifft(y, n=n, axis=1) * w[dim]).real
- y[:,0] = y[:,0]/np.sqrt(2)
- y = y.reshape(siz)
-
- #end
- #end
-
- return y.reshape(shape0), w
-
-def idctn(y, axis=None, w=None):
- '''
- IDCTN N-D inverse discrete cosine transform.
- X = IDCTN(Y) inverts the N-D DCT transform, returning the original
- array if Y was obtained using Y = DCTN(X).
-
- IDCTN(X,DIM) applies the IDCTN operation across the dimension DIM.
-
- Class Support
- -------------
- Input array can be numeric or logical. The returned array is of class
- double.
-
- Reference
- ---------
- Narasimha M. et al, On the computation of the discrete cosine
- transform, IEEE Trans Comm, 26, 6, 1978, pp 934-936.
-
- Example
- -------
- RGB = imread('autumn.tif');
- I = rgb2gray(RGB);
- J = dctn(I);
- imshow(log(abs(J)),[]), colormap(jet), colorbar
-
- The commands below set values less than magnitude 10 in the DCT matrix
- to zero, then reconstruct the image using the inverse DCT.
-
- J(abs(J)<10) = 0;
- K = idctn(J);
- figure, imshow(I)
- figure, imshow(K,[0 255])
-
- See also
- --------
- dctn, idct, dct
-
- -- Damien Garcia -- 2009/04, revised 2009/11
- website: www.BiomeCardio.com
-
- ----------
- [Y,W] = IDCTN(X,DIM,W) uses and returns the weights which are used by
- the program. If IDCTN is required for several large arrays of same
- size, the weights can be reused to make the algorithm faster. A typical
- syntax is the following:
- w = [];
- for k = 1:10
- [y{k},w] = idctn(x{k},[],w);
- end
- The weights (w) are calculated during the first call of IDCTN then
- reused in the next calls.
- '''
-
- y = np.atleast_1d(y)
- shape0 = y.shape
+ y = _dct(y, type, norm=norm)
+ return y.reshape(shape0)
+def idctn(x, type=2, axis=None, norm='ortho'):
+ y = np.atleast_1d(x)
+ shape0 = y.shape
if axis is None:
y = y.squeeze() # Working across singleton dimensions is useless
-
- dimy = y.ndim
- if dimy==1:
- y = np.atleast_2d(y)
- y = y.T
- # Some modifications are required if Y is a vector
-# if isvector(y):
-# if y.shape[0]==1:
-# if axis==0:
-# return y, None
-# elif axis==1:
-# axis=0
-# y = y.T
-# elif axis==1:
-# return y, None
-#
-
-
- if w is None:
- w = [0,] * dimy
- for dim in range(dimy):
- if axis is not None and dim!=axis:
- continue
- n = (dimy==1)*y.size + (dimy>1)*shape0[dim]
- #w{dim} = exp(1i*(0:n-1)'*pi/2/n);
- w[dim] = np.exp(1j * np.arange(n) * np.pi / (2 * n))
- # --- IDCT algorithm ---
+ ndim = y.ndim
+ isvector = max(shape0)==y.size
+ if isvector:
+ if ndim==1:
+ y = np.atleast_2d(y)
+ y = y.T
+ elif y.shape[0]==1:
+ if axis==0:
+ return x
+ elif axis==1:
+ axis=0
+ y = y.T
+ elif axis==1:
+ return y
+
if np.iscomplex(y).any():
- y = np.complex(idctn(np.real(y),axis,w),idctn(np.imag(y),axis,w))
+ y = idctn(y.real, type, axis, norm) + 1j*idctn(y.imag, type, axis, norm)
else:
- for dim in range(dimy):
- y = shiftdim(y,1)
+ y = np.asfarray(y)
+ for dim in range(ndim):
+ y = y.transpose(np.roll(range(y.ndim), -1))
+ #y = shiftdim(y,1)
if axis is not None and dim!=axis:
- #y = shiftdim(y, 1)
continue
- siz = y.shape
- n = siz[-1]
-
- y = y.reshape((-1,n)) * w[dim]
- y[:,0] = y[:,0]/np.sqrt(2)
- y = (np.fft.ifft(y, n=n, axis=1)).real
- y = y * np.sqrt(2*n)
-
- I = np.empty(n,dtype=int)
- I.put(np.r_[0:n:2],np.r_[0:int(n//2)+np.remainder(n,2)])
- I.put(np.r_[1:n:2],np.r_[n-1:int(n//2)-1:-1])
- y = y[:,I]
-
- y = y.reshape(siz)
-
-
- y = y.reshape(shape0);
- return y, w
+ y = _idct(y, type, norm=norm)
+ return y.reshape(shape0)
+
+#def dct(x, n=None):
+# """
+# Discrete Cosine Transform
+#
+# N-1
+# y[k] = 2* sum x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N.
+# n=0
+#
+# Examples
+# --------
+# >>> import numpy as np
+# >>> x = np.arange(5)
+# >>> np.abs(x-idct(dct(x)))<1e-14
+# array([ True, True, True, True, True], dtype=bool)
+# >>> np.abs(x-dct(idct(x)))<1e-14
+# array([ True, True, True, True, True], dtype=bool)
+#
+# Reference
+# ---------
+# http://en.wikipedia.org/wiki/Discrete_cosine_transform
+# http://users.ece.utexas.edu/~bevans/courses/ee381k/lectures/
+# """
+# fft = np.fft.fft
+# x = np.atleast_1d(x)
+#
+# if n is None:
+# n = x.shape[-1]
+#
+# if x.shape[-1] < n:
+# n_shape = x.shape[:-1] + (n - x.shape[-1],)
+# xx = np.hstack((x, np.zeros(n_shape)))
+# else:
+# xx = x[..., :n]
+#
+# real_x = np.all(np.isreal(xx))
+# if (real_x and (np.remainder(n, 2) == 0)):
+# xp = 2 * fft(np.hstack((xx[..., ::2], xx[..., ::-2])))
+# else:
+# xp = fft(np.hstack((xx, xx[..., ::-1])))
+# xp = xp[..., :n]
+#
+# w = np.exp(-1j * np.arange(n) * np.pi / (2 * n))
+#
+# y = xp * w
+#
+# if real_x:
+# return y.real
+# else:
+# return y
+#
+#def idct(x, n=None):
+# """
+# Inverse Discrete Cosine Transform
+#
+# N-1
+# x[k] = 1/N sum w[n]*y[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N.
+# n=0
+#
+# w(0) = 1/2
+# w(n) = 1 for n>0
+#
+# Examples
+# --------
+# >>> import numpy as np
+# >>> x = np.arange(5)
+# >>> np.abs(x-idct(dct(x)))<1e-14
+# array([ True, True, True, True, True], dtype=bool)
+# >>> np.abs(x-dct(idct(x)))<1e-14
+# array([ True, True, True, True, True], dtype=bool)
+#
+# Reference
+# ---------
+# http://en.wikipedia.org/wiki/Discrete_cosine_transform
+# http://users.ece.utexas.edu/~bevans/courses/ee381k/lectures/
+# """
+#
+# ifft = np.fft.ifft
+# x = np.atleast_1d(x)
+#
+# if n is None:
+# n = x.shape[-1]
+#
+# w = np.exp(1j * np.arange(n) * np.pi / (2 * n))
+#
+# if x.shape[-1] < n:
+# n_shape = x.shape[:-1] + (n - x.shape[-1],)
+# xx = np.hstack((x, np.zeros(n_shape))) * w
+# else:
+# xx = x[..., :n] * w
+#
+# real_x = np.all(np.isreal(x))
+# if (real_x and (np.remainder(n, 2) == 0)):
+# xx[..., 0] = xx[..., 0] * 0.5
+# yp = ifft(xx)
+# y = np.zeros(xx.shape, dtype=complex)
+# y[..., ::2] = yp[..., :n / 2]
+# y[..., ::-2] = yp[..., n / 2::]
+# else:
+# yp = ifft(np.hstack((xx, np.zeros_like(xx[..., 0]), np.conj(xx[..., :0:-1]))))
+# y = yp[..., :n]
+#
+# if real_x:
+# return y.real
+# else:
+# return y
+#
+#def dctn(y, axis=None, w=None):
+# '''
+# DCTN N-D discrete cosine transform.
+#
+# Y = DCTN(X) returns the discrete cosine transform of X. The array Y is
+# the same size as X and contains the discrete cosine transform
+# coefficients. This transform can be inverted using IDCTN.
+#
+# DCTN(X,axis) applies the DCTN operation across the dimension axis.
+#
+# Class Support
+# -------------
+# Input array can be numeric or logical. The returned array is of class
+# double.
+#
+# Reference
+# ---------
+# Narasimha M. et al, On the computation of the discrete cosine
+# transform, IEEE Trans Comm, 26, 6, 1978, pp 934-936.
+#
+# Example
+# -------
+# RGB = imread('autumn.tif');
+# I = rgb2gray(RGB);
+# J = dctn(I);
+# imshow(log(abs(J)),[]), colormap(jet), colorbar
+#
+# The commands below set values less than magnitude 10 in the DCT matrix
+# to zero, then reconstruct the image using the inverse DCT.
+#
+# J(abs(J)<10) = 0;
+# K = idctn(J);
+# figure, imshow(I)
+# figure, imshow(K,[0 255])
+#
+# See also
+# --------
+# idctn, dct, idct
+# '''
+#
+# y = np.atleast_1d(y)
+# shape0 = y.shape
+#
+#
+# if axis is None:
+# y = y.squeeze() # Working across singleton dimensions is useless
+# dimy = y.ndim
+# if dimy==1:
+# y = np.atleast_2d(y)
+# y = y.T
+# # Some modifications are required if Y is a vector
+## if isvector(y):
+## if y.shape[0]==1:
+## if axis==0:
+## return y, None
+## elif axis==1:
+## axis=0
+## y = y.T
+## elif axis==1:
+## return y, None
+#
+# if w is None:
+# w = [0,] * dimy
+# for dim in range(dimy):
+# if axis is not None and dim!=axis:
+# continue
+# n = (dimy==1)*y.size + (dimy>1)*shape0[dim]
+# #w{dim} = exp(1i*(0:n-1)'*pi/2/n);
+# w[dim] = np.exp(1j * np.arange(n) * np.pi / (2 * n))
+#
+# # --- DCT algorithm ---
+# if np.iscomplex(y).any():
+# y = dctn(np.real(y),axis,w) + 1j*dctn(np.imag(y),axis,w)
+# else:
+# for dim in range(dimy):
+# y = shiftdim(y,1)
+# if axis is not None and dim!=axis:
+# #y = shiftdim(y, 1)
+# continue
+# siz = y.shape
+# n = siz[-1]
+# y = y[...,np.r_[0:n:2, 2*int(n//2)-1:0:-2]]
+# y = y.reshape((-1,n))
+# y = y*np.sqrt(2*n);
+# y = (np.fft.ifft(y, n=n, axis=1) * w[dim]).real
+# y[:,0] = y[:,0]/np.sqrt(2)
+# y = y.reshape(siz)
+#
+# #end
+# #end
+#
+# return y.reshape(shape0), w
+#
+#def idctn(y, axis=None, w=None):
+# '''
+# IDCTN N-D inverse discrete cosine transform.
+# X = IDCTN(Y) inverts the N-D DCT transform, returning the original
+# array if Y was obtained using Y = DCTN(X).
+#
+# IDCTN(X,DIM) applies the IDCTN operation across the dimension DIM.
+#
+# Class Support
+# -------------
+# Input array can be numeric or logical. The returned array is of class
+# double.
+#
+# Reference
+# ---------
+# Narasimha M. et al, On the computation of the discrete cosine
+# transform, IEEE Trans Comm, 26, 6, 1978, pp 934-936.
+#
+# Example
+# -------
+# RGB = imread('autumn.tif');
+# I = rgb2gray(RGB);
+# J = dctn(I);
+# imshow(log(abs(J)),[]), colormap(jet), colorbar
+#
+# The commands below set values less than magnitude 10 in the DCT matrix
+# to zero, then reconstruct the image using the inverse DCT.
+#
+# J(abs(J)<10) = 0;
+# K = idctn(J);
+# figure, imshow(I)
+# figure, imshow(K,[0 255])
+#
+# See also
+# --------
+# dctn, idct, dct
+#
+# -- Damien Garcia -- 2009/04, revised 2009/11
+# website: www.BiomeCardio.com
+#
+# ----------
+# [Y,W] = IDCTN(X,DIM,W) uses and returns the weights which are used by
+# the program. If IDCTN is required for several large arrays of same
+# size, the weights can be reused to make the algorithm faster. A typical
+# syntax is the following:
+# w = [];
+# for k = 1:10
+# [y{k},w] = idctn(x{k},[],w);
+# end
+# The weights (w) are calculated during the first call of IDCTN then
+# reused in the next calls.
+# '''
+#
+# y = np.atleast_1d(y)
+# shape0 = y.shape
+#
+# if axis is None:
+# y = y.squeeze() # Working across singleton dimensions is useless
+#
+# dimy = y.ndim
+# if dimy==1:
+# y = np.atleast_2d(y)
+# y = y.T
+# # Some modifications are required if Y is a vector
+## if isvector(y):
+## if y.shape[0]==1:
+## if axis==0:
+## return y, None
+## elif axis==1:
+## axis=0
+## y = y.T
+## elif axis==1:
+## return y, None
+##
+#
+#
+# if w is None:
+# w = [0,] * dimy
+# for dim in range(dimy):
+# if axis is not None and dim!=axis:
+# continue
+# n = (dimy==1)*y.size + (dimy>1)*shape0[dim]
+# #w{dim} = exp(1i*(0:n-1)'*pi/2/n);
+# w[dim] = np.exp(1j * np.arange(n) * np.pi / (2 * n))
+# # --- IDCT algorithm ---
+# if np.iscomplex(y).any():
+# y = np.complex(idctn(np.real(y),axis,w),idctn(np.imag(y),axis,w))
+# else:
+# for dim in range(dimy):
+# y = shiftdim(y,1)
+# if axis is not None and dim!=axis:
+# #y = shiftdim(y, 1)
+# continue
+# siz = y.shape
+# n = siz[-1]
+#
+# y = y.reshape((-1,n)) * w[dim]
+# y[:,0] = y[:,0]/np.sqrt(2)
+# y = (np.fft.ifft(y, n=n, axis=1)).real
+# y = y * np.sqrt(2*n)
+#
+# I = np.empty(n,dtype=int)
+# I.put(np.r_[0:n:2],np.r_[0:int(n//2)+np.remainder(n,2)])
+# I.put(np.r_[1:n:2],np.r_[n-1:int(n//2)-1:-1])
+# y = y[:,I]
+#
+# y = y.reshape(siz)
+#
+#
+# y = y.reshape(shape0);
+# return y, w
@@ -323,33 +570,73 @@ def no_leading_ones(x):
return x[first:]
def shiftdim(x, n=None):
+ '''
+ Shift dimensions
+
+ Parameters
+ ----------
+ x : array
+ n : int
+
+ Notes
+ -----
+ Shiftdim is handy for functions that intend to work along the first
+ non-singleton dimension of the array.
+
+ If n is None returns the array with the same number of elements as X but
+ with any leading singleton dimensions removed.
+
+ When n is positive, shiftdim shifts the dimensions to the left and wraps
+ the n leading dimensions to the end.
+
+ When n is negative, shiftdim shifts the dimensions to the right and pads
+ with singletons.
+
+ See also
+ --------
+ reshape, squeeze
+ '''
if n is None:
- # returns the array B with the same number of
- # elements as X but with any leading singleton
- # dimensions removed.
return x.reshape(no_leading_ones(x.shape))
elif n>=0:
- # When n is positive, shiftdim shifts the dimensions
- # to the left and wraps the n leading dimensions to the end.
return x.transpose(np.roll(range(x.ndim), -n))
else:
- # When n is negative, shiftdim shifts the dimensions
- # to the right and pads with singletons.
return x.reshape((1,)*-n+x.shape)
def test_dctn():
- a = np.arange(12).reshape((3,-1))
- #y = dct(a)
- #x = idct(y)
- #print(y)
- #print(x)
+ a = np.arange(12) #.reshape((3,-1))
+ print('a = ', a)
+ print(' ')
+ y = dct(a)
+ x = idct(y)
+ print('y = dct(a)')
+ print(y)
+ print('x = idct(y)')
+ print(x)
+ print(' ')
- print(a)
- yn = dctn(a)[0]
- xn = idctn(yn)[0]
+# y1 = dct1(a)
+# x1 = idct1(y)
+# print('y1 = dct1(a)')
+# print(y1)
+# print('x1 = idct1(y)')
+# print(x1)
+# print(' ')
+ yn = dctn(a)
+ xn = idctn(yn)
+ print('yn = dctn(a)')
print(yn)
+ print('xn = idctn(yn)')
print(xn)
+ print(' ')
+
+# yn1 = dctn1(a)
+# xn1 = idctn1(yn1)
+# print('yn1 = dctn1(a)')
+# print(yn1)
+# print('xn1 = idctn1(yn)')
+# print(xn1)
diff --git a/pywafo/src/wafo/kdetools.py b/pywafo/src/wafo/kdetools.py
index 3e2d11e..64a3f80 100644
--- a/pywafo/src/wafo/kdetools.py
+++ b/pywafo/src/wafo/kdetools.py
@@ -15,9 +15,12 @@ from misc import tranproc #, trangood
from numpy import pi, sqrt, atleast_2d, exp, newaxis #@UnresolvedImport
from scipy import interpolate, linalg, sparse
from scipy.special import gamma
+from scipy.optimize import brentq
from wafo.misc import meshgrid
from wafo.wafodata import WafoData
+from dctpack import dct
+
import copy
import numpy as np
import scipy
@@ -1034,7 +1037,7 @@ class Kernel(object):
'Density estimation for statistics and data analysis'
Chapman and Hall, pp 31, 103, 175
'''
- def __init__(self, name, fun='hns'):
+ def __init__(self, name, fun='hisj'): #'hns'):
self.kernel = _MKERNEL_DICT[name[:4]]
#self.name = self.kernel.__name__.replace('mkernel_', '').title()
try:
@@ -1363,6 +1366,99 @@ class Kernel(object):
h[dim] = h1
#end % for dim loop
return h
+ def hisj(self, data, inc=128):
+ '''
+ HISJ Improved Sheather-Jones estimate of smoothing parameter.
+
+ Unlike many other implementations, this one is immune to problems
+ caused by multimodal densities with widely separated modes. The
+ estimation does not deteriorate for multimodal densities, because
+ it do not assume a parametric model for the data.
+
+ Parameters
+ ----------
+ data - a vector of data from which the density estimate is constructed;
+ n - the number of mesh points used in the uniform discretization
+
+ Returns
+ -------
+ bandwidth - the optimal bandwidth
+
+ Reference
+ ---------
+ Kernel density estimation via diffusion
+ Z. I. Botev, J. F. Grotowski, and D. P. Kroese (2010)
+ Annals of Statistics, Volume 38, Number 5, pages 2916-2957.
+ '''
+ A = np.atleast_2d(data)
+ d, n = A.shape
+
+ # R= int(mkernel(x)^2), mu2= int(x^2*mkernel(x))
+ mu2, R, unusedRdd = self.stats()
+ STEconstant = R / (mu2 ** (2) * n)
+
+ amin = A.min(axis=1) # Find the minimum value of A.
+ amax = A.max(axis=1) # Find the maximum value of A.
+ arange = amax - amin # Find the range of A.
+
+ #% xa holds the x 'axis' vector, defining a grid of x values where
+ #% the k.d. function will be evaluated.
+
+ ax1 = amin - arange / 8.0
+ bx1 = amax + arange / 8.0
+
+ kernel2 = Kernel('gauss')
+ mu2, R, unusedRdd = kernel2.stats()
+ STEconstant2 = R / (mu2 ** (2) * n)
+
+ def fixed_point(t, N, I, a2):
+ ''' this implements the function t-zeta*gamma^[L](t)'''
+
+ prod = np.prod
+ L = 7
+ logI = np.log(I)
+ f = 2 * pi**(2 * L) * (a2 * exp(L * logI -I * pi ** 2 * t)).sum()
+ for s in range(L - 1, 1, -1):
+ K0 = prod(np.r_[1:2*s:2]) / sqrt(2 * pi);
+ const = (1 + (1. / 2) ** (s + 1. / 2)) / 3
+ time = (2 * const * K0 / N / f) ** (2. / (3 + 2 * s))
+ f = 2 * pi ** (2 * s) * (a2 * exp(s * logI-I * pi ** 2 * time)).sum()
+ return t - (2 * N * sqrt(pi) * f) ** (-2. / 5)
+
+ h = np.empty(d)
+ for dim in range(d):
+ ax = ax1[dim]
+ bx = bx1[dim]
+ xa = np.linspace(ax, bx, inc)
+ R = bx-ax
+
+ c = gridcount(A[dim], xa)
+ N = len(set(A[dim]))
+ a = dct(c/c.sum(), norm=None)
+
+ #% now compute the optimal bandwidth^2 using the referenced method
+ I = np.asfarray(np.arange(1, inc))**2
+ a2 = (a[1:]/2)**2
+ fun = lambda t: fixed_point(t, N, I, a2)
+ # use fzero to solve the equation t=zeta*gamma^[5](t)
+ try:
+ t_star = brentq(fun, a=0, b=.1)
+ except:
+ t_star = 0.28*N**(-2./5)
+
+ # smooth the discrete cosine transform of initial data using t_star
+ # a_t = a*exp(-np.arange(inc)**2*pi**2*t_star/2)
+ # now apply the inverse discrete cosine transform
+ #density = idct(a_t)/R;
+
+ # take the rescaling of the data into account
+ bandwidth = sqrt(t_star)*R
+
+ # Kernel other than Gaussian scale bandwidth
+ h[dim] = bandwidth * (STEconstant / STEconstant2) ** (1.0 / 5)
+ #end % for dim loop
+ return h
+
def hstt(self, data, h0=None, inc=128, maxit=100, releps=0.01, abseps=0.0):
'''HSTT Scott-Tapia-Thompson estimate of smoothing parameter.
@@ -2461,7 +2557,7 @@ def kde_demo1():
pylab.axis([x.min(), x.max(), 0, 0.5])
def kde_demo2():
- '''Demonstrate the difference between transformation- and ordinary-KDE
+ '''Demonstrate the difference between transformation- and ordinary-KDE in 1D
KDEDEMO2 shows that the transformation KDE is a better estimate for
Rayleigh distributed data around 0 than the ordinary KDE.
@@ -2469,10 +2565,10 @@ def kde_demo2():
import scipy.stats as st
data = st.rayleigh.rvs(scale=1, size=300)
- x = np.linspace(1.5e-3, 5, 55);
+ x = np.linspace(1.5e-2, 5, 55);
kde = KDE(data)
- f = kde(output='plot', title='Ordinary KDE')
+ f = kde(output='plot', title='Ordinary KDE (hs=%g)' % kde.hs)
pylab.figure(0)
f.plot()
@@ -2481,7 +2577,7 @@ def kde_demo2():
#plotnorm((data).^(L2)) % gives a straight line => L2 = 0.5 reasonable
tkde = TKDE(data, L2=0.5)
- ft = tkde(x, output='plot', title='Transformation KDE')
+ ft = tkde(x, output='plot', title='Transformation KDE (hs=%g)' % tkde.tkde.hs)
pylab.figure(1)
ft.plot()
@@ -2490,7 +2586,7 @@ def kde_demo2():
pylab.figure(0)
def kde_demo3():
- '''Demonstrate the difference between and ordinary-KDE
+ '''Demonstrate the difference between transformation and ordinary-KDE in 2D
KDEDEMO3 shows that the transformation KDE is a better estimate for
Rayleigh distributed data around 0 than the ordinary KDE.
@@ -2521,6 +2617,9 @@ def kde_demo3():
pylab.figure(0)
def kde_demo4(hs=None, fast=False):
+ '''
+
+ '''
N = 100
#ei = np.random.normal(loc=0, scale=0.075, size=(N,))
ei = np.array([-0.08508516, 0.10462496, 0.07694448, -0.03080661, 0.05777525,
@@ -2568,6 +2667,36 @@ def kde_demo4(hs=None, fast=False):
print(kreg.tkde.tkde.inv_hs)
print(kreg.tkde.tkde.hs)
+
+def kde_demo5(N=50):
+ '''Demonstrate that the improved Sheather-Jones plug-in (hisj) is superior
+ for multimodal distributions
+
+ KDEDEMO5 shows that the improved Sheather-Jones plug-in smoothing is a better
+ compared to normal reference rules (in this case the hns)
+ '''
+ import scipy.stats as st
+
+ data = np.hstack((st.norm.rvs(loc=5, scale=1, size=(N,)), st.norm.rvs(loc=-5, scale=1, size=(N,))))
+
+ #x = np.linspace(1.5e-3, 5, 55)
+
+ kde = KDE(data, kernel=Kernel('gauss', 'hns'))
+ f = kde(output='plot', title='Ordinary KDE', plotflag=1)
+
+
+
+ kde1 = KDE(data, kernel=Kernel('gauss', 'hisj'))
+ f1 = kde1(output='plot', label='Ordinary KDE', plotflag=1)
+
+ pylab.figure(0)
+ f.plot('r', label='hns=%g' % kde.hs)
+ #pylab.figure(2)
+ f1.plot('b', label='hisj=%g' % kde1.hs)
+ x = np.linspace(-4,4)
+ for loc in [-5,5]:
+ pylab.plot(x + loc, st.norm.pdf(x, 0, scale=1)/2, 'k:', label='True density')
+ pylab.legend()
def test_docstrings():
import doctest
@@ -2575,4 +2704,4 @@ def test_docstrings():
if __name__ == '__main__':
#test_docstrings()
- kde_demo4()
\ No newline at end of file
+ kde_demo2()
\ No newline at end of file
diff --git a/pywafo/src/wafo/test/test_kdetools.py b/pywafo/src/wafo/test/test_kdetools.py
index 4a603d1..bdd5180 100644
--- a/pywafo/src/wafo/test/test_kdetools.py
+++ b/pywafo/src/wafo/test/test_kdetools.py
@@ -212,6 +212,16 @@ def test_smooth_params():
>>> gauss.hldpi(data)
array([ 0.1732289 , 0.33159097, 0.3107633 ])
+
+
+ >>> gauss.hisj(data)
+ array([ 0.29400043, 0.74277133, 0.51251583])
+
+ >>> data = np.array([0.753557920000000, 0.727791940000000, 0.941491690000000,
+ ... 0.078411190000000, 2.322918870000000, 1.104199950000000, 0.770551140000000,
+ ... 0.602882730000000, 1.368836350000000, 1.747543260000000, 1.095475610000000,
+ ... 1.016711330000000, 0.732111430000000, 0.618917190000000, 0.759034870000000,
+ ... 1.891946900000000, 0.724338080000000, 1.929730940000000, 0.447498380000000, 1.365084520000000])
'''
def test_gridcount_1D():
'''