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

1from numpy import (logical_and, asarray, pi, zeros_like, 

2 piecewise, array, arctan2, tan, zeros, arange, floor) 

3from numpy.core.umath import (sqrt, exp, greater, less, cos, add, sin, 

4 less_equal, greater_equal) 

5 

6# From splinemodule.c 

7from .spline import cspline2d, sepfir2d 

8 

9from scipy.special import comb, gamma 

10 

11__all__ = ['spline_filter', 'bspline', 'gauss_spline', 'cubic', 'quadratic', 

12 'cspline1d', 'qspline1d', 'cspline1d_eval', 'qspline1d_eval'] 

13 

14 

15def factorial(n): 

16 return gamma(n + 1) 

17 

18 

19def spline_filter(Iin, lmbda=5.0): 

20 """Smoothing spline (cubic) filtering of a rank-2 array. 

21 

22 Filter an input data set, `Iin`, using a (cubic) smoothing spline of 

23 fall-off `lmbda`. 

24 """ 

25 intype = Iin.dtype.char 

26 hcol = array([1.0, 4.0, 1.0], 'f') / 6.0 

27 if intype in ['F', 'D']: 

28 Iin = Iin.astype('F') 

29 ckr = cspline2d(Iin.real, lmbda) 

30 cki = cspline2d(Iin.imag, lmbda) 

31 outr = sepfir2d(ckr, hcol, hcol) 

32 outi = sepfir2d(cki, hcol, hcol) 

33 out = (outr + 1j * outi).astype(intype) 

34 elif intype in ['f', 'd']: 

35 ckr = cspline2d(Iin, lmbda) 

36 out = sepfir2d(ckr, hcol, hcol) 

37 out = out.astype(intype) 

38 else: 

39 raise TypeError("Invalid data type for Iin") 

40 return out 

41 

42 

43_splinefunc_cache = {} 

44 

45 

46def _bspline_piecefunctions(order): 

47 """Returns the function defined over the left-side pieces for a bspline of 

48 a given order. 

49 

50 The 0th piece is the first one less than 0. The last piece is a function 

51 identical to 0 (returned as the constant 0). (There are order//2 + 2 total 

52 pieces). 

53 

54 Also returns the condition functions that when evaluated return boolean 

55 arrays for use with `numpy.piecewise`. 

56 """ 

57 try: 

58 return _splinefunc_cache[order] 

59 except KeyError: 

60 pass 

61 

62 def condfuncgen(num, val1, val2): 

63 if num == 0: 

64 return lambda x: logical_and(less_equal(x, val1), 

65 greater_equal(x, val2)) 

66 elif num == 2: 

67 return lambda x: less_equal(x, val2) 

68 else: 

69 return lambda x: logical_and(less(x, val1), 

70 greater_equal(x, val2)) 

71 

72 last = order // 2 + 2 

73 if order % 2: 

74 startbound = -1.0 

75 else: 

76 startbound = -0.5 

77 condfuncs = [condfuncgen(0, 0, startbound)] 

78 bound = startbound 

79 for num in range(1, last - 1): 

80 condfuncs.append(condfuncgen(1, bound, bound - 1)) 

81 bound = bound - 1 

82 condfuncs.append(condfuncgen(2, 0, -(order + 1) / 2.0)) 

83 

84 # final value of bound is used in piecefuncgen below 

85 

86 # the functions to evaluate are taken from the left-hand side 

87 # in the general expression derived from the central difference 

88 # operator (because they involve fewer terms). 

89 

90 fval = factorial(order) 

91 

92 def piecefuncgen(num): 

93 Mk = order // 2 - num 

94 if (Mk < 0): 

95 return 0 # final function is 0 

96 coeffs = [(1 - 2 * (k % 2)) * float(comb(order + 1, k, exact=1)) / fval 

97 for k in range(Mk + 1)] 

98 shifts = [-bound - k for k in range(Mk + 1)] 

99 

100 def thefunc(x): 

101 res = 0.0 

102 for k in range(Mk + 1): 

103 res += coeffs[k] * (x + shifts[k]) ** order 

104 return res 

105 return thefunc 

106 

107 funclist = [piecefuncgen(k) for k in range(last)] 

108 

109 _splinefunc_cache[order] = (funclist, condfuncs) 

110 

111 return funclist, condfuncs 

112 

113 

114def bspline(x, n): 

115 """B-spline basis function of order n. 

116 

117 Notes 

118 ----- 

119 Uses numpy.piecewise and automatic function-generator. 

120 

121 """ 

122 ax = -abs(asarray(x)) 

123 # number of pieces on the left-side is (n+1)/2 

