python自动重采样数据

前端之家收集整理的这篇文章主要介绍了python自动重采样数据 前端之家小编觉得挺不错的,现在分享给大家,也给大家做个参考。

假设我有以下数据(测量值):

enter image description here

如您所见,有很多尖点(即坡度变化很大的地方).因此,最好在这些点周围再进行一些测量.为此,我编写了一个脚本:

>我计算了三个连续点的曲率:
Menger曲率:https://en.wikipedia.org/wiki/Menger_curvature#Definition
>然后,根据曲率决定应重新采样哪些值.

…然后迭代直到平均曲率下降…但是它不起作用,因为它上升了.你知道为什么吗 ?

这是完整的代码(在x值的长度达到60后停止它):

  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3. def curvature(A,B,C):
  4. """Calculates the Menger curvature fro three Points,given as numpy arrays.
  5. Sources:
  6. Menger curvature: https://en.wikipedia.org/wiki/Menger_curvature#Definition
  7. Area of a triangle given 3 points: https://math.stackexchange.com/questions/516219/finding-out-the-area-of-a-triangle-if-the-coordinates-of-the-three-vertices-are
  8. """
  9. # Pre-check: Making sure that the input points are all numpy arrays
  10. if any(x is not np.ndarray for x in [type(A),type(B),type(C)]):
  11. print("The input points need to be a numpy array,currently it is a ",type(A))
  12. # Augment Columns
  13. A_aug = np.append(A,1)
  14. B_aug = np.append(B,1)
  15. C_aug = np.append(C,1)
  16. # Caclulate Area of Triangle
  17. matrix = np.column_stack((A_aug,B_aug,C_aug))
  18. area = 1/2*np.linalg.det(matrix)
  19. # Special case: Two or more points are equal
  20. if np.all(A == B) or np.all(B == C):
  21. curvature = 0
  22. else:
  23. curvature = 4*area/(np.linalg.norm(A-B)*np.linalg.norm(B-C)*np.linalg.norm(C-A))
  24. # Return Menger curvature
  25. return curvature
  26. def values_to_calulate(x,curvature_list,max_curvature):
  27. """Calculates the new x values which need to be calculated
  28. Middle point between the three points that were used to calculate the curvature """
  29. i = 0
  30. new_x = np.empty(0)
  31. for curvature in curvature_list:
  32. if curvature > max_curvature:
  33. new_x = np.append(new_x,x[i]+(x[i+2]-x[i])/3 )
  34. i = i+1
  35. return new_x
  36. def plot(x,y,title,xLabel,yLabel):
  37. """Just to visualize"""
  38. # Plot
  39. plt.scatter(x,y)
  40. plt.plot(x,'-o')
  41. # Give a title for the sine wave plot
  42. plt.title(title)
  43. # Give x axis label for the sine wave plot
  44. plt.xlabel(xLabel)
  45. # Give y axis label for the sine wave plot
  46. plt.ylabel(yLabel)
  47. plt.grid(True,which='both')
  48. plt.axhline(y=0,color='k')
  49. # Display the sine wave
  50. plt.show
  51. plt.pause(0.05)
  52. ### STARTS HERE
  53. # Get x values of the sine wave
  54. x = np.arange(0,10,1);
  55. # Amplitude of the sine wave is sine of a variable like time
  56. def function(x):
  57. return 1+np.sin(x)*np.cos(x)**2
  58. y = function(x)
  59. # Plot it
  60. plot(x,title='Data',xLabel='Time',yLabel='Amplitude')
  61. continue_Loop = True
  62. while continue_Loop == True :
  63. curvature_list = np.empty(0)
  64. for i in range(len(x)-2):
  65. # Get the three points
  66. A = np.array([x[i],y[i]])
  67. B = np.array([x[i+1],y[i+1]])
  68. C = np.array([x[i+2],y[i+2]])
  69. # Calculate the curvature
  70. curvature_value = abs(curvature(A,C))
  71. curvature_list = np.append(curvature_list,curvature_value)
  72. print("len: ",len(x) )
  73. print("average curvature: ",np.average(curvature_list))
  74. # Calculate the points that need to be added
  75. x_new = values_to_calulate(x,max_curvature=0.3)
  76. # Add those values to the current x list:
  77. x = np.sort(np.append(x,x_new))
  78. # STOPED IT AFTER len(x) == 60
  79. if len(x) >= 60:
  80. continue_Loop = False
  81. # Amplitude of the sine wave is sine of a variable like time
  82. y = function(x)
  83. # Plot it
  84. plot(x,yLabel='Amplitude')

它应该是这样的:

enter image description here

编辑:

如果让它继续运行…:

enter image description here

最佳答案
因此,总结以上我的评论

>您正在计算曲线的平均曲率,没有理由将其设为0.在每个点上,无论您的点有多近,圆半径都会收敛到该点处的曲率,而不是0.
>一种替代方法是使用两点之间的绝对导数变化:保持采样直到abs(d(df / dx))<. some_threshold其中d(df / dx)=(df / dx)[n]-(df / dx)[n-1]

猜你在找的Python相关文章