Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/misc/common.py : 10%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2Functions which are common and require SciPy Base and Level 1 SciPy
3(special, linalg)
4"""
6from numpy import arange, newaxis, hstack, prod, array, frombuffer, load
8__all__ = ['central_diff_weights', 'derivative', 'ascent', 'face',
9 'electrocardiogram']
12def central_diff_weights(Np, ndiv=1):
13 """
14 Return weights for an Np-point central derivative.
16 Assumes equally-spaced function points.
18 If weights are in the vector w, then
19 derivative is w[0] * f(x-ho*dx) + ... + w[-1] * f(x+h0*dx)
21 Parameters
22 ----------
23 Np : int
24 Number of points for the central derivative.
25 ndiv : int, optional
26 Number of divisions. Default is 1.
28 Returns
29 -------
30 w : ndarray
31 Weights for an Np-point central derivative. Its size is `Np`.
33 Notes
34 -----
35 Can be inaccurate for a large number of points.
37 Examples
38 --------
39 We can calculate a derivative value of a function.
41 >>> from scipy.misc import central_diff_weights
42 >>> def f(x):
43 ... return 2 * x**2 + 3
44 >>> x = 3.0 # derivative point
45 >>> h = 0.1 # differential step
46 >>> Np = 3 # point number for central derivative
47 >>> weights = central_diff_weights(Np) # weights for first derivative
48 >>> vals = [f(x + (i - Np/2) * h) for i in range(Np)]
49 >>> sum(w * v for (w, v) in zip(weights, vals))/h
50 11.79999999999998
52 This value is close to the analytical solution:
53 f'(x) = 4x, so f'(3) = 12
55 References
56 ----------
57 .. [1] https://en.wikipedia.org/wiki/Finite_difference
59 """
60 if Np < ndiv + 1:
61 raise ValueError("Number of points must be at least the derivative order + 1.")
62 if Np % 2 == 0:
63 raise ValueError("The number of points must be odd.")
64 from scipy import linalg
65 ho = Np >> 1
66 x = arange(-ho,ho+1.0)
67 x = x[:,newaxis]
68 X = x**0.0
69 for k in range(1,Np):
70 X = hstack([X,x**k])
71 w = prod(arange(1,ndiv+1),axis=0)*linalg.inv(X)[ndiv]
72 return w
75def derivative(func, x0, dx=1.0, n=1, args=(), order=3):
76 """
77 Find the nth derivative of a function at a point.
79 Given a function, use a central difference formula with spacing `dx` to
80 compute the nth derivative at `x0`.
82 Parameters
83 ----------
84 func : function
85 Input function.
86 x0 : float
87 The point at which the nth derivative is found.
88 dx : float, optional
89 Spacing.
90 n : int, optional
91 Order of the derivative. Default is 1.
92 args : tuple, optional
93 Arguments
94 order : int, optional
95 Number of points to use, must be odd.
97 Notes
98 -----
99 Decreasing the step size too small can result in round-off error.
101 Examples
102 --------
103 >>> from scipy.misc import derivative
104 >>> def f(x):
105 ... return x**3 + x**2
106 >>> derivative(f, 1.0, dx=1e-6)
107 4.9999999999217337
109 """
110 if order < n + 1:
111 raise ValueError("'order' (the number of points used to compute the derivative), "
112 "must be at least the derivative order 'n' + 1.")
113 if order % 2 == 0:
114 raise ValueError("'order' (the number of points used to compute the derivative) "
115 "must be odd.")
116 # pre-computed for n=1 and 2 and low-order for speed.
117 if n == 1:
118 if order == 3:
119 weights = array([-1,0,1])/2.0
120 elif order == 5:
121 weights = array([1,-8,0,8,-1])/12.0
122 elif order == 7:
123 weights = array([-1,9,-45,0,45,-9,1])/60.0
124 elif order == 9:
125 weights = array([3,-32,168,-672,0,672,-168,32,-3])/840.0
126 else:
127 weights = central_diff_weights(order,1)
128 elif n == 2:
129 if order == 3:
130 weights = array([1,-2.0,1])
131 elif order == 5:
132 weights = array([-1,16,-30,16,-1])/12.0
133 elif order == 7:
134 weights = array([2,-27,270,-490,270,-27,2])/180.0
135 elif order == 9:
136 weights = array([-9,128,-1008,8064,-14350,8064,-1008,128,-9])/5040.0
137 else:
138 weights = central_diff_weights(order,2)
139 else:
140 weights = central_diff_weights(order, n)
141 val = 0.0
142 ho = order >> 1
143 for k in range(order):
144 val += weights[k]*func(x0+(k-ho)*dx,*args)
145 return val / prod((dx,)*n,axis=0)
148def ascent():
149 """
150 Get an 8-bit grayscale bit-depth, 512 x 512 derived image for easy use in demos
152 The image is derived from accent-to-the-top.jpg at
153 http://www.public-domain-image.com/people-public-domain-images-pictures/
155 Parameters
156 ----------
157 None
159 Returns
160 -------
161 ascent : ndarray
162 convenient image to use for testing and demonstration
164 Examples
165 --------
166 >>> import scipy.misc
167 >>> ascent = scipy.misc.ascent()
168 >>> ascent.shape
169 (512, 512)
170 >>> ascent.max()
171 255
173 >>> import matplotlib.pyplot as plt
174 >>> plt.gray()
175 >>> plt.imshow(ascent)
176 >>> plt.show()
178 """
179 import pickle
180 import os
181 fname = os.path.join(os.path.dirname(__file__),'ascent.dat')
182 with open(fname, 'rb') as f:
183 ascent = array(pickle.load(f))
184 return ascent
187def face(gray=False):
188 """
189 Get a 1024 x 768, color image of a raccoon face.
191 raccoon-procyon-lotor.jpg at http://www.public-domain-image.com
193 Parameters
194 ----------
195 gray : bool, optional
196 If True return 8-bit grey-scale image, otherwise return a color image
198 Returns
199 -------
200 face : ndarray
201 image of a racoon face
203 Examples
204 --------
205 >>> import scipy.misc
206 >>> face = scipy.misc.face()
207 >>> face.shape
208 (768, 1024, 3)
209 >>> face.max()
210 255
211 >>> face.dtype
212 dtype('uint8')
214 >>> import matplotlib.pyplot as plt
215 >>> plt.gray()
216 >>> plt.imshow(face)
217 >>> plt.show()
219 """
220 import bz2
221 import os
222 with open(os.path.join(os.path.dirname(__file__), 'face.dat'), 'rb') as f:
223 rawdata = f.read()
224 data = bz2.decompress(rawdata)
225 face = frombuffer(data, dtype='uint8')
226 face.shape = (768, 1024, 3)
227 if gray is True:
228 face = (0.21 * face[:,:,0] + 0.71 * face[:,:,1] + 0.07 * face[:,:,2]).astype('uint8')
229 return face
232def electrocardiogram():
233 """
234 Load an electrocardiogram as an example for a 1-D signal.
236 The returned signal is a 5 minute long electrocardiogram (ECG), a medical
237 recording of the heart's electrical activity, sampled at 360 Hz.
239 Returns
240 -------
241 ecg : ndarray
242 The electrocardiogram in millivolt (mV) sampled at 360 Hz.
244 Notes
245 -----
246 The provided signal is an excerpt (19:35 to 24:35) from the `record 208`_
247 (lead MLII) provided by the MIT-BIH Arrhythmia Database [1]_ on
248 PhysioNet [2]_. The excerpt includes noise induced artifacts, typical
249 heartbeats as well as pathological changes.
251 .. _record 208: https://physionet.org/physiobank/database/html/mitdbdir/records.htm#208
253 .. versionadded:: 1.1.0
255 References
256 ----------
257 .. [1] Moody GB, Mark RG. The impact of the MIT-BIH Arrhythmia Database.
258 IEEE Eng in Med and Biol 20(3):45-50 (May-June 2001).
259 (PMID: 11446209); :doi:`10.13026/C2F305`
260 .. [2] Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh,
261 Mark RG, Mietus JE, Moody GB, Peng C-K, Stanley HE. PhysioBank,
262 PhysioToolkit, and PhysioNet: Components of a New Research Resource
263 for Complex Physiologic Signals. Circulation 101(23):e215-e220;
264 :doi:`10.1161/01.CIR.101.23.e215`
266 Examples
267 --------
268 >>> from scipy.misc import electrocardiogram
269 >>> ecg = electrocardiogram()
270 >>> ecg
271 array([-0.245, -0.215, -0.185, ..., -0.405, -0.395, -0.385])
272 >>> ecg.shape, ecg.mean(), ecg.std()
273 ((108000,), -0.16510875, 0.5992473991177294)
275 As stated the signal features several areas with a different morphology.
276 E.g., the first few seconds show the electrical activity of a heart in
277 normal sinus rhythm as seen below.
279 >>> import matplotlib.pyplot as plt
280 >>> fs = 360
281 >>> time = np.arange(ecg.size) / fs
282 >>> plt.plot(time, ecg)
283 >>> plt.xlabel("time in s")
284 >>> plt.ylabel("ECG in mV")
285 >>> plt.xlim(9, 10.2)
286 >>> plt.ylim(-1, 1.5)
287 >>> plt.show()
289 After second 16, however, the first premature ventricular contractions, also
290 called extrasystoles, appear. These have a different morphology compared to
291 typical heartbeats. The difference can easily be observed in the following
292 plot.
294 >>> plt.plot(time, ecg)
295 >>> plt.xlabel("time in s")
296 >>> plt.ylabel("ECG in mV")
297 >>> plt.xlim(46.5, 50)
298 >>> plt.ylim(-2, 1.5)
299 >>> plt.show()
301 At several points large artifacts disturb the recording, e.g.:
303 >>> plt.plot(time, ecg)
304 >>> plt.xlabel("time in s")
305 >>> plt.ylabel("ECG in mV")
306 >>> plt.xlim(207, 215)
307 >>> plt.ylim(-2, 3.5)
308 >>> plt.show()
310 Finally, examining the power spectrum reveals that most of the biosignal is
311 made up of lower frequencies. At 60 Hz the noise induced by the mains
312 electricity can be clearly observed.
314 >>> from scipy.signal import welch
315 >>> f, Pxx = welch(ecg, fs=fs, nperseg=2048, scaling="spectrum")
316 >>> plt.semilogy(f, Pxx)
317 >>> plt.xlabel("Frequency in Hz")
318 >>> plt.ylabel("Power spectrum of the ECG in mV**2")
319 >>> plt.xlim(f[[0, -1]])
320 >>> plt.show()
321 """
322 import os
323 file_path = os.path.join(os.path.dirname(__file__), "ecg.dat")
324 with load(file_path) as file:
325 ecg = file["ecg"].astype(int) # np.uint16 -> int
326 # Convert raw output of ADC to mV: (ecg - adc_zero) / adc_gain
327 ecg = (ecg - 1024) / 200.0
328 return ecg