124 funclist, condfuncs = _bspline_piecefunctions(n) 

125 condlist = [func(ax) for func in condfuncs] 

126 return piecewise(ax, condlist, funclist) 

127 

128 

129def gauss_spline(x, n): 

130 """Gaussian approximation to B-spline basis function of order n. 

131 

132 Parameters 

133 ---------- 

134 n : int 

135 The order of the spline. Must be nonnegative, i.e., n >= 0 

136 

137 References 

138 ---------- 

139 .. [1] Bouma H., Vilanova A., Bescos J.O., ter Haar Romeny B.M., Gerritsen 

140 F.A. (2007) Fast and Accurate Gaussian Derivatives Based on B-Splines. In: 

141 Sgallari F., Murli A., Paragios N. (eds) Scale Space and Variational 

142 Methods in Computer Vision. SSVM 2007. Lecture Notes in Computer 

143 Science, vol 4485. Springer, Berlin, Heidelberg 

144 """ 

145 signsq = (n + 1) / 12.0 

146 return 1 / sqrt(2 * pi * signsq) * exp(-x ** 2 / 2 / signsq) 

147 

148 

149def cubic(x): 

150 """A cubic B-spline. 

151 

152 This is a special case of `bspline`, and equivalent to ``bspline(x, 3)``. 

153 """ 

154 ax = abs(asarray(x)) 

155 res = zeros_like(ax) 

156 cond1 = less(ax, 1) 

157 if cond1.any(): 

158 ax1 = ax[cond1] 

159 res[cond1] = 2.0 / 3 - 1.0 / 2 * ax1 ** 2 * (2 - ax1) 

160 cond2 = ~cond1 & less(ax, 2) 

161 if cond2.any(): 

162 ax2 = ax[cond2] 

163 res[cond2] = 1.0 / 6 * (2 - ax2) ** 3 

164 return res 

165 

166 

167def quadratic(x): 

168 """A quadratic B-spline. 

169 

170 This is a special case of `bspline`, and equivalent to ``bspline(x, 2)``. 

171 """ 

172 ax = abs(asarray(x)) 

173 res = zeros_like(ax) 

174 cond1 = less(ax, 0.5) 

175 if cond1.any(): 

176 ax1 = ax[cond1] 

177 res[cond1] = 0.75 - ax1 ** 2 

178 cond2 = ~cond1 & less(ax, 1.5) 

179 if cond2.any(): 

180 ax2 = ax[cond2] 

181 res[cond2] = (ax2 - 1.5) ** 2 / 2.0 

182 return res 

183 

184 

185def _coeff_smooth(lam): 

186 xi = 1 - 96 * lam + 24 * lam * sqrt(3 + 144 * lam) 

187 omeg = arctan2(sqrt(144 * lam - 1), sqrt(xi)) 

188 rho = (24 * lam - 1 - sqrt(xi)) / (24 * lam) 

189 rho = rho * sqrt((48 * lam + 24 * lam * sqrt(3 + 144 * lam)) / xi) 

190 return rho, omeg 

191 

192 

193def _hc(k, cs, rho, omega): 

194 return (cs / sin(omega) * (rho ** k) * sin(omega * (k + 1)) * 

195 greater(k, -1)) 

196 

197 

198def _hs(k, cs, rho, omega): 

199 c0 = (cs * cs * (1 + rho * rho) / (1 - rho * rho) / 

200 (1 - 2 * rho * rho * cos(2 * omega) + rho ** 4)) 

201 gamma = (1 - rho * rho) / (1 + rho * rho) / tan(omega) 

202 ak = abs(k) 

203 return c0 * rho ** ak * (cos(omega * ak) + gamma * sin(omega * ak)) 

204 

205 

206def _cubic_smooth_coeff(signal, lamb): 

207 rho, omega = _coeff_smooth(lamb) 

208 cs = 1 - 2 * rho * cos(omega) + rho * rho 

209 K = len(signal) 

210 yp = zeros((K,), signal.dtype.char) 

211 k = arange(K) 

212 yp[0] = (_hc(0, cs, rho, omega) * signal[0] + 

213 add.reduce(_hc(k + 1, cs, rho, omega) * signal)) 

214 

215 yp[1] = (_hc(0, cs, rho, omega) * signal[0] + 

216 _hc(1, cs, rho, omega) * signal[1] + 

217 add.reduce(_hc(k + 2, cs, rho, omega) * signal)) 

218 

219 for n in range(2, K): 

220 yp[n] = (cs * signal[n] + 2 * rho * cos(omega) * yp[n - 1] - 

221 rho * rho * yp[n - 2]) 

