_realtransforms.py 26 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717
  1. from ._basic import _dispatch
  2. from scipy._lib.uarray import Dispatchable
  3. from scipy._lib._array_api import xp_capabilities
  4. import numpy as np
  5. __all__ = ['dct', 'idct', 'dst', 'idst', 'dctn', 'idctn', 'dstn', 'idstn']
  6. @xp_capabilities(cpu_only=True, allow_dask_compute=True)
  7. @_dispatch
  8. def dctn(x, type=2, s=None, axes=None, norm=None, overwrite_x=False,
  9. workers=None, *, orthogonalize=None):
  10. """
  11. Return multidimensional Discrete Cosine Transform along the specified axes.
  12. Parameters
  13. ----------
  14. x : array_like
  15. The input array.
  16. type : {1, 2, 3, 4}, optional
  17. Type of the DCT (see Notes). Default type is 2.
  18. s : int or array_like of ints or None, optional
  19. The shape of the result. If both `s` and `axes` (see below) are None,
  20. `s` is ``x.shape``; if `s` is None but `axes` is not None, then `s` is
  21. ``numpy.take(x.shape, axes, axis=0)``.
  22. If ``s[i] > x.shape[i]``, the ith dimension of the input is padded with zeros.
  23. If ``s[i] < x.shape[i]``, the ith dimension of the input is truncated to length
  24. ``s[i]``.
  25. If any element of `s` is -1, the size of the corresponding dimension of
  26. `x` is used.
  27. axes : int or array_like of ints or None, optional
  28. Axes over which the DCT is computed. If not given, the last ``len(s)``
  29. axes are used, or all axes if `s` is also not specified.
  30. norm : {"backward", "ortho", "forward"}, optional
  31. Normalization mode (see Notes). Default is "backward".
  32. overwrite_x : bool, optional
  33. If True, the contents of `x` can be destroyed; the default is False.
  34. workers : int, optional
  35. Maximum number of workers to use for parallel computation. If negative,
  36. the value wraps around from ``os.cpu_count()``.
  37. See :func:`~scipy.fft.fft` for more details.
  38. orthogonalize : bool, optional
  39. Whether to use the orthogonalized DCT variant (see Notes).
  40. Defaults to ``True`` when ``norm="ortho"`` and ``False`` otherwise.
  41. .. versionadded:: 1.8.0
  42. Returns
  43. -------
  44. y : ndarray of real
  45. The transformed input array.
  46. See Also
  47. --------
  48. idctn : Inverse multidimensional DCT
  49. Notes
  50. -----
  51. For full details of the DCT types and normalization modes, as well as
  52. references, see `dct`.
  53. Examples
  54. --------
  55. >>> import numpy as np
  56. >>> from scipy.fft import dctn, idctn
  57. >>> rng = np.random.default_rng()
  58. >>> y = rng.standard_normal((16, 16))
  59. >>> np.allclose(y, idctn(dctn(y)))
  60. True
  61. """
  62. return (Dispatchable(x, np.ndarray),)
  63. @xp_capabilities(cpu_only=True, allow_dask_compute=True)
  64. @_dispatch
  65. def idctn(x, type=2, s=None, axes=None, norm=None, overwrite_x=False,
  66. workers=None, orthogonalize=None):
  67. """
  68. Return multidimensional Inverse Discrete Cosine Transform along the specified axes.
  69. Parameters
  70. ----------
  71. x : array_like
  72. The input array.
  73. type : {1, 2, 3, 4}, optional
  74. Type of the DCT (see Notes). Default type is 2.
  75. s : int or array_like of ints or None, optional
  76. The shape of the result. If both `s` and `axes` (see below) are
  77. None, `s` is ``x.shape``; if `s` is None but `axes` is
  78. not None, then `s` is ``numpy.take(x.shape, axes, axis=0)``.
  79. If ``s[i] > x.shape[i]``, the ith dimension of the input is padded with zeros.
  80. If ``s[i] < x.shape[i]``, the ith dimension of the input is truncated to length
  81. ``s[i]``.
  82. If any element of `s` is -1, the size of the corresponding dimension of
  83. `x` is used.
  84. axes : int or array_like of ints or None, optional
  85. Axes over which the IDCT is computed. If not given, the last ``len(s)``
  86. axes are used, or all axes if `s` is also not specified.
  87. norm : {"backward", "ortho", "forward"}, optional
  88. Normalization mode (see Notes). Default is "backward".
  89. overwrite_x : bool, optional
  90. If True, the contents of `x` can be destroyed; the default is False.
  91. workers : int, optional
  92. Maximum number of workers to use for parallel computation. If negative,
  93. the value wraps around from ``os.cpu_count()``.
  94. See :func:`~scipy.fft.fft` for more details.
  95. orthogonalize : bool, optional
  96. Whether to use the orthogonalized IDCT variant (see Notes).
  97. Defaults to ``True`` when ``norm="ortho"`` and ``False`` otherwise.
  98. .. versionadded:: 1.8.0
  99. Returns
  100. -------
  101. y : ndarray of real
  102. The transformed input array.
  103. See Also
  104. --------
  105. dctn : multidimensional DCT
  106. Notes
  107. -----
  108. For full details of the IDCT types and normalization modes, as well as
  109. references, see `idct`.
  110. Examples
  111. --------
  112. >>> import numpy as np
  113. >>> from scipy.fft import dctn, idctn
  114. >>> rng = np.random.default_rng()
  115. >>> y = rng.standard_normal((16, 16))
  116. >>> np.allclose(y, idctn(dctn(y)))
  117. True
  118. """
  119. return (Dispatchable(x, np.ndarray),)
  120. @xp_capabilities(cpu_only=True, allow_dask_compute=True)
  121. @_dispatch
  122. def dstn(x, type=2, s=None, axes=None, norm=None, overwrite_x=False,
  123. workers=None, orthogonalize=None):
  124. """
  125. Return multidimensional Discrete Sine Transform along the specified axes.
  126. Parameters
  127. ----------
  128. x : array_like
  129. The input array.
  130. type : {1, 2, 3, 4}, optional
  131. Type of the DST (see Notes). Default type is 2.
  132. s : int or array_like of ints or None, optional
  133. The shape of the result. If both `s` and `axes` (see below) are None,
  134. `s` is ``x.shape``; if `s` is None but `axes` is not None, then `s` is
  135. ``numpy.take(x.shape, axes, axis=0)``.
  136. If ``s[i] > x.shape[i]``, the ith dimension of the input is padded with zeros.
  137. If ``s[i] < x.shape[i]``, the ith dimension of the input is truncated to length
  138. ``s[i]``.
  139. If any element of `shape` is -1, the size of the corresponding dimension
  140. of `x` is used.
  141. axes : int or array_like of ints or None, optional
  142. Axes over which the DST is computed. If not given, the last ``len(s)``
  143. axes are used, or all axes if `s` is also not specified.
  144. norm : {"backward", "ortho", "forward"}, optional
  145. Normalization mode (see Notes). Default is "backward".
  146. overwrite_x : bool, optional
  147. If True, the contents of `x` can be destroyed; the default is False.
  148. workers : int, optional
  149. Maximum number of workers to use for parallel computation. If negative,
  150. the value wraps around from ``os.cpu_count()``.
  151. See :func:`~scipy.fft.fft` for more details.
  152. orthogonalize : bool, optional
  153. Whether to use the orthogonalized DST variant (see Notes).
  154. Defaults to ``True`` when ``norm="ortho"`` and ``False`` otherwise.
  155. .. versionadded:: 1.8.0
  156. Returns
  157. -------
  158. y : ndarray of real
  159. The transformed input array.
  160. See Also
  161. --------
  162. idstn : Inverse multidimensional DST
  163. Notes
  164. -----
  165. For full details of the DST types and normalization modes, as well as
  166. references, see `dst`.
  167. Examples
  168. --------
  169. >>> import numpy as np
  170. >>> from scipy.fft import dstn, idstn
  171. >>> rng = np.random.default_rng()
  172. >>> y = rng.standard_normal((16, 16))
  173. >>> np.allclose(y, idstn(dstn(y)))
  174. True
  175. """
  176. return (Dispatchable(x, np.ndarray),)
  177. @xp_capabilities(cpu_only=True, allow_dask_compute=True)
  178. @_dispatch
  179. def idstn(x, type=2, s=None, axes=None, norm=None, overwrite_x=False,
  180. workers=None, orthogonalize=None):
  181. """
  182. Return multidimensional Inverse Discrete Sine Transform along the specified axes.
  183. Parameters
  184. ----------
  185. x : array_like
  186. The input array.
  187. type : {1, 2, 3, 4}, optional
  188. Type of the DST (see Notes). Default type is 2.
  189. s : int or array_like of ints or None, optional
  190. The shape of the result. If both `s` and `axes` (see below) are None,
  191. `s` is ``x.shape``; if `s` is None but `axes` is not None, then `s` is
  192. ``numpy.take(x.shape, axes, axis=0)``.
  193. If ``s[i] > x.shape[i]``, the ith dimension of the input is padded with zeros.
  194. If ``s[i] < x.shape[i]``, the ith dimension of the input is truncated to length
  195. ``s[i]``.
  196. If any element of `s` is -1, the size of the corresponding dimension of
  197. `x` is used.
  198. axes : int or array_like of ints or None, optional
  199. Axes over which the IDST is computed. If not given, the last ``len(s)``
  200. axes are used, or all axes if `s` is also not specified.
  201. norm : {"backward", "ortho", "forward"}, optional
  202. Normalization mode (see Notes). Default is "backward".
  203. overwrite_x : bool, optional
  204. If True, the contents of `x` can be destroyed; the default is False.
  205. workers : int, optional
  206. Maximum number of workers to use for parallel computation. If negative,
  207. the value wraps around from ``os.cpu_count()``.
  208. See :func:`~scipy.fft.fft` for more details.
  209. orthogonalize : bool, optional
  210. Whether to use the orthogonalized IDST variant (see Notes).
  211. Defaults to ``True`` when ``norm="ortho"`` and ``False`` otherwise.
  212. .. versionadded:: 1.8.0
  213. Returns
  214. -------
  215. y : ndarray of real
  216. The transformed input array.
  217. See Also
  218. --------
  219. dstn : multidimensional DST
  220. Notes
  221. -----
  222. For full details of the IDST types and normalization modes, as well as
  223. references, see `idst`.
  224. Examples
  225. --------
  226. >>> import numpy as np
  227. >>> from scipy.fft import dstn, idstn
  228. >>> rng = np.random.default_rng()
  229. >>> y = rng.standard_normal((16, 16))
  230. >>> np.allclose(y, idstn(dstn(y)))
  231. True
  232. """
  233. return (Dispatchable(x, np.ndarray),)
  234. @xp_capabilities(cpu_only=True, allow_dask_compute=True,
  235. skip_backends=[('jax.numpy', 'XXX large tolerance violations')])
  236. @_dispatch
  237. def dct(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False, workers=None,
  238. orthogonalize=None):
  239. r"""Return the Discrete Cosine Transform of arbitrary type sequence x.
  240. Parameters
  241. ----------
  242. x : array_like
  243. The input array.
  244. type : {1, 2, 3, 4}, optional
  245. Type of the DCT (see Notes). Default type is 2.
  246. n : int, optional
  247. Length of the transform. If ``n < x.shape[axis]``, `x` is
  248. truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
  249. default results in ``n = x.shape[axis]``.
  250. axis : int, optional
  251. Axis along which the dct is computed; the default is over the
  252. last axis (i.e., ``axis=-1``).
  253. norm : {"backward", "ortho", "forward"}, optional
  254. Normalization mode (see Notes). Default is "backward".
  255. overwrite_x : bool, optional
  256. If True, the contents of `x` can be destroyed; the default is False.
  257. workers : int, optional
  258. Maximum number of workers to use for parallel computation. If negative,
  259. the value wraps around from ``os.cpu_count()``.
  260. See :func:`~scipy.fft.fft` for more details.
  261. orthogonalize : bool, optional
  262. Whether to use the orthogonalized DCT variant (see Notes).
  263. Defaults to ``True`` when ``norm="ortho"`` and ``False`` otherwise.
  264. .. versionadded:: 1.8.0
  265. Returns
  266. -------
  267. y : ndarray of real
  268. The transformed input array.
  269. See Also
  270. --------
  271. idct : Inverse DCT
  272. Notes
  273. -----
  274. For a single dimension array ``x``, ``dct(x, norm='ortho')`` is equal to
  275. MATLAB ``dct(x)``.
  276. .. warning:: For ``type in {1, 2, 3}``, ``norm="ortho"`` breaks the direct
  277. correspondence with the direct Fourier transform. To recover
  278. it you must specify ``orthogonalize=False``.
  279. For ``norm="ortho"`` both the `dct` and `idct` are scaled by the same
  280. overall factor in both directions. By default, the transform is also
  281. orthogonalized which for types 1, 2 and 3 means the transform definition is
  282. modified to give orthogonality of the DCT matrix (see below).
  283. For ``norm="backward"``, there is no scaling on `dct` and the `idct` is
  284. scaled by ``1/N`` where ``N`` is the "logical" size of the DCT. For
  285. ``norm="forward"`` the ``1/N`` normalization is applied to the forward
  286. `dct` instead and the `idct` is unnormalized.
  287. There are, theoretically, 8 types of the DCT, only the first 4 types are
  288. implemented in SciPy.'The' DCT generally refers to DCT type 2, and 'the'
  289. Inverse DCT generally refers to DCT type 3.
  290. **Type I**
  291. There are several definitions of the DCT-I; we use the following
  292. (for ``norm="backward"``)
  293. .. math::
  294. y_k = x_0 + (-1)^k x_{N-1} + 2 \sum_{n=1}^{N-2} x_n \cos\left(
  295. \frac{\pi k n}{N-1} \right)
  296. If ``orthogonalize=True``, ``x[0]`` and ``x[N-1]`` are multiplied by a
  297. scaling factor of :math:`\sqrt{2}`, and ``y[0]`` and ``y[N-1]`` are divided
  298. by :math:`\sqrt{2}`. When combined with ``norm="ortho"``, this makes the
  299. corresponding matrix of coefficients orthonormal (``O @ O.T = np.eye(N)``).
  300. .. note::
  301. The DCT-I is only supported for input size > 1.
  302. **Type II**
  303. There are several definitions of the DCT-II; we use the following
  304. (for ``norm="backward"``)
  305. .. math::
  306. y_k = 2 \sum_{n=0}^{N-1} x_n \cos\left(\frac{\pi k(2n+1)}{2N} \right)
  307. If ``orthogonalize=True``, ``y[0]`` is divided by :math:`\sqrt{2}` which,
  308. when combined with ``norm="ortho"``, makes the corresponding matrix of
  309. coefficients orthonormal (``O @ O.T = np.eye(N)``).
  310. **Type III**
  311. There are several definitions, we use the following (for
  312. ``norm="backward"``)
  313. .. math::
  314. y_k = x_0 + 2 \sum_{n=1}^{N-1} x_n \cos\left(\frac{\pi(2k+1)n}{2N}\right)
  315. If ``orthogonalize=True``, ``x[0]`` terms are multiplied by
  316. :math:`\sqrt{2}` which, when combined with ``norm="ortho"``, makes the
  317. corresponding matrix of coefficients orthonormal (``O @ O.T = np.eye(N)``).
  318. The (unnormalized) DCT-III is the inverse of the (unnormalized) DCT-II, up
  319. to a factor `2N`. The orthonormalized DCT-III is exactly the inverse of
  320. the orthonormalized DCT-II.
  321. **Type IV**
  322. There are several definitions of the DCT-IV; we use the following
  323. (for ``norm="backward"``)
  324. .. math::
  325. y_k = 2 \sum_{n=0}^{N-1} x_n \cos\left(\frac{\pi(2k+1)(2n+1)}{4N} \right)
  326. ``orthogonalize`` has no effect here, as the DCT-IV matrix is already
  327. orthogonal up to a scale factor of ``2N``.
  328. References
  329. ----------
  330. .. [1] 'A Fast Cosine Transform in One and Two Dimensions', by J.
  331. Makhoul, `IEEE Transactions on acoustics, speech and signal
  332. processing` vol. 28(1), pp. 27-34,
  333. :doi:`10.1109/TASSP.1980.1163351` (1980).
  334. .. [2] Wikipedia, "Discrete cosine transform",
  335. https://en.wikipedia.org/wiki/Discrete_cosine_transform
  336. Examples
  337. --------
  338. The Type 1 DCT is equivalent to the FFT (though faster) for real,
  339. even-symmetrical inputs. The output is also real and even-symmetrical.
  340. Half of the FFT input is used to generate half of the FFT output:
  341. >>> from scipy.fft import fft, dct
  342. >>> import numpy as np
  343. >>> fft(np.array([4., 3., 5., 10., 5., 3.])).real
  344. array([ 30., -8., 6., -2., 6., -8.])
  345. >>> dct(np.array([4., 3., 5., 10.]), 1)
  346. array([ 30., -8., 6., -2.])
  347. """
  348. return (Dispatchable(x, np.ndarray),)
  349. @xp_capabilities(cpu_only=True, allow_dask_compute=True)
  350. @_dispatch
  351. def idct(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False,
  352. workers=None, orthogonalize=None):
  353. """
  354. Return the Inverse Discrete Cosine Transform of an arbitrary type sequence.
  355. Parameters
  356. ----------
  357. x : array_like
  358. The input array.
  359. type : {1, 2, 3, 4}, optional
  360. Type of the DCT (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 idct is computed; the default is over the
  367. last axis (i.e., ``axis=-1``).
  368. norm : {"backward", "ortho", "forward"}, optional
  369. Normalization mode (see Notes). Default is "backward".
  370. overwrite_x : bool, optional
  371. If True, the contents of `x` can be destroyed; the default is False.
  372. workers : int, optional
  373. Maximum number of workers to use for parallel computation. If negative,
  374. the value wraps around from ``os.cpu_count()``.
  375. See :func:`~scipy.fft.fft` for more details.
  376. orthogonalize : bool, optional
  377. Whether to use the orthogonalized IDCT variant (see Notes).
  378. Defaults to ``True`` when ``norm="ortho"`` and ``False`` otherwise.
  379. .. versionadded:: 1.8.0
  380. Returns
  381. -------
  382. idct : ndarray of real
  383. The transformed input array.
  384. See Also
  385. --------
  386. dct : Forward DCT
  387. Notes
  388. -----
  389. For a single dimension array `x`, ``idct(x, norm='ortho')`` is equal to
  390. MATLAB ``idct(x)``.
  391. .. warning:: For ``type in {1, 2, 3}``, ``norm="ortho"`` breaks the direct
  392. correspondence with the inverse direct Fourier transform. To
  393. recover it you must specify ``orthogonalize=False``.
  394. For ``norm="ortho"`` both the `dct` and `idct` are scaled by the same
  395. overall factor in both directions. By default, the transform is also
  396. orthogonalized which for types 1, 2 and 3 means the transform definition is
  397. modified to give orthogonality of the IDCT matrix (see `dct` for the full
  398. definitions).
  399. 'The' IDCT is the IDCT-II, which is the same as the normalized DCT-III.
  400. The IDCT is equivalent to a normal DCT except for the normalization and
  401. type. DCT type 1 and 4 are their own inverse and DCTs 2 and 3 are each
  402. other's inverses.
  403. Examples
  404. --------
  405. The Type 1 DCT is equivalent to the DFT for real, even-symmetrical
  406. inputs. The output is also real and even-symmetrical. Half of the IFFT
  407. input is used to generate half of the IFFT output:
  408. >>> from scipy.fft import ifft, idct
  409. >>> import numpy as np
  410. >>> ifft(np.array([ 30., -8., 6., -2., 6., -8.])).real
  411. array([ 4., 3., 5., 10., 5., 3.])
  412. >>> idct(np.array([ 30., -8., 6., -2.]), 1)
  413. array([ 4., 3., 5., 10.])
  414. """
  415. return (Dispatchable(x, np.ndarray),)
  416. @xp_capabilities(cpu_only=True, allow_dask_compute=True,
  417. skip_backends=[('jax.numpy', 'XXX large tolerance violations')])
  418. @_dispatch
  419. def dst(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False, workers=None,
  420. orthogonalize=None):
  421. r"""
  422. Return the Discrete Sine Transform of arbitrary type sequence x.
  423. Parameters
  424. ----------
  425. x : array_like
  426. The input array.
  427. type : {1, 2, 3, 4}, optional
  428. Type of the DST (see Notes). Default type is 2.
  429. n : int, optional
  430. Length of the transform. If ``n < x.shape[axis]``, `x` is
  431. truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
  432. default results in ``n = x.shape[axis]``.
  433. axis : int, optional
  434. Axis along which the dst is computed; the default is over the
  435. last axis (i.e., ``axis=-1``).
  436. norm : {"backward", "ortho", "forward"}, optional
  437. Normalization mode (see Notes). Default is "backward".
  438. overwrite_x : bool, optional
  439. If True, the contents of `x` can be destroyed; the default is False.
  440. workers : int, optional
  441. Maximum number of workers to use for parallel computation. If negative,
  442. the value wraps around from ``os.cpu_count()``.
  443. See :func:`~scipy.fft.fft` for more details.
  444. orthogonalize : bool, optional
  445. Whether to use the orthogonalized DST variant (see Notes).
  446. Defaults to ``True`` when ``norm="ortho"`` and ``False`` otherwise.
  447. .. versionadded:: 1.8.0
  448. Returns
  449. -------
  450. dst : ndarray of reals
  451. The transformed input array.
  452. See Also
  453. --------
  454. idst : Inverse DST
  455. Notes
  456. -----
  457. .. warning:: For ``type in {2, 3}``, ``norm="ortho"`` breaks the direct
  458. correspondence with the direct Fourier transform. To recover
  459. it you must specify ``orthogonalize=False``.
  460. For ``norm="ortho"`` both the `dst` and `idst` are scaled by the same
  461. overall factor in both directions. By default, the transform is also
  462. orthogonalized which for types 2 and 3 means the transform definition is
  463. modified to give orthogonality of the DST matrix (see below).
  464. For ``norm="backward"``, there is no scaling on the `dst` and the `idst` is
  465. scaled by ``1/N`` where ``N`` is the "logical" size of the DST.
  466. There are, theoretically, 8 types of the DST for different combinations of
  467. even/odd boundary conditions and boundary off sets [1]_, only the first
  468. 4 types are implemented in SciPy.
  469. **Type I**
  470. There are several definitions of the DST-I; we use the following for
  471. ``norm="backward"``. DST-I assumes the input is odd around :math:`n=-1` and
  472. :math:`n=N`.
  473. .. math::
  474. y_k = 2 \sum_{n=0}^{N-1} x_n \sin\left(\frac{\pi(k+1)(n+1)}{N+1}\right)
  475. Note that the DST-I is only supported for input size > 1.
  476. The (unnormalized) DST-I is its own inverse, up to a factor :math:`2(N+1)`.
  477. The orthonormalized DST-I is exactly its own inverse.
  478. ``orthogonalize`` has no effect here, as the DST-I matrix is already
  479. orthogonal up to a scale factor of ``2N``.
  480. **Type II**
  481. There are several definitions of the DST-II; we use the following for
  482. ``norm="backward"``. DST-II assumes the input is odd around :math:`n=-1/2` and
  483. :math:`n=N-1/2`; the output is odd around :math:`k=-1` and even around :math:`k=N-1`
  484. .. math::
  485. y_k = 2 \sum_{n=0}^{N-1} x_n \sin\left(\frac{\pi(k+1)(2n+1)}{2N}\right)
  486. If ``orthogonalize=True``, ``y[-1]`` is divided :math:`\sqrt{2}` which, when
  487. combined with ``norm="ortho"``, makes the corresponding matrix of
  488. coefficients orthonormal (``O @ O.T = np.eye(N)``).
  489. **Type III**
  490. There are several definitions of the DST-III, we use the following (for
  491. ``norm="backward"``). DST-III assumes the input is odd around :math:`n=-1` and
  492. even around :math:`n=N-1`
  493. .. math::
  494. y_k = (-1)^k x_{N-1} + 2 \sum_{n=0}^{N-2} x_n \sin\left(
  495. \frac{\pi(2k+1)(n+1)}{2N}\right)
  496. If ``orthogonalize=True``, ``x[-1]`` is multiplied by :math:`\sqrt{2}`
  497. which, when combined with ``norm="ortho"``, makes the corresponding matrix
  498. of coefficients orthonormal (``O @ O.T = np.eye(N)``).
  499. The (unnormalized) DST-III is the inverse of the (unnormalized) DST-II, up
  500. to a factor :math:`2N`. The orthonormalized DST-III is exactly the inverse of the
  501. orthonormalized DST-II.
  502. **Type IV**
  503. There are several definitions of the DST-IV, we use the following (for
  504. ``norm="backward"``). DST-IV assumes the input is odd around :math:`n=-0.5` and
  505. even around :math:`n=N-0.5`
  506. .. math::
  507. y_k = 2 \sum_{n=0}^{N-1} x_n \sin\left(\frac{\pi(2k+1)(2n+1)}{4N}\right)
  508. ``orthogonalize`` has no effect here, as the DST-IV matrix is already
  509. orthogonal up to a scale factor of ``2N``.
  510. The (unnormalized) DST-IV is its own inverse, up to a factor :math:`2N`. The
  511. orthonormalized DST-IV is exactly its own inverse.
  512. Examples
  513. --------
  514. Compute the DST of a simple 1D array:
  515. >>> import numpy as np
  516. >>> from scipy.fft import dst
  517. >>> x = np.array([1, -1, 1, -1])
  518. >>> dst(x, type=2)
  519. array([0., 0., 0., 8.])
  520. This computes the Discrete Sine Transform (DST) of type-II for the input array.
  521. The output contains the transformed values corresponding to the given input sequence
  522. References
  523. ----------
  524. .. [1] Wikipedia, "Discrete sine transform",
  525. https://en.wikipedia.org/wiki/Discrete_sine_transform
  526. """
  527. return (Dispatchable(x, np.ndarray),)
  528. @xp_capabilities(cpu_only=True, allow_dask_compute=True)
  529. @_dispatch
  530. def idst(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False,
  531. workers=None, orthogonalize=None):
  532. """
  533. Return the Inverse Discrete Sine Transform of an arbitrary type sequence.
  534. Parameters
  535. ----------
  536. x : array_like
  537. The input array.
  538. type : {1, 2, 3, 4}, optional
  539. Type of the DST (see Notes). Default type is 2.
  540. n : int, optional
  541. Length of the transform. If ``n < x.shape[axis]``, `x` is
  542. truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
  543. default results in ``n = x.shape[axis]``.
  544. axis : int, optional
  545. Axis along which the idst is computed; the default is over the
  546. last axis (i.e., ``axis=-1``).
  547. norm : {"backward", "ortho", "forward"}, optional
  548. Normalization mode (see Notes). Default is "backward".
  549. overwrite_x : bool, optional
  550. If True, the contents of `x` can be destroyed; the default is False.
  551. workers : int, optional
  552. Maximum number of workers to use for parallel computation. If negative,
  553. the value wraps around from ``os.cpu_count()``.
  554. See :func:`~scipy.fft.fft` for more details.
  555. orthogonalize : bool, optional
  556. Whether to use the orthogonalized IDST variant (see Notes).
  557. Defaults to ``True`` when ``norm="ortho"`` and ``False`` otherwise.
  558. .. versionadded:: 1.8.0
  559. Returns
  560. -------
  561. idst : ndarray of real
  562. The transformed input array.
  563. See Also
  564. --------
  565. dst : Forward DST
  566. Notes
  567. -----
  568. .. warning:: For ``type in {2, 3}``, ``norm="ortho"`` breaks the direct
  569. correspondence with the inverse direct Fourier transform.
  570. For ``norm="ortho"`` both the `dst` and `idst` are scaled by the same
  571. overall factor in both directions. By default, the transform is also
  572. orthogonalized which for types 2 and 3 means the transform definition is
  573. modified to give orthogonality of the DST matrix (see `dst` for the full
  574. definitions).
  575. 'The' IDST is the IDST-II, which is the same as the normalized DST-III.
  576. The IDST is equivalent to a normal DST except for the normalization and
  577. type. DST type 1 and 4 are their own inverse and DSTs 2 and 3 are each
  578. other's inverses.
  579. """
  580. return (Dispatchable(x, np.ndarray),)