Hide keyboard shortcuts

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""" 

5 

6from numpy import arange, newaxis, hstack, prod, array, frombuffer, load 

7 

8__all__ = ['central_diff_weights', 'derivative', 'ascent', 'face', 

9 'electrocardiogram'] 

10 

11 

12def central_diff_weights(Np, ndiv=1): 

13 """ 

14 Return weights for an Np-point central derivative. 

15 

16 Assumes equally-spaced function points. 

17 

18 If weights are in the vector w, then 

19 derivative is w[0] * f(x-ho*dx) + ... + w[-1] * f(x+h0*dx) 

20 

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. 

27 

28 Returns 

29 ------- 

30 w : ndarray 

31 Weights for an Np-point central derivative. Its size is `Np`. 

32 

33 Notes 

34 ----- 

35 Can be inaccurate for a large number of points. 

36 

37 Examples 

38 -------- 

39 We can calculate a derivative value of a function. 

40 

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 

51 

52 This value is close to the analytical solution: 

53 f'(x) = 4x, so f'(3) = 12 

54 

55 References 

56 ---------- 

57 .. [1] https://en.wikipedia.org/wiki/Finite_difference 

58 

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 

73 

74 

75def derivative(func, x0, dx=1.0, n=1, args=(), order=3): 

76 """ 

77 Find the nth derivative of a function at a point. 

78 

79 Given a function, use a central difference formula with spacing `dx` to 

80 compute the nth derivative at `x0`. 

81 

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. 

96 

97 Notes 

98 ----- 

99 Decreasing the step size too small can result in round-off error. 

100 

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 

108 

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) 

146 

147 

148def ascent(): 

149 """ 

150 Get an 8-bit grayscale bit-depth, 512 x 512 derived image for easy use in demos 

151 

152 The image is derived from accent-to-the-top.jpg at 

153 http://www.public-domain-image.com/people-public-domain-images-pictures/ 

154 

155 Parameters 

156 ---------- 

157 None 

158 

159 Returns 

160 ------- 

161 ascent : ndarray 

162 convenient image to use for testing and demonstration 

163 

164 Examples 

165 -------- 

166 >>> import scipy.misc 

167 >>> ascent = scipy.misc.ascent() 

168 >>> ascent.shape 

169 (512, 512) 

170 >>> ascent.max() 

171 255 

172 

173 >>> import matplotlib.pyplot as plt 

174 >>> plt.gray() 

175 >>> plt.imshow(ascent) 

176 >>> plt.show() 

177 

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 

185 

186 

187def face(gray=False): 

188 """ 

189 Get a 1024 x 768, color image of a raccoon face. 

190 

191 raccoon-procyon-lotor.jpg at http://www.public-domain-image.com 

192 

193 Parameters 

194 ---------- 

195 gray : bool, optional 

196 If True return 8-bit grey-scale image, otherwise return a color image 

197 

198 Returns 

199 ------- 

200 face : ndarray 

201 image of a racoon face 

202 

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') 

213 

214 >>> import matplotlib.pyplot as plt 

215 >>> plt.gray() 

216 >>> plt.imshow(face) 

217 >>> plt.show() 

218 

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 

230 

231 

232def electrocardiogram(): 

233 """ 

234 Load an electrocardiogram as an example for a 1-D signal. 

235 

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. 

238 

239 Returns 

240 ------- 

241 ecg : ndarray 

242 The electrocardiogram in millivolt (mV) sampled at 360 Hz. 

243 

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. 

250 

251 .. _record 208: https://physionet.org/physiobank/database/html/mitdbdir/records.htm#208 

252 

253 .. versionadded:: 1.1.0 

254 

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` 

265 

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) 

274 

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. 

278 

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() 

288 

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. 

293 

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() 

300 

301 At several points large artifacts disturb the recording, e.g.: 

302 

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() 

309 

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. 

313 

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