222 

223 y = zeros((K,), signal.dtype.char) 

224 

225 y[K - 1] = add.reduce((_hs(k, cs, rho, omega) + 

226 _hs(k + 1, cs, rho, omega)) * signal[::-1]) 

227 y[K - 2] = add.reduce((_hs(k - 1, cs, rho, omega) + 

228 _hs(k + 2, cs, rho, omega)) * signal[::-1]) 

229 

230 for n in range(K - 3, -1, -1): 

231 y[n] = (cs * yp[n] + 2 * rho * cos(omega) * y[n + 1] - 

232 rho * rho * y[n + 2]) 

233 

234 return y 

235 

236 

237def _cubic_coeff(signal): 

238 zi = -2 + sqrt(3) 

239 K = len(signal) 

240 yplus = zeros((K,), signal.dtype.char) 

241 powers = zi ** arange(K) 

242 yplus[0] = signal[0] + zi * add.reduce(powers * signal) 

243 for k in range(1, K): 

244 yplus[k] = signal[k] + zi * yplus[k - 1] 

245 output = zeros((K,), signal.dtype) 

246 output[K - 1] = zi / (zi - 1) * yplus[K - 1] 

247 for k in range(K - 2, -1, -1): 

248 output[k] = zi * (output[k + 1] - yplus[k]) 

249 return output * 6.0 

250 

251 

252def _quadratic_coeff(signal): 

253 zi = -3 + 2 * sqrt(2.0) 

254 K = len(signal) 

255 yplus = zeros((K,), signal.dtype.char) 

256 powers = zi ** arange(K) 

257 yplus[0] = signal[0] + zi * add.reduce(powers * signal) 

258 for k in range(1, K): 

259 yplus[k] = signal[k] + zi * yplus[k - 1] 

260 output = zeros((K,), signal.dtype.char) 

261 output[K - 1] = zi / (zi - 1) * yplus[K - 1] 

262 for k in range(K - 2, -1, -1): 

263 output[k] = zi * (output[k + 1] - yplus[k]) 

264 return output * 8.0 

265 

266 

267def cspline1d(signal, lamb=0.0): 

268 """ 

269 Compute cubic spline coefficients for rank-1 array. 

270 

271 Find the cubic spline coefficients for a 1-D signal assuming 

272 mirror-symmetric boundary conditions. To obtain the signal back from the 

273 spline representation mirror-symmetric-convolve these coefficients with a 

274 length 3 FIR window [1.0, 4.0, 1.0]/ 6.0 . 

275 

276 Parameters 

277 ---------- 

278 signal : ndarray 

279 A rank-1 array representing samples of a signal. 

280 lamb : float, optional 

281 Smoothing coefficient, default is 0.0. 

282 

283 Returns 

284 ------- 

285 c : ndarray 

286 Cubic spline coefficients. 

287 

288 """ 

289 if lamb != 0.0: 

290 return _cubic_smooth_coeff(signal, lamb) 

291 else: 

292 return _cubic_coeff(signal) 

293 

294 

295def qspline1d(signal, lamb=0.0): 

296 """Compute quadratic spline coefficients for rank-1 array. 

297 

298 Parameters 

299 ---------- 

300 signal : ndarray 

301 A rank-1 array representing samples of a signal. 

302 lamb : float, optional 

303 Smoothing coefficient (must be zero for now). 

304 

305 Returns 

306 ------- 

307 c : ndarray 

308 Quadratic spline coefficients. 

309 

310 See Also 

311 -------- 

312 qspline1d_eval : Evaluate a quadratic spline at the new set of points. 

313 

314 Notes 

315 ----- 

316 Find the quadratic spline coefficients for a 1-D signal assuming 

317 mirror-symmetric boundary conditions. To obtain the signal back from the 

318 spline representation mirror-symmetric-convolve these coefficients with a 

319 length 3 FIR window [1.0, 6.0, 1.0]/ 8.0 . 

320 

321 Examples 

322 -------- 

323 We can filter a signal to reduce and smooth out high-frequency noise with 

324 a quadratic spline: 

325 

326 >>> import matplotlib.pyplot as plt 

327 >>> from scipy.signal import qspline1d, qspline1d_eval 

328 >>> sig = np.repeat([0., 1., 0.], 100) 

329 >>> sig += np.random.randn(len(sig))*0.05 # add noise 

330 >>> time = np.linspace(0, len(sig)) 

331 >>> filtered = qspline1d_eval(qspline1d(sig), time) 

332 >>> plt.plot(sig, label="signal") 

333 >>> plt.plot(time, filtered, label="filtered") 

334 >>> plt.legend() 

335 >>> plt.show() 

336 

337 """ 

