回答編集履歴
5
説明の追記
answer
CHANGED
|
@@ -25,7 +25,7 @@
|
|
|
25
25
|
|
|
26
26
|
【追記】
|
|
27
27
|
[Savitzky-Golay smoothing filter for not equally spaced data](https://dsp.stackexchange.com/questions/1676/savitzky-golay-smoothing-filter-for-not-equally-spaced-data)
|
|
28
|
-
に、不等間隔でも計算できるSavitzky–Golay法のPythonコードが載ってたので、それを微分にも使えるように修正しました
|
|
28
|
+
に、不等間隔(non-uniform)でも計算できるSavitzky–Golay法のPythonコードが載ってたので、それを微分(differential)にも使えるように修正しました
|
|
29
29
|
よろしければお試しください
|
|
30
30
|
```python
|
|
31
31
|
import numpy as np
|
4
コードのバグ修正
answer
CHANGED
|
@@ -29,6 +29,7 @@
|
|
|
29
29
|
よろしければお試しください
|
|
30
30
|
```python
|
|
31
31
|
import numpy as np
|
|
32
|
+
import math
|
|
32
33
|
|
|
33
34
|
def non_uniform_savgol_der(x, y, window, polynom, der):
|
|
34
35
|
"""
|
|
@@ -50,7 +51,7 @@
|
|
|
50
51
|
polynom : int
|
|
51
52
|
The order of polynom used. Must be smaller than the window size
|
|
52
53
|
der : int
|
|
53
|
-
The order of derivative. Must be positive and
|
|
54
|
+
The order of derivative. Must be positive and smaller than polynom order
|
|
54
55
|
|
|
55
56
|
Returns
|
|
56
57
|
-------
|
|
@@ -75,8 +76,8 @@
|
|
|
75
76
|
if polynom >= window:
|
|
76
77
|
raise ValueError('"polynom" must be less than "window"')
|
|
77
78
|
|
|
78
|
-
if polynom <
|
|
79
|
+
if polynom < 0:
|
|
79
|
-
raise ValueError('"polynom" must be larger than
|
|
80
|
+
raise ValueError('"polynom" must be larger than 0')
|
|
80
81
|
|
|
81
82
|
if type(der) is not int:
|
|
82
83
|
raise TypeError('"der" must be an integer')
|
|
@@ -84,8 +85,8 @@
|
|
|
84
85
|
if der < 0:
|
|
85
86
|
raise TypeError('"der" must be an positive integer')
|
|
86
87
|
|
|
87
|
-
if der >
|
|
88
|
+
if der > polynom:
|
|
88
|
-
raise TypeError('"der" must be less than
|
|
89
|
+
raise TypeError('"der" must be equal or less than "polynom"')
|
|
89
90
|
|
|
90
91
|
half_window = window // 2
|
|
91
92
|
polynom += 1
|
|
@@ -123,7 +124,7 @@
|
|
|
123
124
|
y_smoothed[i] = 0
|
|
124
125
|
for j in range(0, window, 1):
|
|
125
126
|
#y_smoothed[i] += coeffs[0, j] * y[i + j - half_window]
|
|
126
|
-
y_smoothed[i] += coeffs[der, j] * y[i + j - half_window]
|
|
127
|
+
y_smoothed[i] += coeffs[der, j] * y[i + j - half_window] * math.factorial(der)
|
|
127
128
|
|
|
128
129
|
# If at the end or beginning, store all coefficients for the polynom
|
|
129
130
|
if i == half_window:
|
|
@@ -143,7 +144,7 @@
|
|
|
143
144
|
x_i = 1
|
|
144
145
|
#for j in range(0, polynom, 1):
|
|
145
146
|
for j in range(der, polynom, 1):
|
|
146
|
-
y_smoothed[i] += first_coeffs[j] * x_i
|
|
147
|
+
y_smoothed[i] += first_coeffs[j] * x_i * math.factorial(j)
|
|
147
148
|
x_i *= x[i] - x[half_window]
|
|
148
149
|
|
|
149
150
|
# Interpolate the result at the right border
|
|
@@ -152,7 +153,7 @@
|
|
|
152
153
|
x_i = 1
|
|
153
154
|
#for j in range(0, polynom, 1):
|
|
154
155
|
for j in range(der, polynom, 1):
|
|
155
|
-
y_smoothed[i] += last_coeffs[j] * x_i
|
|
156
|
+
y_smoothed[i] += last_coeffs[j] * x_i * math.factorial(j)
|
|
156
157
|
x_i *= x[i] - x[-half_window - 1]
|
|
157
158
|
|
|
158
159
|
return y_smoothed
|
|
@@ -162,8 +163,8 @@
|
|
|
162
163
|
import matplotlib.pyplot as plt
|
|
163
164
|
|
|
164
165
|
# テスト用ダミーデータ
|
|
165
|
-
x = np.arange(0,
|
|
166
|
+
x = np.arange(0, 10, 0.1)
|
|
166
|
-
y =
|
|
167
|
+
y = 2 * x * x + 3 * x + 4 + np.random.randn(100) / 10
|
|
167
168
|
y11 = non_uniform_savgol_der(x, y, 11, 2, 0)
|
|
168
169
|
plt.plot(x, y, x, y11)
|
|
169
170
|
plt.grid()
|
3
コード追加
answer
CHANGED
|
@@ -182,4 +182,27 @@
|
|
|
182
182
|
plt.plot(x, ddy, x, ddy11)
|
|
183
183
|
plt.grid()
|
|
184
184
|
plt.show()
|
|
185
|
+
```
|
|
186
|
+
|
|
187
|
+
|
|
188
|
+
xが順番に並んでないのは、sortすれば直せます
|
|
189
|
+
```python
|
|
190
|
+
# テスト用ダミーデータ
|
|
191
|
+
x = np.array([1, 2, 3, 5, 4, 6])
|
|
192
|
+
y = x * 2
|
|
193
|
+
print(x)
|
|
194
|
+
print(y)
|
|
195
|
+
|
|
196
|
+
# ソート (xの順番でyも)
|
|
197
|
+
tmp = zip(x, y)
|
|
198
|
+
tmp2 = sorted(tmp)
|
|
199
|
+
xx, yy = zip(*tmp2)
|
|
200
|
+
xx = np.array(xx)
|
|
201
|
+
yy = np.array(yy)
|
|
202
|
+
|
|
203
|
+
# 比較
|
|
204
|
+
print(x)
|
|
205
|
+
print(xx)
|
|
206
|
+
print(y)
|
|
207
|
+
print(yy)
|
|
185
208
|
```
|
2
コード追加
answer
CHANGED
|
@@ -21,4 +21,165 @@
|
|
|
21
21
|
【訂正】
|
|
22
22
|
これはxが等間隔じゃないと使えない手法です
|
|
23
23
|
データ見たら等間隔ではないので、使えませんね
|
|
24
|
-
失礼しました
|
|
24
|
+
失礼しました
|
|
25
|
+
|
|
26
|
+
【追記】
|
|
27
|
+
[Savitzky-Golay smoothing filter for not equally spaced data](https://dsp.stackexchange.com/questions/1676/savitzky-golay-smoothing-filter-for-not-equally-spaced-data)
|
|
28
|
+
に、不等間隔でも計算できるSavitzky–Golay法のPythonコードが載ってたので、それを微分にも使えるように修正しました
|
|
29
|
+
よろしければお試しください
|
|
30
|
+
```python
|
|
31
|
+
import numpy as np
|
|
32
|
+
|
|
33
|
+
def non_uniform_savgol_der(x, y, window, polynom, der):
|
|
34
|
+
"""
|
|
35
|
+
Applies a Savitzky-Golay filter to y with non-uniform spacing
|
|
36
|
+
as defined in x
|
|
37
|
+
|
|
38
|
+
This is based on https://dsp.stackexchange.com/questions/1676/savitzky-golay-smoothing-filter-for-not-equally-spaced-data
|
|
39
|
+
The borders are interpolated like scipy.signal.savgol_filter would do
|
|
40
|
+
|
|
41
|
+
Parameters
|
|
42
|
+
----------
|
|
43
|
+
x : array_like
|
|
44
|
+
List of floats representing the x values of the data
|
|
45
|
+
y : array_like
|
|
46
|
+
List of floats representing the y values. Must have same length
|
|
47
|
+
as x
|
|
48
|
+
window : int (odd)
|
|
49
|
+
Window length of datapoints. Must be odd and smaller than x
|
|
50
|
+
polynom : int
|
|
51
|
+
The order of polynom used. Must be smaller than the window size
|
|
52
|
+
der : int
|
|
53
|
+
The order of derivative. Must be positive and less than 3
|
|
54
|
+
|
|
55
|
+
Returns
|
|
56
|
+
-------
|
|
57
|
+
np.array of float
|
|
58
|
+
The smoothed y values
|
|
59
|
+
"""
|
|
60
|
+
if len(x) != len(y):
|
|
61
|
+
raise ValueError('"x" and "y" must be of the same size')
|
|
62
|
+
|
|
63
|
+
if len(x) < window:
|
|
64
|
+
raise ValueError('The data size must be larger than the window size')
|
|
65
|
+
|
|
66
|
+
if type(window) is not int:
|
|
67
|
+
raise TypeError('"window" must be an integer')
|
|
68
|
+
|
|
69
|
+
if window % 2 == 0:
|
|
70
|
+
raise ValueError('The "window" must be an odd integer')
|
|
71
|
+
|
|
72
|
+
if type(polynom) is not int:
|
|
73
|
+
raise TypeError('"polynom" must be an integer')
|
|
74
|
+
|
|
75
|
+
if polynom >= window:
|
|
76
|
+
raise ValueError('"polynom" must be less than "window"')
|
|
77
|
+
|
|
78
|
+
if polynom < 2:
|
|
79
|
+
raise ValueError('"polynom" must be larger than 1')
|
|
80
|
+
|
|
81
|
+
if type(der) is not int:
|
|
82
|
+
raise TypeError('"der" must be an integer')
|
|
83
|
+
|
|
84
|
+
if der < 0:
|
|
85
|
+
raise TypeError('"der" must be an positive integer')
|
|
86
|
+
|
|
87
|
+
if der > 2:
|
|
88
|
+
raise TypeError('"der" must be less than 3')
|
|
89
|
+
|
|
90
|
+
half_window = window // 2
|
|
91
|
+
polynom += 1
|
|
92
|
+
|
|
93
|
+
# Initialize variables
|
|
94
|
+
A = np.empty((window, polynom)) # Matrix
|
|
95
|
+
tA = np.empty((polynom, window)) # Transposed matrix
|
|
96
|
+
t = np.empty(window) # Local x variables
|
|
97
|
+
y_smoothed = np.full(len(y), np.nan)
|
|
98
|
+
|
|
99
|
+
# Start smoothing
|
|
100
|
+
for i in range(half_window, len(x) - half_window, 1):
|
|
101
|
+
# Center a window of x values on x[i]
|
|
102
|
+
for j in range(0, window, 1):
|
|
103
|
+
t[j] = x[i + j - half_window] - x[i]
|
|
104
|
+
|
|
105
|
+
# Create the initial matrix A and its transposed form tA
|
|
106
|
+
for j in range(0, window, 1):
|
|
107
|
+
r = 1.0
|
|
108
|
+
for k in range(0, polynom, 1):
|
|
109
|
+
A[j, k] = r
|
|
110
|
+
tA[k, j] = r
|
|
111
|
+
r *= t[j]
|
|
112
|
+
|
|
113
|
+
# Multiply the two matrices
|
|
114
|
+
tAA = np.matmul(tA, A)
|
|
115
|
+
|
|
116
|
+
# Invert the product of the matrices
|
|
117
|
+
tAA = np.linalg.inv(tAA)
|
|
118
|
+
|
|
119
|
+
# Calculate the pseudoinverse of the design matrix
|
|
120
|
+
coeffs = np.matmul(tAA, tA)
|
|
121
|
+
|
|
122
|
+
# Calculate c0 which is also the y value for y[i]
|
|
123
|
+
y_smoothed[i] = 0
|
|
124
|
+
for j in range(0, window, 1):
|
|
125
|
+
#y_smoothed[i] += coeffs[0, j] * y[i + j - half_window]
|
|
126
|
+
y_smoothed[i] += coeffs[der, j] * y[i + j - half_window]
|
|
127
|
+
|
|
128
|
+
# If at the end or beginning, store all coefficients for the polynom
|
|
129
|
+
if i == half_window:
|
|
130
|
+
first_coeffs = np.zeros(polynom)
|
|
131
|
+
for j in range(0, window, 1):
|
|
132
|
+
for k in range(polynom):
|
|
133
|
+
first_coeffs[k] += coeffs[k, j] * y[j]
|
|
134
|
+
elif i == len(x) - half_window - 1:
|
|
135
|
+
last_coeffs = np.zeros(polynom)
|
|
136
|
+
for j in range(0, window, 1):
|
|
137
|
+
for k in range(polynom):
|
|
138
|
+
last_coeffs[k] += coeffs[k, j] * y[len(y) - window + j]
|
|
139
|
+
|
|
140
|
+
# Interpolate the result at the left border
|
|
141
|
+
for i in range(0, half_window, 1):
|
|
142
|
+
y_smoothed[i] = 0
|
|
143
|
+
x_i = 1
|
|
144
|
+
#for j in range(0, polynom, 1):
|
|
145
|
+
for j in range(der, polynom, 1):
|
|
146
|
+
y_smoothed[i] += first_coeffs[j] * x_i
|
|
147
|
+
x_i *= x[i] - x[half_window]
|
|
148
|
+
|
|
149
|
+
# Interpolate the result at the right border
|
|
150
|
+
for i in range(len(x) - half_window, len(x), 1):
|
|
151
|
+
y_smoothed[i] = 0
|
|
152
|
+
x_i = 1
|
|
153
|
+
#for j in range(0, polynom, 1):
|
|
154
|
+
for j in range(der, polynom, 1):
|
|
155
|
+
y_smoothed[i] += last_coeffs[j] * x_i
|
|
156
|
+
x_i *= x[i] - x[-half_window - 1]
|
|
157
|
+
|
|
158
|
+
return y_smoothed
|
|
159
|
+
|
|
160
|
+
if __name__ == '__main__':
|
|
161
|
+
import numpy as np
|
|
162
|
+
import matplotlib.pyplot as plt
|
|
163
|
+
|
|
164
|
+
# テスト用ダミーデータ
|
|
165
|
+
x = np.arange(0, 100)
|
|
166
|
+
y = np.arange(10, 110) + np.random.randn(100) * 2
|
|
167
|
+
y11 = non_uniform_savgol_der(x, y, 11, 2, 0)
|
|
168
|
+
plt.plot(x, y, x, y11)
|
|
169
|
+
plt.grid()
|
|
170
|
+
plt.show()
|
|
171
|
+
|
|
172
|
+
# 微分
|
|
173
|
+
dy = np.gradient(y, x)
|
|
174
|
+
dy11 = non_uniform_savgol_der(x, y, 11, 2, 1)
|
|
175
|
+
plt.plot(x, dy, x, dy11)
|
|
176
|
+
plt.grid()
|
|
177
|
+
plt.show()
|
|
178
|
+
|
|
179
|
+
# 二階微分
|
|
180
|
+
ddy = np.gradient(dy, x)
|
|
181
|
+
ddy11 = non_uniform_savgol_der(x, y, 11, 2, 2)
|
|
182
|
+
plt.plot(x, ddy, x, ddy11)
|
|
183
|
+
plt.grid()
|
|
184
|
+
plt.show()
|
|
185
|
+
```
|
1
誤りの告知
answer
CHANGED
|
@@ -16,4 +16,9 @@
|
|
|
16
16
|
plt.grid()
|
|
17
17
|
plt.show()
|
|
18
18
|
```
|
|
19
|
-
点数が多いほど平滑化が強くなります
|
|
19
|
+
点数が多いほど平滑化が強くなります
|
|
20
|
+
|
|
21
|
+
【訂正】
|
|
22
|
+
これはxが等間隔じゃないと使えない手法です
|
|
23
|
+
データ見たら等間隔ではないので、使えませんね
|
|
24
|
+
失礼しました
|