realtransforms.py 23 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739
  1. """
  2. Real spectrum transforms (DCT, DST, MDCT)
  3. """
  4. from __future__ import division, print_function, absolute_import
  5. __all__ = ['dct', 'idct', 'dst', 'idst', 'dctn', 'idctn', 'dstn', 'idstn']
  6. import numpy as np
  7. from scipy.fftpack import _fftpack
  8. from scipy.fftpack.basic import _datacopied, _fix_shape, _asfarray
  9. from scipy.fftpack.helper import _init_nd_shape_and_axes
  10. import atexit
  11. atexit.register(_fftpack.destroy_ddct1_cache)
  12. atexit.register(_fftpack.destroy_ddct2_cache)
  13. atexit.register(_fftpack.destroy_ddct4_cache)
  14. atexit.register(_fftpack.destroy_dct1_cache)
  15. atexit.register(_fftpack.destroy_dct2_cache)
  16. atexit.register(_fftpack.destroy_dct4_cache)
  17. atexit.register(_fftpack.destroy_ddst1_cache)
  18. atexit.register(_fftpack.destroy_ddst2_cache)
  19. atexit.register(_fftpack.destroy_dst1_cache)
  20. atexit.register(_fftpack.destroy_dst2_cache)
  21. def dctn(x, type=2, shape=None, axes=None, norm=None, overwrite_x=False):
  22. """
  23. Return multidimensional Discrete Cosine Transform along the specified axes.
  24. Parameters
  25. ----------
  26. x : array_like
  27. The input array.
  28. type : {1, 2, 3, 4}, optional
  29. Type of the DCT (see Notes). Default type is 2.
  30. shape : int or array_like of ints or None, optional
  31. The shape of the result. If both `shape` and `axes` (see below) are
  32. None, `shape` is ``x.shape``; if `shape` is None but `axes` is
  33. not None, then `shape` is ``scipy.take(x.shape, axes, axis=0)``.
  34. If ``shape[i] > x.shape[i]``, the i-th dimension is padded with zeros.
  35. If ``shape[i] < x.shape[i]``, the i-th dimension is truncated to
  36. length ``shape[i]``.
  37. If any element of `shape` is -1, the size of the corresponding
  38. dimension of `x` is used.
  39. axes : int or array_like of ints or None, optional
  40. Axes along which the DCT is computed.
  41. The default is over all axes.
  42. norm : {None, 'ortho'}, optional
  43. Normalization mode (see Notes). Default is None.
  44. overwrite_x : bool, optional
  45. If True, the contents of `x` can be destroyed; the default is False.
  46. Returns
  47. -------
  48. y : ndarray of real
  49. The transformed input array.
  50. See Also
  51. --------
  52. idctn : Inverse multidimensional DCT
  53. Notes
  54. -----
  55. For full details of the DCT types and normalization modes, as well as
  56. references, see `dct`.
  57. Examples
  58. --------
  59. >>> from scipy.fftpack import dctn, idctn
  60. >>> y = np.random.randn(16, 16)
  61. >>> np.allclose(y, idctn(dctn(y, norm='ortho'), norm='ortho'))
  62. True
  63. """
  64. x = np.asanyarray(x)
  65. shape, axes = _init_nd_shape_and_axes(x, shape, axes)
  66. for n, ax in zip(shape, axes):
  67. x = dct(x, type=type, n=n, axis=ax, norm=norm, overwrite_x=overwrite_x)
  68. return x
  69. def idctn(x, type=2, shape=None, axes=None, norm=None, overwrite_x=False):
  70. """
  71. Return multidimensional Discrete Cosine Transform along the specified axes.
  72. Parameters
  73. ----------
  74. x : array_like
  75. The input array.
  76. type : {1, 2, 3, 4}, optional
  77. Type of the DCT (see Notes). Default type is 2.
  78. shape : int or array_like of ints or None, optional
  79. The shape of the result. If both `shape` and `axes` (see below) are
  80. None, `shape` is ``x.shape``; if `shape` is None but `axes` is
  81. not None, then `shape` is ``scipy.take(x.shape, axes, axis=0)``.
  82. If ``shape[i] > x.shape[i]``, the i-th dimension is padded with zeros.
  83. If ``shape[i] < x.shape[i]``, the i-th dimension is truncated to
  84. length ``shape[i]``.
  85. If any element of `shape` is -1, the size of the corresponding
  86. dimension of `x` is used.
  87. axes : int or array_like of ints or None, optional
  88. Axes along which the IDCT is computed.
  89. The default is over all axes.
  90. norm : {None, 'ortho'}, optional
  91. Normalization mode (see Notes). Default is None.
  92. overwrite_x : bool, optional
  93. If True, the contents of `x` can be destroyed; the default is False.
  94. Returns
  95. -------
  96. y : ndarray of real
  97. The transformed input array.
  98. See Also
  99. --------
  100. dctn : multidimensional DCT
  101. Notes
  102. -----
  103. For full details of the IDCT types and normalization modes, as well as
  104. references, see `idct`.
  105. Examples
  106. --------
  107. >>> from scipy.fftpack import dctn, idctn
  108. >>> y = np.random.randn(16, 16)
  109. >>> np.allclose(y, idctn(dctn(y, norm='ortho'), norm='ortho'))
  110. True
  111. """
  112. x = np.asanyarray(x)
  113. shape, axes = _init_nd_shape_and_axes(x, shape, axes)
  114. for n, ax in zip(shape, axes):
  115. x = idct(x, type=type, n=n, axis=ax, norm=norm,
  116. overwrite_x=overwrite_x)
  117. return x
  118. def dstn(x, type=2, shape=None, axes=None, norm=None, overwrite_x=False):
  119. """
  120. Return multidimensional Discrete Sine Transform along the specified axes.
  121. Parameters
  122. ----------
  123. x : array_like
  124. The input array.
  125. type : {1, 2, 3, 4}, optional
  126. Type of the DCT (see Notes). Default type is 2.
  127. shape : int or array_like of ints or None, optional
  128. The shape of the result. If both `shape` and `axes` (see below) are
  129. None, `shape` is ``x.shape``; if `shape` is None but `axes` is
  130. not None, then `shape` is ``scipy.take(x.shape, axes, axis=0)``.
  131. If ``shape[i] > x.shape[i]``, the i-th dimension is padded with zeros.
  132. If ``shape[i] < x.shape[i]``, the i-th dimension is truncated to
  133. length ``shape[i]``.
  134. If any element of `shape` is -1, the size of the corresponding
  135. dimension of `x` is used.
  136. axes : int or array_like of ints or None, optional
  137. Axes along which the DCT is computed.
  138. The default is over all axes.
  139. norm : {None, 'ortho'}, optional
  140. Normalization mode (see Notes). Default is None.
  141. overwrite_x : bool, optional
  142. If True, the contents of `x` can be destroyed; the default is False.
  143. Returns
  144. -------
  145. y : ndarray of real
  146. The transformed input array.
  147. See Also
  148. --------
  149. idstn : Inverse multidimensional DST
  150. Notes
  151. -----
  152. For full details of the DST types and normalization modes, as well as
  153. references, see `dst`.
  154. Examples
  155. --------
  156. >>> from scipy.fftpack import dstn, idstn
  157. >>> y = np.random.randn(16, 16)
  158. >>> np.allclose(y, idstn(dstn(y, norm='ortho'), norm='ortho'))
  159. True
  160. """
  161. x = np.asanyarray(x)
  162. shape, axes = _init_nd_shape_and_axes(x, shape, axes)
  163. for n, ax in zip(shape, axes):
  164. x = dst(x, type=type, n=n, axis=ax, norm=norm, overwrite_x=overwrite_x)
  165. return x
  166. def idstn(x, type=2, shape=None, axes=None, norm=None, overwrite_x=False):
  167. """
  168. Return multidimensional Discrete Sine Transform along the specified axes.
  169. Parameters
  170. ----------
  171. x : array_like
  172. The input array.
  173. type : {1, 2, 3, 4}, optional
  174. Type of the DCT (see Notes). Default type is 2.
  175. shape : int or array_like of ints or None, optional
  176. The shape of the result. If both `shape` and `axes` (see below) are
  177. None, `shape` is ``x.shape``; if `shape` is None but `axes` is
  178. not None, then `shape` is ``scipy.take(x.shape, axes, axis=0)``.
  179. If ``shape[i] > x.shape[i]``, the i-th dimension is padded with zeros.
  180. If ``shape[i] < x.shape[i]``, the i-th dimension is truncated to
  181. length ``shape[i]``.
  182. If any element of `shape` is -1, the size of the corresponding
  183. dimension of `x` is used.
  184. axes : int or array_like of ints or None, optional
  185. Axes along which the IDCT is computed.
  186. The default is over all axes.
  187. norm : {None, 'ortho'}, optional
  188. Normalization mode (see Notes). Default is None.
  189. overwrite_x : bool, optional
  190. If True, the contents of `x` can be destroyed; the default is False.
  191. Returns
  192. -------
  193. y : ndarray of real
  194. The transformed input array.
  195. See Also
  196. --------
  197. dctn : multidimensional DST
  198. Notes
  199. -----
  200. For full details of the IDST types and normalization modes, as well as
  201. references, see `idst`.
  202. Examples
  203. --------
  204. >>> from scipy.fftpack import dstn, idstn
  205. >>> y = np.random.randn(16, 16)
  206. >>> np.allclose(y, idstn(dstn(y, norm='ortho'), norm='ortho'))
  207. True
  208. """
  209. x = np.asanyarray(x)
  210. shape, axes = _init_nd_shape_and_axes(x, shape, axes)
  211. for n, ax in zip(shape, axes):
  212. x = idst(x, type=type, n=n, axis=ax, norm=norm,
  213. overwrite_x=overwrite_x)
  214. return x
  215. def dct(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False):
  216. """
  217. Return the Discrete Cosine Transform of arbitrary type sequence x.
  218. Parameters
  219. ----------
  220. x : array_like
  221. The input array.
  222. type : {1, 2, 3, 4}, optional
  223. Type of the DCT (see Notes). Default type is 2.
  224. n : int, optional
  225. Length of the transform. If ``n < x.shape[axis]``, `x` is
  226. truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
  227. default results in ``n = x.shape[axis]``.
  228. axis : int, optional
  229. Axis along which the dct is computed; the default is over the
  230. last axis (i.e., ``axis=-1``).
  231. norm : {None, 'ortho'}, optional
  232. Normalization mode (see Notes). Default is None.
  233. overwrite_x : bool, optional
  234. If True, the contents of `x` can be destroyed; the default is False.
  235. Returns
  236. -------
  237. y : ndarray of real
  238. The transformed input array.
  239. See Also
  240. --------
  241. idct : Inverse DCT
  242. Notes
  243. -----
  244. For a single dimension array ``x``, ``dct(x, norm='ortho')`` is equal to
  245. MATLAB ``dct(x)``.
  246. There are theoretically 8 types of the DCT, only the first 4 types are
  247. implemented in scipy. 'The' DCT generally refers to DCT type 2, and 'the'
  248. Inverse DCT generally refers to DCT type 3.
  249. **Type I**
  250. There are several definitions of the DCT-I; we use the following
  251. (for ``norm=None``)::
  252. N-2
  253. y[k] = x[0] + (-1)**k x[N-1] + 2 * sum x[n]*cos(pi*k*n/(N-1))
  254. n=1
  255. If ``norm='ortho'``, ``x[0]`` and ``x[N-1]`` are multiplied by
  256. a scaling factor of ``sqrt(2)``, and ``y[k]`` is multiplied by a
  257. scaling factor `f`::
  258. f = 0.5*sqrt(1/(N-1)) if k = 0 or N-1,
  259. f = 0.5*sqrt(2/(N-1)) otherwise.
  260. .. versionadded:: 1.2.0
  261. Orthonormalization in DCT-I.
  262. .. note::
  263. The DCT-I is only supported for input size > 1.
  264. **Type II**
  265. There are several definitions of the DCT-II; we use the following
  266. (for ``norm=None``)::
  267. N-1
  268. y[k] = 2* sum x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N.
  269. n=0
  270. If ``norm='ortho'``, ``y[k]`` is multiplied by a scaling factor `f`::
  271. f = sqrt(1/(4*N)) if k = 0,
  272. f = sqrt(1/(2*N)) otherwise.
  273. Which makes the corresponding matrix of coefficients orthonormal
  274. (``OO' = Id``).
  275. **Type III**
  276. There are several definitions, we use the following
  277. (for ``norm=None``)::
  278. N-1
  279. y[k] = x[0] + 2 * sum x[n]*cos(pi*(k+0.5)*n/N), 0 <= k < N.
  280. n=1
  281. or, for ``norm='ortho'`` and 0 <= k < N::
  282. N-1
  283. y[k] = x[0] / sqrt(N) + sqrt(2/N) * sum x[n]*cos(pi*(k+0.5)*n/N)
  284. n=1
  285. The (unnormalized) DCT-III is the inverse of the (unnormalized) DCT-II, up
  286. to a factor `2N`. The orthonormalized DCT-III is exactly the inverse of
  287. the orthonormalized DCT-II.
  288. **Type IV**
  289. There are several definitions of the DCT-IV; we use the following
  290. (for ``norm=None``)::
  291. N-1
  292. y[k] = 2* sum x[n]*cos(pi*(2k+1)*(2n+1)/(4*N)), 0 <= k < N.
  293. n=0
  294. If ``norm='ortho'``, ``y[k]`` is multiplied by a scaling factor `f`::
  295. f = 0.5*sqrt(2/N)
  296. .. versionadded:: 1.2.0
  297. Support for DCT-IV.
  298. References
  299. ----------
  300. .. [1] 'A Fast Cosine Transform in One and Two Dimensions', by J.
  301. Makhoul, `IEEE Transactions on acoustics, speech and signal
  302. processing` vol. 28(1), pp. 27-34,
  303. :doi:`10.1109/TASSP.1980.1163351` (1980).
  304. .. [2] Wikipedia, "Discrete cosine transform",
  305. https://en.wikipedia.org/wiki/Discrete_cosine_transform
  306. Examples
  307. --------
  308. The Type 1 DCT is equivalent to the FFT (though faster) for real,
  309. even-symmetrical inputs. The output is also real and even-symmetrical.
  310. Half of the FFT input is used to generate half of the FFT output:
  311. >>> from scipy.fftpack import fft, dct
  312. >>> fft(np.array([4., 3., 5., 10., 5., 3.])).real
  313. array([ 30., -8., 6., -2., 6., -8.])
  314. >>> dct(np.array([4., 3., 5., 10.]), 1)
  315. array([ 30., -8., 6., -2.])
  316. """
  317. return _dct(x, type, n, axis, normalize=norm, overwrite_x=overwrite_x)
  318. def idct(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False):
  319. """
  320. Return the Inverse Discrete Cosine Transform of an arbitrary type sequence.
  321. Parameters
  322. ----------
  323. x : array_like
  324. The input array.
  325. type : {1, 2, 3, 4}, optional
  326. Type of the DCT (see Notes). Default type is 2.
  327. n : int, optional
  328. Length of the transform. If ``n < x.shape[axis]``, `x` is
  329. truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
  330. default results in ``n = x.shape[axis]``.
  331. axis : int, optional
  332. Axis along which the idct is computed; the default is over the
  333. last axis (i.e., ``axis=-1``).
  334. norm : {None, 'ortho'}, optional
  335. Normalization mode (see Notes). Default is None.
  336. overwrite_x : bool, optional
  337. If True, the contents of `x` can be destroyed; the default is False.
  338. Returns
  339. -------
  340. idct : ndarray of real
  341. The transformed input array.
  342. See Also
  343. --------
  344. dct : Forward DCT
  345. Notes
  346. -----
  347. For a single dimension array `x`, ``idct(x, norm='ortho')`` is equal to
  348. MATLAB ``idct(x)``.
  349. 'The' IDCT is the IDCT of type 2, which is the same as DCT of type 3.
  350. IDCT of type 1 is the DCT of type 1, IDCT of type 2 is the DCT of type
  351. 3, and IDCT of type 3 is the DCT of type 2. IDCT of type 4 is the DCT
  352. of type 4. For the definition of these types, see `dct`.
  353. Examples
  354. --------
  355. The Type 1 DCT is equivalent to the DFT for real, even-symmetrical
  356. inputs. The output is also real and even-symmetrical. Half of the IFFT
  357. input is used to generate half of the IFFT output:
  358. >>> from scipy.fftpack import ifft, idct
  359. >>> ifft(np.array([ 30., -8., 6., -2., 6., -8.])).real
  360. array([ 4., 3., 5., 10., 5., 3.])
  361. >>> idct(np.array([ 30., -8., 6., -2.]), 1) / 6
  362. array([ 4., 3., 5., 10.])
  363. """
  364. # Inverse/forward type table
  365. _TP = {1:1, 2:3, 3:2, 4:4}
  366. return _dct(x, _TP[type], n, axis, normalize=norm, overwrite_x=overwrite_x)
  367. def _get_dct_fun(type, dtype):
  368. try:
  369. name = {'float64':'ddct%d', 'float32':'dct%d'}[dtype.name]
  370. except KeyError:
  371. raise ValueError("dtype %s not supported" % dtype)
  372. try:
  373. f = getattr(_fftpack, name % type)
  374. except AttributeError as e:
  375. raise ValueError(str(e) + ". Type %d not understood" % type)
  376. return f
  377. def _get_norm_mode(normalize):
  378. try:
  379. nm = {None:0, 'ortho':1}[normalize]
  380. except KeyError:
  381. raise ValueError("Unknown normalize mode %s" % normalize)
  382. return nm
  383. def __fix_shape(x, n, axis, dct_or_dst):
  384. tmp = _asfarray(x)
  385. copy_made = _datacopied(tmp, x)
  386. if n is None:
  387. n = tmp.shape[axis]
  388. elif n != tmp.shape[axis]:
  389. tmp, copy_made2 = _fix_shape(tmp, n, axis)
  390. copy_made = copy_made or copy_made2
  391. if n < 1:
  392. raise ValueError("Invalid number of %s data points "
  393. "(%d) specified." % (dct_or_dst, n))
  394. return tmp, n, copy_made
  395. def _raw_dct(x0, type, n, axis, nm, overwrite_x):
  396. f = _get_dct_fun(type, x0.dtype)
  397. return _eval_fun(f, x0, n, axis, nm, overwrite_x)
  398. def _raw_dst(x0, type, n, axis, nm, overwrite_x):
  399. f = _get_dst_fun(type, x0.dtype)
  400. return _eval_fun(f, x0, n, axis, nm, overwrite_x)
  401. def _eval_fun(f, tmp, n, axis, nm, overwrite_x):
  402. if axis == -1 or axis == len(tmp.shape) - 1:
  403. return f(tmp, n, nm, overwrite_x)
  404. tmp = np.swapaxes(tmp, axis, -1)
  405. tmp = f(tmp, n, nm, overwrite_x)
  406. return np.swapaxes(tmp, axis, -1)
  407. def _dct(x, type, n=None, axis=-1, overwrite_x=False, normalize=None):
  408. """
  409. Return Discrete Cosine Transform of arbitrary type sequence x.
  410. Parameters
  411. ----------
  412. x : array_like
  413. input array.
  414. n : int, optional
  415. Length of the transform. If ``n < x.shape[axis]``, `x` is
  416. truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
  417. default results in ``n = x.shape[axis]``.
  418. axis : int, optional
  419. Axis along which the dct is computed; the default is over the
  420. last axis (i.e., ``axis=-1``).
  421. overwrite_x : bool, optional
  422. If True, the contents of `x` can be destroyed; the default is False.
  423. Returns
  424. -------
  425. z : ndarray
  426. """
  427. x0, n, copy_made = __fix_shape(x, n, axis, 'DCT')
  428. if type == 1 and n < 2:
  429. raise ValueError("DCT-I is not defined for size < 2")
  430. overwrite_x = overwrite_x or copy_made
  431. nm = _get_norm_mode(normalize)
  432. if np.iscomplexobj(x0):
  433. return (_raw_dct(x0.real, type, n, axis, nm, overwrite_x) + 1j *
  434. _raw_dct(x0.imag, type, n, axis, nm, overwrite_x))
  435. else:
  436. return _raw_dct(x0, type, n, axis, nm, overwrite_x)
  437. def dst(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False):
  438. """
  439. Return the Discrete Sine Transform of arbitrary type sequence x.
  440. Parameters
  441. ----------
  442. x : array_like
  443. The input array.
  444. type : {1, 2, 3, 4}, optional
  445. Type of the DST (see Notes). Default type is 2.
  446. n : int, optional
  447. Length of the transform. If ``n < x.shape[axis]``, `x` is
  448. truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
  449. default results in ``n = x.shape[axis]``.
  450. axis : int, optional
  451. Axis along which the dst is computed; the default is over the
  452. last axis (i.e., ``axis=-1``).
  453. norm : {None, 'ortho'}, optional
  454. Normalization mode (see Notes). Default is None.
  455. overwrite_x : bool, optional
  456. If True, the contents of `x` can be destroyed; the default is False.
  457. Returns
  458. -------
  459. dst : ndarray of reals
  460. The transformed input array.
  461. See Also
  462. --------
  463. idst : Inverse DST
  464. Notes
  465. -----
  466. For a single dimension array ``x``.
  467. There are theoretically 8 types of the DST for different combinations of
  468. even/odd boundary conditions and boundary off sets [1]_, only the first
  469. 3 types are implemented in scipy.
  470. **Type I**
  471. There are several definitions of the DST-I; we use the following
  472. for ``norm=None``. DST-I assumes the input is odd around n=-1 and n=N. ::
  473. N-1
  474. y[k] = 2 * sum x[n]*sin(pi*(k+1)*(n+1)/(N+1))
  475. n=0
  476. Note that the DST-I is only supported for input size > 1
  477. The (unnormalized) DST-I is its own inverse, up to a factor `2(N+1)`.
  478. The orthonormalized DST-I is exactly its own inverse.
  479. **Type II**
  480. There are several definitions of the DST-II; we use the following
  481. for ``norm=None``. DST-II assumes the input is odd around n=-1/2 and
  482. n=N-1/2; the output is odd around k=-1 and even around k=N-1 ::
  483. N-1
  484. y[k] = 2* sum x[n]*sin(pi*(k+1)*(n+0.5)/N), 0 <= k < N.
  485. n=0
  486. if ``norm='ortho'``, ``y[k]`` is multiplied by a scaling factor `f` ::
  487. f = sqrt(1/(4*N)) if k == 0
  488. f = sqrt(1/(2*N)) otherwise.
  489. **Type III**
  490. There are several definitions of the DST-III, we use the following
  491. (for ``norm=None``). DST-III assumes the input is odd around n=-1
  492. and even around n=N-1 ::
  493. N-2
  494. y[k] = x[N-1]*(-1)**k + 2* sum x[n]*sin(pi*(k+0.5)*(n+1)/N), 0 <= k < N.
  495. n=0
  496. The (unnormalized) DST-III is the inverse of the (unnormalized) DST-II, up
  497. to a factor `2N`. The orthonormalized DST-III is exactly the inverse of
  498. the orthonormalized DST-II.
  499. .. versionadded:: 0.11.0
  500. **Type IV**
  501. There are several definitions of the DST-IV, we use the following
  502. (for ``norm=None``). DST-IV assumes the input is odd around n=-0.5
  503. and even around n=N-0.5 ::
  504. N-1
  505. y[k] = 2* sum x[n]*sin(pi*(k+0.5)*(n+0.5)/N), 0 <= k < N.
  506. n=0
  507. The (unnormalized) DST-IV is its own inverse, up
  508. to a factor `2N`. The orthonormalized DST-IV is exactly its own inverse.
  509. .. versionadded:: 1.2.0
  510. Support for DST-IV.
  511. References
  512. ----------
  513. .. [1] Wikipedia, "Discrete sine transform",
  514. https://en.wikipedia.org/wiki/Discrete_sine_transform
  515. """
  516. return _dst(x, type, n, axis, normalize=norm, overwrite_x=overwrite_x)
  517. def idst(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False):
  518. """
  519. Return the Inverse Discrete Sine Transform of an arbitrary type sequence.
  520. Parameters
  521. ----------
  522. x : array_like
  523. The input array.
  524. type : {1, 2, 3, 4}, optional
  525. Type of the DST (see Notes). Default type is 2.
  526. n : int, optional
  527. Length of the transform. If ``n < x.shape[axis]``, `x` is
  528. truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
  529. default results in ``n = x.shape[axis]``.
  530. axis : int, optional
  531. Axis along which the idst is computed; the default is over the
  532. last axis (i.e., ``axis=-1``).
  533. norm : {None, 'ortho'}, optional
  534. Normalization mode (see Notes). Default is None.
  535. overwrite_x : bool, optional
  536. If True, the contents of `x` can be destroyed; the default is False.
  537. Returns
  538. -------
  539. idst : ndarray of real
  540. The transformed input array.
  541. See Also
  542. --------
  543. dst : Forward DST
  544. Notes
  545. -----
  546. 'The' IDST is the IDST of type 2, which is the same as DST of type 3.
  547. IDST of type 1 is the DST of type 1, IDST of type 2 is the DST of type
  548. 3, and IDST of type 3 is the DST of type 2. For the definition of these
  549. types, see `dst`.
  550. .. versionadded:: 0.11.0
  551. """
  552. # Inverse/forward type table
  553. _TP = {1:1, 2:3, 3:2, 4:4}
  554. return _dst(x, _TP[type], n, axis, normalize=norm, overwrite_x=overwrite_x)
  555. def _get_dst_fun(type, dtype):
  556. try:
  557. name = {'float64':'ddst%d', 'float32':'dst%d'}[dtype.name]
  558. except KeyError:
  559. raise ValueError("dtype %s not supported" % dtype)
  560. try:
  561. f = getattr(_fftpack, name % type)
  562. except AttributeError as e:
  563. raise ValueError(str(e) + ". Type %d not understood" % type)
  564. return f
  565. def _dst(x, type, n=None, axis=-1, overwrite_x=False, normalize=None):
  566. """
  567. Return Discrete Sine Transform of arbitrary type sequence x.
  568. Parameters
  569. ----------
  570. x : array_like
  571. input array.
  572. n : int, optional
  573. Length of the transform.
  574. axis : int, optional
  575. Axis along which the dst is computed. (default=-1)
  576. overwrite_x : bool, optional
  577. If True the contents of x can be destroyed. (default=False)
  578. Returns
  579. -------
  580. z : real ndarray
  581. """
  582. x0, n, copy_made = __fix_shape(x, n, axis, 'DST')
  583. if type == 1 and n < 2:
  584. raise ValueError("DST-I is not defined for size < 2")
  585. overwrite_x = overwrite_x or copy_made
  586. nm = _get_norm_mode(normalize)
  587. if np.iscomplexobj(x0):
  588. return (_raw_dst(x0.real, type, n, axis, nm, overwrite_x) + 1j *
  589. _raw_dst(x0.imag, type, n, axis, nm, overwrite_x))
  590. else:
  591. return _raw_dst(x0, type, n, axis, nm, overwrite_x)