_realtransforms.py 19 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598
  1. """
  2. Real spectrum transforms (DCT, DST, MDCT)
  3. """
  4. __all__ = ['dct', 'idct', 'dst', 'idst', 'dctn', 'idctn', 'dstn', 'idstn']
  5. from scipy.fft import _pocketfft
  6. from ._helper import _good_shape
  7. _inverse_typemap = {1: 1, 2: 3, 3: 2, 4: 4}
  8. def dctn(x, type=2, shape=None, axes=None, norm=None, overwrite_x=False):
  9. """
  10. Return multidimensional Discrete Cosine Transform along the specified axes.
  11. Parameters
  12. ----------
  13. x : array_like
  14. The input array.
  15. type : {1, 2, 3, 4}, optional
  16. Type of the DCT (see Notes). Default type is 2.
  17. shape : int or array_like of ints or None, optional
  18. The shape of the result. If both `shape` and `axes` (see below) are
  19. None, `shape` is ``x.shape``; if `shape` is None but `axes` is
  20. not None, then `shape` is ``numpy.take(x.shape, axes, axis=0)``.
  21. If ``shape[i] > x.shape[i]``, the ith dimension is padded with zeros.
  22. If ``shape[i] < x.shape[i]``, the ith dimension is truncated to
  23. length ``shape[i]``.
  24. If any element of `shape` is -1, the size of the corresponding
  25. dimension of `x` is used.
  26. axes : int or array_like of ints or None, optional
  27. Axes along which the DCT is computed.
  28. The default is over all axes.
  29. norm : {None, 'ortho'}, optional
  30. Normalization mode (see Notes). Default is None.
  31. overwrite_x : bool, optional
  32. If True, the contents of `x` can be destroyed; the default is False.
  33. Returns
  34. -------
  35. y : ndarray of real
  36. The transformed input array.
  37. See Also
  38. --------
  39. idctn : Inverse multidimensional DCT
  40. Notes
  41. -----
  42. For full details of the DCT types and normalization modes, as well as
  43. references, see `dct`.
  44. Examples
  45. --------
  46. >>> import numpy as np
  47. >>> from scipy.fftpack import dctn, idctn
  48. >>> rng = np.random.default_rng()
  49. >>> y = rng.standard_normal((16, 16))
  50. >>> np.allclose(y, idctn(dctn(y, norm='ortho'), norm='ortho'))
  51. True
  52. """
  53. shape = _good_shape(x, shape, axes)
  54. return _pocketfft.dctn(x, type, shape, axes, norm, overwrite_x)
  55. def idctn(x, type=2, shape=None, axes=None, norm=None, overwrite_x=False):
  56. """
  57. Return multidimensional Discrete Cosine Transform along the specified axes.
  58. Parameters
  59. ----------
  60. x : array_like
  61. The input array.
  62. type : {1, 2, 3, 4}, optional
  63. Type of the DCT (see Notes). Default type is 2.
  64. shape : int or array_like of ints or None, optional
  65. The shape of the result. If both `shape` and `axes` (see below) are
  66. None, `shape` is ``x.shape``; if `shape` is None but `axes` is
  67. not None, then `shape` is ``numpy.take(x.shape, axes, axis=0)``.
  68. If ``shape[i] > x.shape[i]``, the ith dimension is padded with zeros.
  69. If ``shape[i] < x.shape[i]``, the ith dimension is truncated to
  70. length ``shape[i]``.
  71. If any element of `shape` is -1, the size of the corresponding
  72. dimension of `x` is used.
  73. axes : int or array_like of ints or None, optional
  74. Axes along which the IDCT is computed.
  75. The default is over all axes.
  76. norm : {None, 'ortho'}, optional
  77. Normalization mode (see Notes). Default is None.
  78. overwrite_x : bool, optional
  79. If True, the contents of `x` can be destroyed; the default is False.
  80. Returns
  81. -------
  82. y : ndarray of real
  83. The transformed input array.
  84. See Also
  85. --------
  86. dctn : multidimensional DCT
  87. Notes
  88. -----
  89. For full details of the IDCT types and normalization modes, as well as
  90. references, see `idct`.
  91. Examples
  92. --------
  93. >>> import numpy as np
  94. >>> from scipy.fftpack import dctn, idctn
  95. >>> rng = np.random.default_rng()
  96. >>> y = rng.standard_normal((16, 16))
  97. >>> np.allclose(y, idctn(dctn(y, norm='ortho'), norm='ortho'))
  98. True
  99. """
  100. type = _inverse_typemap[type]
  101. shape = _good_shape(x, shape, axes)
  102. return _pocketfft.dctn(x, type, shape, axes, norm, overwrite_x)
  103. def dstn(x, type=2, shape=None, axes=None, norm=None, overwrite_x=False):
  104. """
  105. Return multidimensional Discrete Sine Transform along the specified axes.
  106. Parameters
  107. ----------
  108. x : array_like
  109. The input array.
  110. type : {1, 2, 3, 4}, optional
  111. Type of the DST (see Notes). Default type is 2.
  112. shape : int or array_like of ints or None, optional
  113. The shape of the result. If both `shape` and `axes` (see below) are
  114. None, `shape` is ``x.shape``; if `shape` is None but `axes` is
  115. not None, then `shape` is ``numpy.take(x.shape, axes, axis=0)``.
  116. If ``shape[i] > x.shape[i]``, the ith dimension is padded with zeros.
  117. If ``shape[i] < x.shape[i]``, the ith dimension is truncated to
  118. length ``shape[i]``.
  119. If any element of `shape` is -1, the size of the corresponding
  120. dimension of `x` is used.
  121. axes : int or array_like of ints or None, optional
  122. Axes along which the DCT is computed.
  123. The default is over all axes.
  124. norm : {None, 'ortho'}, optional
  125. Normalization mode (see Notes). Default is None.
  126. overwrite_x : bool, optional
  127. If True, the contents of `x` can be destroyed; the default is False.
  128. Returns
  129. -------
  130. y : ndarray of real
  131. The transformed input array.
  132. See Also
  133. --------
  134. idstn : Inverse multidimensional DST
  135. Notes
  136. -----
  137. For full details of the DST types and normalization modes, as well as
  138. references, see `dst`.
  139. Examples
  140. --------
  141. >>> import numpy as np
  142. >>> from scipy.fftpack import dstn, idstn
  143. >>> rng = np.random.default_rng()
  144. >>> y = rng.standard_normal((16, 16))
  145. >>> np.allclose(y, idstn(dstn(y, norm='ortho'), norm='ortho'))
  146. True
  147. """
  148. shape = _good_shape(x, shape, axes)
  149. return _pocketfft.dstn(x, type, shape, axes, norm, overwrite_x)
  150. def idstn(x, type=2, shape=None, axes=None, norm=None, overwrite_x=False):
  151. """
  152. Return multidimensional Discrete Sine Transform along the specified axes.
  153. Parameters
  154. ----------
  155. x : array_like
  156. The input array.
  157. type : {1, 2, 3, 4}, optional
  158. Type of the DST (see Notes). Default type is 2.
  159. shape : int or array_like of ints or None, optional
  160. The shape of the result. If both `shape` and `axes` (see below) are
  161. None, `shape` is ``x.shape``; if `shape` is None but `axes` is
  162. not None, then `shape` is ``numpy.take(x.shape, axes, axis=0)``.
  163. If ``shape[i] > x.shape[i]``, the ith dimension is padded with zeros.
  164. If ``shape[i] < x.shape[i]``, the ith dimension is truncated to
  165. length ``shape[i]``.
  166. If any element of `shape` is -1, the size of the corresponding
  167. dimension of `x` is used.
  168. axes : int or array_like of ints or None, optional
  169. Axes along which the IDST is computed.
  170. The default is over all axes.
  171. norm : {None, 'ortho'}, optional
  172. Normalization mode (see Notes). Default is None.
  173. overwrite_x : bool, optional
  174. If True, the contents of `x` can be destroyed; the default is False.
  175. Returns
  176. -------
  177. y : ndarray of real
  178. The transformed input array.
  179. See Also
  180. --------
  181. dstn : multidimensional DST
  182. Notes
  183. -----
  184. For full details of the IDST types and normalization modes, as well as
  185. references, see `idst`.
  186. Examples
  187. --------
  188. >>> import numpy as np
  189. >>> from scipy.fftpack import dstn, idstn
  190. >>> rng = np.random.default_rng()
  191. >>> y = rng.standard_normal((16, 16))
  192. >>> np.allclose(y, idstn(dstn(y, norm='ortho'), norm='ortho'))
  193. True
  194. """
  195. type = _inverse_typemap[type]
  196. shape = _good_shape(x, shape, axes)
  197. return _pocketfft.dstn(x, type, shape, axes, norm, overwrite_x)
  198. def dct(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False):
  199. r"""
  200. Return the Discrete Cosine Transform of arbitrary type sequence x.
  201. Parameters
  202. ----------
  203. x : array_like
  204. The input array.
  205. type : {1, 2, 3, 4}, optional
  206. Type of the DCT (see Notes). Default type is 2.
  207. n : int, optional
  208. Length of the transform. If ``n < x.shape[axis]``, `x` is
  209. truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
  210. default results in ``n = x.shape[axis]``.
  211. axis : int, optional
  212. Axis along which the dct is computed; the default is over the
  213. last axis (i.e., ``axis=-1``).
  214. norm : {None, 'ortho'}, optional
  215. Normalization mode (see Notes). Default is None.
  216. overwrite_x : bool, optional
  217. If True, the contents of `x` can be destroyed; the default is False.
  218. Returns
  219. -------
  220. y : ndarray of real
  221. The transformed input array.
  222. See Also
  223. --------
  224. idct : Inverse DCT
  225. Notes
  226. -----
  227. For a single dimension array ``x``, ``dct(x, norm='ortho')`` is equal to
  228. MATLAB ``dct(x)``.
  229. There are, theoretically, 8 types of the DCT, only the first 4 types are
  230. implemented in scipy. 'The' DCT generally refers to DCT type 2, and 'the'
  231. Inverse DCT generally refers to DCT type 3.
  232. **Type I**
  233. There are several definitions of the DCT-I; we use the following
  234. (for ``norm=None``)
  235. .. math::
  236. y_k = x_0 + (-1)^k x_{N-1} + 2 \sum_{n=1}^{N-2} x_n \cos\left(
  237. \frac{\pi k n}{N-1} \right)
  238. If ``norm='ortho'``, ``x[0]`` and ``x[N-1]`` are multiplied by a scaling
  239. factor of :math:`\sqrt{2}`, and ``y[k]`` is multiplied by a scaling factor
  240. ``f``
  241. .. math::
  242. f = \begin{cases}
  243. \frac{1}{2}\sqrt{\frac{1}{N-1}} & \text{if }k=0\text{ or }N-1, \\
  244. \frac{1}{2}\sqrt{\frac{2}{N-1}} & \text{otherwise} \end{cases}
  245. .. versionadded:: 1.2.0
  246. Orthonormalization in DCT-I.
  247. .. note::
  248. The DCT-I is only supported for input size > 1.
  249. **Type II**
  250. There are several definitions of the DCT-II; we use the following
  251. (for ``norm=None``)
  252. .. math::
  253. y_k = 2 \sum_{n=0}^{N-1} x_n \cos\left(\frac{\pi k(2n+1)}{2N} \right)
  254. If ``norm='ortho'``, ``y[k]`` is multiplied by a scaling factor ``f``
  255. .. math::
  256. f = \begin{cases}
  257. \sqrt{\frac{1}{4N}} & \text{if }k=0, \\
  258. \sqrt{\frac{1}{2N}} & \text{otherwise} \end{cases}
  259. which makes the corresponding matrix of coefficients orthonormal
  260. (``O @ O.T = np.eye(N)``).
  261. **Type III**
  262. There are several definitions, we use the following (for ``norm=None``)
  263. .. math::
  264. y_k = x_0 + 2 \sum_{n=1}^{N-1} x_n \cos\left(\frac{\pi(2k+1)n}{2N}\right)
  265. or, for ``norm='ortho'``
  266. .. math::
  267. y_k = \frac{x_0}{\sqrt{N}} + \sqrt{\frac{2}{N}} \sum_{n=1}^{N-1} x_n
  268. \cos\left(\frac{\pi(2k+1)n}{2N}\right)
  269. The (unnormalized) DCT-III is the inverse of the (unnormalized) DCT-II, up
  270. to a factor `2N`. The orthonormalized DCT-III is exactly the inverse of
  271. the orthonormalized DCT-II.
  272. **Type IV**
  273. There are several definitions of the DCT-IV; we use the following
  274. (for ``norm=None``)
  275. .. math::
  276. y_k = 2 \sum_{n=0}^{N-1} x_n \cos\left(\frac{\pi(2k+1)(2n+1)}{4N} \right)
  277. If ``norm='ortho'``, ``y[k]`` is multiplied by a scaling factor ``f``
  278. .. math::
  279. f = \frac{1}{\sqrt{2N}}
  280. .. versionadded:: 1.2.0
  281. Support for DCT-IV.
  282. References
  283. ----------
  284. .. [1] 'A Fast Cosine Transform in One and Two Dimensions', by J.
  285. Makhoul, `IEEE Transactions on acoustics, speech and signal
  286. processing` vol. 28(1), pp. 27-34,
  287. :doi:`10.1109/TASSP.1980.1163351` (1980).
  288. .. [2] Wikipedia, "Discrete cosine transform",
  289. https://en.wikipedia.org/wiki/Discrete_cosine_transform
  290. Examples
  291. --------
  292. The Type 1 DCT is equivalent to the FFT (though faster) for real,
  293. even-symmetrical inputs. The output is also real and even-symmetrical.
  294. Half of the FFT input is used to generate half of the FFT output:
  295. >>> from scipy.fftpack import fft, dct
  296. >>> import numpy as np
  297. >>> fft(np.array([4., 3., 5., 10., 5., 3.])).real
  298. array([ 30., -8., 6., -2., 6., -8.])
  299. >>> dct(np.array([4., 3., 5., 10.]), 1)
  300. array([ 30., -8., 6., -2.])
  301. """
  302. return _pocketfft.dct(x, type, n, axis, norm, overwrite_x)
  303. def idct(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False):
  304. """
  305. Return the Inverse Discrete Cosine Transform of an arbitrary type sequence.
  306. Parameters
  307. ----------
  308. x : array_like
  309. The input array.
  310. type : {1, 2, 3, 4}, optional
  311. Type of the DCT (see Notes). Default type is 2.
  312. n : int, optional
  313. Length of the transform. If ``n < x.shape[axis]``, `x` is
  314. truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
  315. default results in ``n = x.shape[axis]``.
  316. axis : int, optional
  317. Axis along which the idct is computed; the default is over the
  318. last axis (i.e., ``axis=-1``).
  319. norm : {None, 'ortho'}, optional
  320. Normalization mode (see Notes). Default is None.
  321. overwrite_x : bool, optional
  322. If True, the contents of `x` can be destroyed; the default is False.
  323. Returns
  324. -------
  325. idct : ndarray of real
  326. The transformed input array.
  327. See Also
  328. --------
  329. dct : Forward DCT
  330. Notes
  331. -----
  332. For a single dimension array `x`, ``idct(x, norm='ortho')`` is equal to
  333. MATLAB ``idct(x)``.
  334. 'The' IDCT is the IDCT of type 2, which is the same as DCT of type 3.
  335. IDCT of type 1 is the DCT of type 1, IDCT of type 2 is the DCT of type
  336. 3, and IDCT of type 3 is the DCT of type 2. IDCT of type 4 is the DCT
  337. of type 4. For the definition of these types, see `dct`.
  338. Examples
  339. --------
  340. The Type 1 DCT is equivalent to the DFT for real, even-symmetrical
  341. inputs. The output is also real and even-symmetrical. Half of the IFFT
  342. input is used to generate half of the IFFT output:
  343. >>> from scipy.fftpack import ifft, idct
  344. >>> import numpy as np
  345. >>> ifft(np.array([ 30., -8., 6., -2., 6., -8.])).real
  346. array([ 4., 3., 5., 10., 5., 3.])
  347. >>> idct(np.array([ 30., -8., 6., -2.]), 1) / 6
  348. array([ 4., 3., 5., 10.])
  349. """
  350. type = _inverse_typemap[type]
  351. return _pocketfft.dct(x, type, n, axis, norm, overwrite_x)
  352. def dst(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False):
  353. r"""
  354. Return the Discrete Sine Transform of arbitrary type sequence x.
  355. Parameters
  356. ----------
  357. x : array_like
  358. The input array.
  359. type : {1, 2, 3, 4}, optional
  360. Type of the DST (see Notes). Default type is 2.
  361. n : int, optional
  362. Length of the transform. If ``n < x.shape[axis]``, `x` is
  363. truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
  364. default results in ``n = x.shape[axis]``.
  365. axis : int, optional
  366. Axis along which the dst is computed; the default is over the
  367. last axis (i.e., ``axis=-1``).
  368. norm : {None, 'ortho'}, optional
  369. Normalization mode (see Notes). Default is None.
  370. overwrite_x : bool, optional
  371. If True, the contents of `x` can be destroyed; the default is False.
  372. Returns
  373. -------
  374. dst : ndarray of reals
  375. The transformed input array.
  376. See Also
  377. --------
  378. idst : Inverse DST
  379. Notes
  380. -----
  381. For a single dimension array ``x``.
  382. There are, theoretically, 8 types of the DST for different combinations of
  383. even/odd boundary conditions and boundary off sets [1]_, only the first
  384. 4 types are implemented in scipy.
  385. **Type I**
  386. There are several definitions of the DST-I; we use the following
  387. for ``norm=None``. DST-I assumes the input is odd around `n=-1` and `n=N`.
  388. .. math::
  389. y_k = 2 \sum_{n=0}^{N-1} x_n \sin\left(\frac{\pi(k+1)(n+1)}{N+1}\right)
  390. Note that the DST-I is only supported for input size > 1.
  391. The (unnormalized) DST-I is its own inverse, up to a factor `2(N+1)`.
  392. The orthonormalized DST-I is exactly its own inverse.
  393. **Type II**
  394. There are several definitions of the DST-II; we use the following for
  395. ``norm=None``. DST-II assumes the input is odd around `n=-1/2` and
  396. `n=N-1/2`; the output is odd around :math:`k=-1` and even around `k=N-1`
  397. .. math::
  398. y_k = 2 \sum_{n=0}^{N-1} x_n \sin\left(\frac{\pi(k+1)(2n+1)}{2N}\right)
  399. if ``norm='ortho'``, ``y[k]`` is multiplied by a scaling factor ``f``
  400. .. math::
  401. f = \begin{cases}
  402. \sqrt{\frac{1}{4N}} & \text{if }k = 0, \\
  403. \sqrt{\frac{1}{2N}} & \text{otherwise} \end{cases}
  404. **Type III**
  405. There are several definitions of the DST-III, we use the following (for
  406. ``norm=None``). DST-III assumes the input is odd around `n=-1` and even
  407. around `n=N-1`
  408. .. math::
  409. y_k = (-1)^k x_{N-1} + 2 \sum_{n=0}^{N-2} x_n \sin\left(
  410. \frac{\pi(2k+1)(n+1)}{2N}\right)
  411. The (unnormalized) DST-III is the inverse of the (unnormalized) DST-II, up
  412. to a factor `2N`. The orthonormalized DST-III is exactly the inverse of the
  413. orthonormalized DST-II.
  414. .. versionadded:: 0.11.0
  415. **Type IV**
  416. There are several definitions of the DST-IV, we use the following (for
  417. ``norm=None``). DST-IV assumes the input is odd around `n=-0.5` and even
  418. around `n=N-0.5`
  419. .. math::
  420. y_k = 2 \sum_{n=0}^{N-1} x_n \sin\left(\frac{\pi(2k+1)(2n+1)}{4N}\right)
  421. The (unnormalized) DST-IV is its own inverse, up to a factor `2N`. The
  422. orthonormalized DST-IV is exactly its own inverse.
  423. .. versionadded:: 1.2.0
  424. Support for DST-IV.
  425. References
  426. ----------
  427. .. [1] Wikipedia, "Discrete sine transform",
  428. https://en.wikipedia.org/wiki/Discrete_sine_transform
  429. """
  430. return _pocketfft.dst(x, type, n, axis, norm, overwrite_x)
  431. def idst(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False):
  432. """
  433. Return the Inverse Discrete Sine Transform of an arbitrary type sequence.
  434. Parameters
  435. ----------
  436. x : array_like
  437. The input array.
  438. type : {1, 2, 3, 4}, optional
  439. Type of the DST (see Notes). Default type is 2.
  440. n : int, optional
  441. Length of the transform. If ``n < x.shape[axis]``, `x` is
  442. truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
  443. default results in ``n = x.shape[axis]``.
  444. axis : int, optional
  445. Axis along which the idst is computed; the default is over the
  446. last axis (i.e., ``axis=-1``).
  447. norm : {None, 'ortho'}, optional
  448. Normalization mode (see Notes). Default is None.
  449. overwrite_x : bool, optional
  450. If True, the contents of `x` can be destroyed; the default is False.
  451. Returns
  452. -------
  453. idst : ndarray of real
  454. The transformed input array.
  455. See Also
  456. --------
  457. dst : Forward DST
  458. Notes
  459. -----
  460. 'The' IDST is the IDST of type 2, which is the same as DST of type 3.
  461. IDST of type 1 is the DST of type 1, IDST of type 2 is the DST of type
  462. 3, and IDST of type 3 is the DST of type 2. For the definition of these
  463. types, see `dst`.
  464. .. versionadded:: 0.11.0
  465. """
  466. type = _inverse_typemap[type]
  467. return _pocketfft.dst(x, type, n, axis, norm, overwrite_x)