test_bsplines.py 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224
  1. # pylint: disable=missing-docstring
  2. from __future__ import division, print_function, absolute_import
  3. import numpy as np
  4. from numpy import array
  5. from numpy.testing import (TestCase, assert_equal, run_module_suite,
  6. assert_allclose, assert_array_equal,
  7. assert_almost_equal)
  8. from pytest import raises
  9. import scipy.signal.bsplines as bsp
  10. class TestBSplines(object):
  11. """Test behaviors of bsplines. The values tested against were returned as of
  12. scipy 1.1.0 and are included for regression testing purposes"""
  13. def test_factorial(self):
  14. # can't all be zero state
  15. assert_equal(bsp.factorial(1), 1)
  16. def test_spline_filter(self):
  17. np.random.seed(12457)
  18. # Test the type-error branch
  19. raises(TypeError, bsp.spline_filter, array([0]), 0)
  20. # Test the complex branch
  21. data_array_complex = np.random.rand(7, 7) + np.random.rand(7, 7)*1j
  22. # make the magnitude exceed 1, and make some negative
  23. data_array_complex = 10*(1+1j-2*data_array_complex)
  24. result_array_complex = array(
  25. [[-4.61489230e-01-1.92994022j, 8.33332443+6.25519943j,
  26. 6.96300745e-01-9.05576038j, 5.28294849+3.97541356j,
  27. 5.92165565+7.68240595j, 6.59493160-1.04542804j,
  28. 9.84503460-5.85946894j],
  29. [-8.78262329-8.4295969j, 7.20675516+5.47528982j,
  30. -8.17223072+2.06330729j, -4.38633347-8.65968037j,
  31. 9.89916801-8.91720295j, 2.67755103+8.8706522j,
  32. 6.24192142+3.76879835j],
  33. [-3.15627527+2.56303072j, 9.87658501-0.82838702j,
  34. -9.96930313+8.72288895j, 3.17193985+6.42474651j,
  35. -4.50919819-6.84576082j, 5.75423431+9.94723988j,
  36. 9.65979767+6.90665293j],
  37. [-8.28993416-6.61064005j, 9.71416473e-01-9.44907284j,
  38. -2.38331890+9.25196648j, -7.08868170-0.77403212j,
  39. 4.89887714+7.05371094j, -1.37062311-2.73505688j,
  40. 7.70705748+2.5395329j],
  41. [2.51528406-1.82964492j, 3.65885472+2.95454836j,
  42. 5.16786575-1.66362023j, -8.77737999e-03+5.72478867j,
  43. 4.10533333-3.10287571j, 9.04761887+1.54017115j,
  44. -5.77960968e-01-7.87758923j],
  45. [9.86398506-3.98528528j, -4.71444130-2.44316983j,
  46. -1.68038976-1.12708664j, 2.84695053+1.01725709j,
  47. 1.14315915-8.89294529j, -3.17127085-5.42145538j,
  48. 1.91830420-6.16370344j],
  49. [7.13875294+2.91851187j, -5.35737514+9.64132309j,
  50. -9.66586399+0.70250005j, -9.87717438-2.0262239j,
  51. 9.93160629+1.5630846j, 4.71948051-2.22050714j,
  52. 9.49550819+7.8995142j]])
  53. # FIXME: for complex types, the computations are done in
  54. # single precision (reason unclear). When this is changed,
  55. # this test needs updating.
  56. assert_allclose(bsp.spline_filter(data_array_complex, 0),
  57. result_array_complex, rtol=1e-6)
  58. # Test the real branch
  59. np.random.seed(12457)
  60. data_array_real = np.random.rand(12, 12)
  61. # make the magnitude exceed 1, and make some negative
  62. data_array_real = 10*(1-2*data_array_real)
  63. result_array_real = array(
  64. [[-.463312621, 8.33391222, .697290949, 5.28390836,
  65. 5.92066474, 6.59452137, 9.84406950, -8.78324188,
  66. 7.20675750, -8.17222994, -4.38633345, 9.89917069],
  67. [2.67755154, 6.24192170, -3.15730578, 9.87658581,
  68. -9.96930425, 3.17194115, -4.50919947, 5.75423446,
  69. 9.65979824, -8.29066885, .971416087, -2.38331897],
  70. [-7.08868346, 4.89887705, -1.37062289, 7.70705838,
  71. 2.51526461, 3.65885497, 5.16786604, -8.77715342e-03,
  72. 4.10533325, 9.04761993, -.577960351, 9.86382519],
  73. [-4.71444301, -1.68038985, 2.84695116, 1.14315938,
  74. -3.17127091, 1.91830461, 7.13779687, -5.35737482,
  75. -9.66586425, -9.87717456, 9.93160672, 4.71948144],
  76. [9.49551194, -1.92958436, 6.25427993, -9.05582911,
  77. 3.97562282, 7.68232426, -1.04514824, -5.86021443,
  78. -8.43007451, 5.47528997, 2.06330736, -8.65968112],
  79. [-8.91720100, 8.87065356, 3.76879937, 2.56222894,
  80. -.828387146, 8.72288903, 6.42474741, -6.84576083,
  81. 9.94724115, 6.90665380, -6.61084494, -9.44907391],
  82. [9.25196790, -.774032030, 7.05371046, -2.73505725,
  83. 2.53953305, -1.82889155, 2.95454824, -1.66362046,
  84. 5.72478916, -3.10287679, 1.54017123, -7.87759020],
  85. [-3.98464539, -2.44316992, -1.12708657, 1.01725672,
  86. -8.89294671, -5.42145629, -6.16370321, 2.91775492,
  87. 9.64132208, .702499998, -2.02622392, 1.56308431],
  88. [-2.22050773, 7.89951554, 5.98970713, -7.35861835,
  89. 5.45459283, -7.76427957, 3.67280490, -4.05521315,
  90. 4.51967507, -3.22738749, -3.65080177, 3.05630155],
  91. [-6.21240584, -.296796126, -8.34800163, 9.21564563,
  92. -3.61958784, -4.77120006, -3.99454057, 1.05021988e-03,
  93. -6.95982829, 6.04380797, 8.43181250, -2.71653339],
  94. [1.19638037, 6.99718842e-02, 6.72020394, -2.13963198,
  95. 3.75309875, -5.70076744, 5.92143551, -7.22150575,
  96. -3.77114594, -1.11903194, -5.39151466, 3.06620093],
  97. [9.86326886, 1.05134482, -7.75950607, -3.64429655,
  98. 7.81848957, -9.02270373, 3.73399754, -4.71962549,
  99. -7.71144306, 3.78263161, 6.46034818, -4.43444731]])
  100. assert_allclose(bsp.spline_filter(data_array_real, 0),
  101. result_array_real)
  102. def test_bspline(self):
  103. np.random.seed(12458)
  104. assert_allclose(bsp.bspline(np.random.rand(1, 1), 2),
  105. array([[0.73694695]]))
  106. data_array_complex = np.random.rand(4, 4) + np.random.rand(4, 4)*1j
  107. data_array_complex = 0.1*data_array_complex
  108. result_array_complex = array(
  109. [[0.40882362, 0.41021151, 0.40886708, 0.40905103],
  110. [0.40829477, 0.41021230, 0.40966097, 0.40939871],
  111. [0.41036803, 0.40901724, 0.40965331, 0.40879513],
  112. [0.41032862, 0.40925287, 0.41037754, 0.41027477]])
  113. assert_allclose(bsp.bspline(data_array_complex, 10),
  114. result_array_complex)
  115. def test_gauss_spline(self):
  116. np.random.seed(12459)
  117. assert_almost_equal(bsp.gauss_spline(0, 0), 1.381976597885342)
  118. assert_allclose(bsp.gauss_spline(array([1.]), 1), array([0.04865217]))
  119. def test_cubic(self):
  120. np.random.seed(12460)
  121. assert_array_equal(bsp.cubic([0]), array([0]))
  122. data_array_complex = np.random.rand(4, 4) + np.random.rand(4, 4)*1j
  123. data_array_complex = 1+1j-2*data_array_complex
  124. # scaling the magnitude by 10 makes the results close enough to zero,
  125. # that the assertion fails, so just make the elements have a mix of
  126. # positive and negative imaginary components...
  127. result_array_complex = array(
  128. [[0.23056563, 0.38414406, 0.08342987, 0.06904847],
  129. [0.17240848, 0.47055447, 0.63896278, 0.39756424],
  130. [0.12672571, 0.65862632, 0.1116695, 0.09700386],
  131. [0.3544116, 0.17856518, 0.1528841, 0.17285762]])
  132. assert_allclose(bsp.cubic(data_array_complex), result_array_complex)
  133. def test_quadratic(self):
  134. np.random.seed(12461)
  135. assert_array_equal(bsp.quadratic([0]), array([0]))
  136. data_array_complex = np.random.rand(4, 4) + np.random.rand(4, 4)*1j
  137. # scaling the magnitude by 10 makes the results all zero,
  138. # so just make the elements have a mix of positive and negative
  139. # imaginary components...
  140. data_array_complex = (1+1j-2*data_array_complex)
  141. result_array_complex = array(
  142. [[0.23062746, 0.06338176, 0.34902312, 0.31944105],
  143. [0.14701256, 0.13277773, 0.29428615, 0.09814697],
  144. [0.52873842, 0.06484157, 0.09517566, 0.46420389],
  145. [0.09286829, 0.09371954, 0.1422526, 0.16007024]])
  146. assert_allclose(bsp.quadratic(data_array_complex),
  147. result_array_complex)
  148. def test_cspline1d(self):
  149. np.random.seed(12462)
  150. assert_array_equal(bsp.cspline1d(array([0])), [0.])
  151. c1d = array([1.21037185, 1.86293902, 2.98834059, 4.11660378,
  152. 4.78893826])
  153. # test lamda != 0
  154. assert_allclose(bsp.cspline1d(array([1., 2, 3, 4, 5]), 1), c1d)
  155. c1d0 = array([0.78683946, 2.05333735, 2.99981113, 3.94741812,
  156. 5.21051638])
  157. assert_allclose(bsp.cspline1d(array([1., 2, 3, 4, 5])), c1d0)
  158. def test_qspline1d(self):
  159. np.random.seed(12463)
  160. assert_array_equal(bsp.qspline1d(array([0])), [0.])
  161. # test lamda != 0
  162. raises(ValueError, bsp.qspline1d, array([1., 2, 3, 4, 5]), 1.)
  163. raises(ValueError, bsp.qspline1d, array([1., 2, 3, 4, 5]), -1.)
  164. q1d0 = array([0.85350007, 2.02441743, 2.99999534, 3.97561055,
  165. 5.14634135])
  166. assert_allclose(bsp.qspline1d(array([1., 2, 3, 4, 5])), q1d0)
  167. def test_cspline1d_eval(self):
  168. np.random.seed(12464)
  169. assert_allclose(bsp.cspline1d_eval(array([0., 0]), [0.]), array([0.]))
  170. assert_array_equal(bsp.cspline1d_eval(array([1., 0, 1]), []),
  171. array([]))
  172. x = [-3, -2, -1, 0, 1, 2, 3, 4, 5, 6]
  173. dx = x[1]-x[0]
  174. newx = [-6., -5.5, -5., -4.5, -4., -3.5, -3., -2.5, -2., -1.5, -1.,
  175. -0.5, 0., 0.5, 1., 1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5., 5.5, 6.,
  176. 6.5, 7., 7.5, 8., 8.5, 9., 9.5, 10., 10.5, 11., 11.5, 12.,
  177. 12.5]
  178. y = array([4.216, 6.864, 3.514, 6.203, 6.759, 7.433, 7.874, 5.879,
  179. 1.396, 4.094])
  180. cj = bsp.cspline1d(y)
  181. newy = array([6.203, 4.41570658, 3.514, 5.16924703, 6.864, 6.04643068,
  182. 4.21600281, 6.04643068, 6.864, 5.16924703, 3.514,
  183. 4.41570658, 6.203, 6.80717667, 6.759, 6.98971173, 7.433,
  184. 7.79560142, 7.874, 7.41525761, 5.879, 3.18686814, 1.396,
  185. 2.24889482, 4.094, 2.24889482, 1.396, 3.18686814, 5.879,
  186. 7.41525761, 7.874, 7.79560142, 7.433, 6.98971173, 6.759,
  187. 6.80717667, 6.203, 4.41570658])
  188. assert_allclose(bsp.cspline1d_eval(cj, newx, dx=dx, x0=x[0]), newy)
  189. def test_qspline1d_eval(self):
  190. np.random.seed(12465)
  191. assert_allclose(bsp.qspline1d_eval(array([0., 0]), [0.]), array([0.]))
  192. assert_array_equal(bsp.qspline1d_eval(array([1., 0, 1]), []),
  193. array([]))
  194. x = [-3, -2, -1, 0, 1, 2, 3, 4, 5, 6]
  195. dx = x[1]-x[0]
  196. newx = [-6., -5.5, -5., -4.5, -4., -3.5, -3., -2.5, -2., -1.5, -1.,
  197. -0.5, 0., 0.5, 1., 1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5., 5.5, 6.,
  198. 6.5, 7., 7.5, 8., 8.5, 9., 9.5, 10., 10.5, 11., 11.5, 12.,
  199. 12.5]
  200. y = array([4.216, 6.864, 3.514, 6.203, 6.759, 7.433, 7.874, 5.879,
  201. 1.396, 4.094])
  202. cj = bsp.qspline1d(y)
  203. newy = array([6.203, 4.49418159, 3.514, 5.18390821, 6.864, 5.91436915,
  204. 4.21600002, 5.91436915, 6.864, 5.18390821, 3.514,
  205. 4.49418159, 6.203, 6.71900226, 6.759, 7.03980488, 7.433,
  206. 7.81016848, 7.874, 7.32718426, 5.879, 3.23872593, 1.396,
  207. 2.34046013, 4.094, 2.34046013, 1.396, 3.23872593, 5.879,
  208. 7.32718426, 7.874, 7.81016848, 7.433, 7.03980488, 6.759,
  209. 6.71900226, 6.203, 4.49418159])
  210. assert_allclose(bsp.qspline1d_eval(cj, newx, dx=dx, x0=x[0]), newy)