回答編集履歴

説明の追記

2020/12/05 02:33

投稿

jbpb0

スコア7651

test CHANGED Viewed

@@ -52,7 +52,7 @@
 [Savitzky-Golay smoothing filter for not equally spaced data](https://dsp.stackexchange.com/questions/1676/savitzky-golay-smoothing-filter-for-not-equally-spaced-data)
-に、不等間隔でも計算できるSavitzky–Golay法のPythonコードが載ってたので、それを微分にも使えるように修正しました
+に、不等間隔(non-uniform)でも計算できるSavitzky–Golay法のPythonコードが載ってたので、それを微分(differential)にも使えるように修正しました
 よろしければお試しください

コードのバグ修正

2020/12/05 02:33

投稿

jbpb0

スコア7651

test CHANGED Viewed

@@ -60,6 +60,8 @@
 import numpy as np
+import math
 def non_uniform_savgol_der(x, y, window, polynom, der):
@@ -102,7 +104,7 @@
     der : int
-        The order of derivative. Must be positive and less than 3
+        The order of derivative. Must be positive and smaller than polynom order
@@ -152,9 +154,9 @@
-    if polynom < 2:
+    if polynom < 0:
-        raise ValueError('"polynom" must be larger than 1')
+        raise ValueError('"polynom" must be larger than 0')
@@ -170,9 +172,9 @@
-    if der > 2:
+    if der > polynom:
-        raise TypeError('"der" must be less than 3')
+        raise TypeError('"der" must be equal or less than "polynom"')
@@ -248,7 +250,7 @@
             #y_smoothed[i] += coeffs[0, j] * y[i + j - half_window]
-            y_smoothed[i] += coeffs[der, j] * y[i + j - half_window]
+            y_smoothed[i] += coeffs[der, j] * y[i + j - half_window] * math.factorial(der)
@@ -288,7 +290,7 @@
         for j in range(der, polynom, 1):
-            y_smoothed[i] += first_coeffs[j] * x_i
+            y_smoothed[i] += first_coeffs[j] * x_i * math.factorial(j)
             x_i *= x[i] - x[half_window]
@@ -306,7 +308,7 @@
         for j in range(der, polynom, 1):
-            y_smoothed[i] += last_coeffs[j] * x_i
+            y_smoothed[i] += last_coeffs[j] * x_i * math.factorial(j)
             x_i *= x[i] - x[-half_window - 1]
@@ -326,9 +328,9 @@
     # テスト用ダミーデータ
-    x = np.arange(0, 100)
+    x = np.arange(0, 10, 0.1)
-    y = np.arange(10, 110) + np.random.randn(100) * 2
+    y = 2 * x * x + 3 * x + 4 + np.random.randn(100) / 10
     y11 = non_uniform_savgol_der(x, y, 11, 2, 0)

コード追加

2020/12/02 14:14

投稿

jbpb0

スコア7651

test CHANGED Viewed

@@ -367,3 +367,49 @@
     plt.show()
 ```
+xが順番に並んでないのは、sortすれば直せます
+```python
+# テスト用ダミーデータ
+x = np.array([1, 2, 3, 5, 4, 6])
+y = x * 2
+print(x)
+print(y)
+# ソート (xの順番でyも)
+tmp = zip(x, y)
+tmp2 = sorted(tmp)
+xx, yy = zip(*tmp2)
+xx = np.array(xx)
+yy = np.array(yy)
+# 比較
+print(x)
+print(xx)
+print(y)
+print(yy)
+```

コード追加

2020/12/02 04:58

投稿

jbpb0

スコア7651

test CHANGED Viewed

@@ -45,3 +45,325 @@
 データ見たら等間隔ではないので、使えませんね
 失礼しました
+【追記】
+[Savitzky-Golay smoothing filter for not equally spaced data](https://dsp.stackexchange.com/questions/1676/savitzky-golay-smoothing-filter-for-not-equally-spaced-data)
+に、不等間隔でも計算できるSavitzky–Golay法のPythonコードが載ってたので、それを微分にも使えるように修正しました
+よろしければお試しください
+```python
+import numpy as np
+def non_uniform_savgol_der(x, y, window, polynom, der):
+    """
+    Applies a Savitzky-Golay filter to y with non-uniform spacing
+    as defined in x
+    This is based on https://dsp.stackexchange.com/questions/1676/savitzky-golay-smoothing-filter-for-not-equally-spaced-data
+    The borders are interpolated like scipy.signal.savgol_filter would do
+    Parameters
+    ----------
+    x : array_like
+        List of floats representing the x values of the data
+    y : array_like
+        List of floats representing the y values. Must have same length
+        as x
+    window : int (odd)
+        Window length of datapoints. Must be odd and smaller than x
+    polynom : int
+        The order of polynom used. Must be smaller than the window size
+    der : int
+        The order of derivative. Must be positive and less than 3
+    Returns
+    -------
+    np.array of float
+        The smoothed y values
+    """
+    if len(x) != len(y):
+        raise ValueError('"x" and "y" must be of the same size')
+    if len(x) < window:
+        raise ValueError('The data size must be larger than the window size')
+    if type(window) is not int:
+        raise TypeError('"window" must be an integer')
+    if window % 2 == 0:
+        raise ValueError('The "window" must be an odd integer')
+    if type(polynom) is not int:
+        raise TypeError('"polynom" must be an integer')
+    if polynom >= window:
+        raise ValueError('"polynom" must be less than "window"')
+    if polynom < 2:
+        raise ValueError('"polynom" must be larger than 1')
+    if type(der) is not int:
+        raise TypeError('"der" must be an integer')
+    if der < 0:
+        raise TypeError('"der" must be an positive integer')
+    if der > 2:
+        raise TypeError('"der" must be less than 3')
+    half_window = window // 2
+    polynom += 1
+    # Initialize variables
+    A = np.empty((window, polynom))     # Matrix
+    tA = np.empty((polynom, window))    # Transposed matrix
+    t = np.empty(window)                # Local x variables
+    y_smoothed = np.full(len(y), np.nan)
+    # Start smoothing
+    for i in range(half_window, len(x) - half_window, 1):
+        # Center a window of x values on x[i]
+        for j in range(0, window, 1):
+            t[j] = x[i + j - half_window] - x[i]
+        # Create the initial matrix A and its transposed form tA
+        for j in range(0, window, 1):
+            r = 1.0
+            for k in range(0, polynom, 1):
+                A[j, k] = r
+                tA[k, j] = r
+                r *= t[j]
+        # Multiply the two matrices
+        tAA = np.matmul(tA, A)
+        # Invert the product of the matrices
+        tAA = np.linalg.inv(tAA)
+        # Calculate the pseudoinverse of the design matrix
+        coeffs = np.matmul(tAA, tA)
+        # Calculate c0 which is also the y value for y[i]
+        y_smoothed[i] = 0
+        for j in range(0, window, 1):
+            #y_smoothed[i] += coeffs[0, j] * y[i + j - half_window]
+            y_smoothed[i] += coeffs[der, j] * y[i + j - half_window]
+        # If at the end or beginning, store all coefficients for the polynom
+        if i == half_window:
+            first_coeffs = np.zeros(polynom)
+            for j in range(0, window, 1):
+                for k in range(polynom):
+                    first_coeffs[k] += coeffs[k, j] * y[j]
+        elif i == len(x) - half_window - 1:
+            last_coeffs = np.zeros(polynom)
+            for j in range(0, window, 1):
+                for k in range(polynom):
+                    last_coeffs[k] += coeffs[k, j] * y[len(y) - window + j]
+    # Interpolate the result at the left border
+    for i in range(0, half_window, 1):
+        y_smoothed[i] = 0
+        x_i = 1
+        #for j in range(0, polynom, 1):
+        for j in range(der, polynom, 1):
+            y_smoothed[i] += first_coeffs[j] * x_i
+            x_i *= x[i] - x[half_window]
+    # Interpolate the result at the right border
+    for i in range(len(x) - half_window, len(x), 1):
+        y_smoothed[i] = 0
+        x_i = 1
+        #for j in range(0, polynom, 1):
+        for j in range(der, polynom, 1):
+            y_smoothed[i] += last_coeffs[j] * x_i
+            x_i *= x[i] - x[-half_window - 1]
+    return y_smoothed
+if __name__ == '__main__':
+    import numpy as np
+    import matplotlib.pyplot as plt
+    # テスト用ダミーデータ
+    x = np.arange(0, 100)
+    y = np.arange(10, 110) + np.random.randn(100) * 2
+    y11 = non_uniform_savgol_der(x, y, 11, 2, 0)
+    plt.plot(x, y, x, y11)
+    plt.grid()
+    plt.show()
+    # 微分
+    dy = np.gradient(y, x)
+    dy11 = non_uniform_savgol_der(x, y, 11, 2, 1)
+    plt.plot(x, dy, x, dy11)
+    plt.grid()
+    plt.show()
+    # 二階微分
+    ddy = np.gradient(dy, x)
+    ddy11 = non_uniform_savgol_der(x, y, 11, 2, 2)
+    plt.plot(x, ddy, x, ddy11)
+    plt.grid()
+    plt.show()
+```

誤りの告知

2020/12/02 04:31

投稿

jbpb0

スコア7651

test CHANGED Viewed

@@ -35,3 +35,13 @@
 ```
 点数が多いほど平滑化が強くなります
+【訂正】
+これはxが等間隔じゃないと使えない手法です
+データ見たら等間隔ではないので、使えませんね
+失礼しました