338 if lamb != 0.0: 

339 raise ValueError("Smoothing quadratic splines not supported yet.") 

340 else: 

341 return _quadratic_coeff(signal) 

342 

343 

344def cspline1d_eval(cj, newx, dx=1.0, x0=0): 

345 """Evaluate a spline at the new set of points. 

346 

347 `dx` is the old sample-spacing while `x0` was the old origin. In 

348 other-words the old-sample points (knot-points) for which the `cj` 

349 represent spline coefficients were at equally-spaced points of: 

350 

351 oldx = x0 + j*dx j=0...N-1, with N=len(cj) 

352 

353 Edges are handled using mirror-symmetric boundary conditions. 

354 

355 """ 

356 newx = (asarray(newx) - x0) / float(dx) 

357 res = zeros_like(newx, dtype=cj.dtype) 

358 if res.size == 0: 

359 return res 

360 N = len(cj) 

361 cond1 = newx < 0 

362 cond2 = newx > (N - 1) 

363 cond3 = ~(cond1 | cond2) 

364 # handle general mirror-symmetry 

365 res[cond1] = cspline1d_eval(cj, -newx[cond1]) 

366 res[cond2] = cspline1d_eval(cj, 2 * (N - 1) - newx[cond2]) 

367 newx = newx[cond3] 

368 if newx.size == 0: 

369 return res 

370 result = zeros_like(newx, dtype=cj.dtype) 

371 jlower = floor(newx - 2).astype(int) + 1 

372 for i in range(4): 

373 thisj = jlower + i 

374 indj = thisj.clip(0, N - 1) # handle edge cases 

375 result += cj[indj] * cubic(newx - thisj) 

376 res[cond3] = result 

377 return res 

378 

379 

380def qspline1d_eval(cj, newx, dx=1.0, x0=0): 

381 """Evaluate a quadratic spline at the new set of points. 

382 

383 Parameters 

384 ---------- 

385 cj : ndarray 

386 Quadratic spline coefficients 

387 newx : ndarray 

388 New set of points. 

389 dx : float, optional 

390 Old sample-spacing, the default value is 1.0. 

391 x0 : int, optional 

392 Old origin, the default value is 0. 

393 

394 Returns 

395 ------- 

396 res : ndarray 

397 Evaluated a quadratic spline points. 

398 

399 See Also 

400 -------- 

401 qspline1d : Compute quadratic spline coefficients for rank-1 array. 

402 

403 Notes 

404 ----- 

405 `dx` is the old sample-spacing while `x0` was the old origin. In 

406 other-words the old-sample points (knot-points) for which the `cj` 

407 represent spline coefficients were at equally-spaced points of:: 

408 

409 oldx = x0 + j*dx j=0...N-1, with N=len(cj) 

410 

411 Edges are handled using mirror-symmetric boundary conditions. 

412 

413 Examples 

414 -------- 

415 We can filter a signal to reduce and smooth out high-frequency noise with 

416 a quadratic spline: 

417 

418 >>> import matplotlib.pyplot as plt 

419 >>> from scipy.signal import qspline1d, qspline1d_eval 

420 >>> sig = np.repeat([0., 1., 0.], 100) 

421 >>> sig += np.random.randn(len(sig))*0.05 # add noise 

422 >>> time = np.linspace(0, len(sig)) 

423 >>> filtered = qspline1d_eval(qspline1d(sig), time) 

424 >>> plt.plot(sig, label="signal") 

425 >>> plt.plot(time, filtered, label="filtered") 

426 >>> plt.legend() 

427 >>> plt.show() 

428 

429 """ 

430 newx = (asarray(newx) - x0) / dx 

431 res = zeros_like(newx) 

432 if res.size == 0: 

433 return res 

434 N = len(cj) 

435 cond1 = newx < 0 

436 cond2 = newx > (N - 1) 

437 cond3 = ~(cond1 | cond2) 

438 # handle general mirror-symmetry 

439 res[cond1] = qspline1d_eval(cj, -newx[cond1]) 

440 res[cond2] = qspline1d_eval(cj, 2 * (N - 1) - newx[cond2]) 

441 newx = newx[cond3] 

442 if newx.size == 0: 

443 return res 

444 result = zeros_like(newx) 

445 jlower = floor(newx - 1.5).astype(int) + 1 

446 for i in range(3): 

447 thisj = jlower + i 

448 indj = thisj.clip(0, N - 1) # handle edge cases 

449 result += cj[indj] * quadratic(newx - thisj) 

450 res[cond3] = result 

451 return res