_common.py 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342
  1. """
  2. Functions which are common and require SciPy Base and Level 1 SciPy
  3. (special, linalg)
  4. """
  5. from scipy._lib.deprecation import _deprecated
  6. from scipy._lib._finite_differences import _central_diff_weights, _derivative
  7. from numpy import array, frombuffer, load
  8. __all__ = ['central_diff_weights', 'derivative', 'ascent', 'face',
  9. 'electrocardiogram']
  10. @_deprecated(msg="scipy.misc.central_diff_weights is deprecated in "
  11. "SciPy v1.10.0; and will be completely removed in "
  12. "SciPy v1.12.0. You may consider using "
  13. "findiff: https://github.com/maroba/findiff or "
  14. "numdifftools: https://github.com/pbrod/numdifftools")
  15. def central_diff_weights(Np, ndiv=1):
  16. """
  17. Return weights for an Np-point central derivative.
  18. Assumes equally-spaced function points.
  19. If weights are in the vector w, then
  20. derivative is w[0] * f(x-ho*dx) + ... + w[-1] * f(x+h0*dx)
  21. .. deprecated:: 1.10.0
  22. `central_diff_weights` has been deprecated from
  23. `scipy.misc.central_diff_weights` in SciPy 1.10.0 and
  24. it will be completely removed in SciPy 1.12.0.
  25. You may consider using
  26. findiff: https://github.com/maroba/findiff or
  27. numdifftools: https://github.com/pbrod/numdifftools
  28. Parameters
  29. ----------
  30. Np : int
  31. Number of points for the central derivative.
  32. ndiv : int, optional
  33. Number of divisions. Default is 1.
  34. Returns
  35. -------
  36. w : ndarray
  37. Weights for an Np-point central derivative. Its size is `Np`.
  38. Notes
  39. -----
  40. Can be inaccurate for a large number of points.
  41. Examples
  42. --------
  43. We can calculate a derivative value of a function.
  44. >>> from scipy.misc import central_diff_weights
  45. >>> def f(x):
  46. ... return 2 * x**2 + 3
  47. >>> x = 3.0 # derivative point
  48. >>> h = 0.1 # differential step
  49. >>> Np = 3 # point number for central derivative
  50. >>> weights = central_diff_weights(Np) # weights for first derivative
  51. >>> vals = [f(x + (i - Np/2) * h) for i in range(Np)]
  52. >>> sum(w * v for (w, v) in zip(weights, vals))/h
  53. 11.79999999999998
  54. This value is close to the analytical solution:
  55. f'(x) = 4x, so f'(3) = 12
  56. References
  57. ----------
  58. .. [1] https://en.wikipedia.org/wiki/Finite_difference
  59. """
  60. return _central_diff_weights(Np, ndiv)
  61. @_deprecated(msg="scipy.misc.derivative is deprecated in "
  62. "SciPy v1.10.0; and will be completely removed in "
  63. "SciPy v1.12.0. You may consider using "
  64. "findiff: https://github.com/maroba/findiff or "
  65. "numdifftools: https://github.com/pbrod/numdifftools")
  66. def derivative(func, x0, dx=1.0, n=1, args=(), order=3):
  67. """
  68. Find the nth derivative of a function at a point.
  69. Given a function, use a central difference formula with spacing `dx` to
  70. compute the nth derivative at `x0`.
  71. .. deprecated:: 1.10.0
  72. `derivative` has been deprecated from `scipy.misc.derivative`
  73. in SciPy 1.10.0 and it will be completely removed in SciPy 1.12.0.
  74. You may consider using
  75. findiff: https://github.com/maroba/findiff or
  76. numdifftools: https://github.com/pbrod/numdifftools
  77. Parameters
  78. ----------
  79. func : function
  80. Input function.
  81. x0 : float
  82. The point at which the nth derivative is found.
  83. dx : float, optional
  84. Spacing.
  85. n : int, optional
  86. Order of the derivative. Default is 1.
  87. args : tuple, optional
  88. Arguments
  89. order : int, optional
  90. Number of points to use, must be odd.
  91. Notes
  92. -----
  93. Decreasing the step size too small can result in round-off error.
  94. Examples
  95. --------
  96. >>> from scipy.misc import derivative
  97. >>> def f(x):
  98. ... return x**3 + x**2
  99. >>> derivative(f, 1.0, dx=1e-6)
  100. 4.9999999999217337
  101. """
  102. return _derivative(func, x0, dx, n, args, order)
  103. @_deprecated(msg="scipy.misc.ascent has been deprecated in SciPy v1.10.0;"
  104. " and will be completely removed in SciPy v1.12.0. "
  105. "Dataset methods have moved into the scipy.datasets "
  106. "module. Use scipy.datasets.ascent instead.")
  107. def ascent():
  108. """
  109. Get an 8-bit grayscale bit-depth, 512 x 512 derived image for easy use in demos
  110. The image is derived from accent-to-the-top.jpg at
  111. http://www.public-domain-image.com/people-public-domain-images-pictures/
  112. .. deprecated:: 1.10.0
  113. `ascent` has been deprecated from `scipy.misc.ascent`
  114. in SciPy 1.10.0 and it will be completely removed in SciPy 1.12.0.
  115. Dataset methods have moved into the `scipy.datasets` module.
  116. Use `scipy.datasets.ascent` instead.
  117. Parameters
  118. ----------
  119. None
  120. Returns
  121. -------
  122. ascent : ndarray
  123. convenient image to use for testing and demonstration
  124. Examples
  125. --------
  126. >>> import scipy.misc
  127. >>> ascent = scipy.misc.ascent()
  128. >>> ascent.shape
  129. (512, 512)
  130. >>> ascent.max()
  131. 255
  132. >>> import matplotlib.pyplot as plt
  133. >>> plt.gray()
  134. >>> plt.imshow(ascent)
  135. >>> plt.show()
  136. """
  137. import pickle
  138. import os
  139. fname = os.path.join(os.path.dirname(__file__),'ascent.dat')
  140. with open(fname, 'rb') as f:
  141. ascent = array(pickle.load(f))
  142. return ascent
  143. @_deprecated(msg="scipy.misc.face has been deprecated in SciPy v1.10.0; "
  144. "and will be completely removed in SciPy v1.12.0. "
  145. "Dataset methods have moved into the scipy.datasets "
  146. "module. Use scipy.datasets.face instead.")
  147. def face(gray=False):
  148. """
  149. Get a 1024 x 768, color image of a raccoon face.
  150. raccoon-procyon-lotor.jpg at http://www.public-domain-image.com
  151. .. deprecated:: 1.10.0
  152. `face` has been deprecated from `scipy.misc.face`
  153. in SciPy 1.10.0 and it will be completely removed in SciPy 1.12.0.
  154. Dataset methods have moved into the `scipy.datasets` module.
  155. Use `scipy.datasets.face` instead.
  156. Parameters
  157. ----------
  158. gray : bool, optional
  159. If True return 8-bit grey-scale image, otherwise return a color image
  160. Returns
  161. -------
  162. face : ndarray
  163. image of a racoon face
  164. Examples
  165. --------
  166. >>> import scipy.misc
  167. >>> face = scipy.misc.face()
  168. >>> face.shape
  169. (768, 1024, 3)
  170. >>> face.max()
  171. 255
  172. >>> face.dtype
  173. dtype('uint8')
  174. >>> import matplotlib.pyplot as plt
  175. >>> plt.gray()
  176. >>> plt.imshow(face)
  177. >>> plt.show()
  178. """
  179. import bz2
  180. import os
  181. with open(os.path.join(os.path.dirname(__file__), 'face.dat'), 'rb') as f:
  182. rawdata = f.read()
  183. data = bz2.decompress(rawdata)
  184. face = frombuffer(data, dtype='uint8')
  185. face.shape = (768, 1024, 3)
  186. if gray is True:
  187. face = (0.21 * face[:,:,0] + 0.71 * face[:,:,1] + 0.07 * face[:,:,2]).astype('uint8')
  188. return face
  189. @_deprecated(msg="scipy.misc.electrocardiogram has been "
  190. "deprecated in SciPy v1.10.0; and will "
  191. "be completely removed in SciPy v1.12.0. "
  192. "Dataset methods have moved into the scipy.datasets "
  193. "module. Use scipy.datasets.electrocardiogram instead.")
  194. def electrocardiogram():
  195. """
  196. Load an electrocardiogram as an example for a 1-D signal.
  197. The returned signal is a 5 minute long electrocardiogram (ECG), a medical
  198. recording of the heart's electrical activity, sampled at 360 Hz.
  199. .. deprecated:: 1.10.0
  200. `electrocardiogram` has been deprecated from
  201. `scipy.misc.electrocardiogram` in SciPy 1.10.0 and it will be
  202. completely removed in SciPy 1.12.0.
  203. Dataset methods have moved into the `scipy.datasets` module.
  204. Use `scipy.datasets.electrocardiogram` instead.
  205. Returns
  206. -------
  207. ecg : ndarray
  208. The electrocardiogram in millivolt (mV) sampled at 360 Hz.
  209. Notes
  210. -----
  211. The provided signal is an excerpt (19:35 to 24:35) from the `record 208`_
  212. (lead MLII) provided by the MIT-BIH Arrhythmia Database [1]_ on
  213. PhysioNet [2]_. The excerpt includes noise induced artifacts, typical
  214. heartbeats as well as pathological changes.
  215. .. _record 208: https://physionet.org/physiobank/database/html/mitdbdir/records.htm#208
  216. .. versionadded:: 1.1.0
  217. References
  218. ----------
  219. .. [1] Moody GB, Mark RG. The impact of the MIT-BIH Arrhythmia Database.
  220. IEEE Eng in Med and Biol 20(3):45-50 (May-June 2001).
  221. (PMID: 11446209); :doi:`10.13026/C2F305`
  222. .. [2] Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh,
  223. Mark RG, Mietus JE, Moody GB, Peng C-K, Stanley HE. PhysioBank,
  224. PhysioToolkit, and PhysioNet: Components of a New Research Resource
  225. for Complex Physiologic Signals. Circulation 101(23):e215-e220;
  226. :doi:`10.1161/01.CIR.101.23.e215`
  227. Examples
  228. --------
  229. >>> from scipy.misc import electrocardiogram
  230. >>> ecg = electrocardiogram()
  231. >>> ecg
  232. array([-0.245, -0.215, -0.185, ..., -0.405, -0.395, -0.385])
  233. >>> ecg.shape, ecg.mean(), ecg.std()
  234. ((108000,), -0.16510875, 0.5992473991177294)
  235. As stated the signal features several areas with a different morphology.
  236. E.g., the first few seconds show the electrical activity of a heart in
  237. normal sinus rhythm as seen below.
  238. >>> import numpy as np
  239. >>> import matplotlib.pyplot as plt
  240. >>> fs = 360
  241. >>> time = np.arange(ecg.size) / fs
  242. >>> plt.plot(time, ecg)
  243. >>> plt.xlabel("time in s")
  244. >>> plt.ylabel("ECG in mV")
  245. >>> plt.xlim(9, 10.2)
  246. >>> plt.ylim(-1, 1.5)
  247. >>> plt.show()
  248. After second 16, however, the first premature ventricular contractions, also
  249. called extrasystoles, appear. These have a different morphology compared to
  250. typical heartbeats. The difference can easily be observed in the following
  251. plot.
  252. >>> plt.plot(time, ecg)
  253. >>> plt.xlabel("time in s")
  254. >>> plt.ylabel("ECG in mV")
  255. >>> plt.xlim(46.5, 50)
  256. >>> plt.ylim(-2, 1.5)
  257. >>> plt.show()
  258. At several points large artifacts disturb the recording, e.g.:
  259. >>> plt.plot(time, ecg)
  260. >>> plt.xlabel("time in s")
  261. >>> plt.ylabel("ECG in mV")
  262. >>> plt.xlim(207, 215)
  263. >>> plt.ylim(-2, 3.5)
  264. >>> plt.show()
  265. Finally, examining the power spectrum reveals that most of the biosignal is
  266. made up of lower frequencies. At 60 Hz the noise induced by the mains
  267. electricity can be clearly observed.
  268. >>> from scipy.signal import welch
  269. >>> f, Pxx = welch(ecg, fs=fs, nperseg=2048, scaling="spectrum")
  270. >>> plt.semilogy(f, Pxx)
  271. >>> plt.xlabel("Frequency in Hz")
  272. >>> plt.ylabel("Power spectrum of the ECG in mV**2")
  273. >>> plt.xlim(f[[0, -1]])
  274. >>> plt.show()
  275. """
  276. import os
  277. file_path = os.path.join(os.path.dirname(__file__), "ecg.dat")
  278. with load(file_path) as file:
  279. ecg = file["ecg"].astype(int) # np.uint16 -> int
  280. # Convert raw output of ADC to mV: (ecg - adc_zero) / adc_gain
  281. ecg = (ecg - 1024) / 200.0
  282. return ecg