回答編集履歴
3
周波数特性グラフ追加
test
CHANGED
|
@@ -27,3 +27,17 @@
|
|
|
27
27
|
axs[1].set_xlim(0, 1500)
|
|
28
28
|
plt.show()
|
|
29
29
|
```
|
|
30
|
+
|
|
31
|
+
|
|
32
|
+
下記のようにして、適当にデータをでっち上げて確認してみました
|
|
33
|
+
|
|
34
|
+
```python
|
|
35
|
+
df = pd.read_csv(in_file, encoding='SHIFT-JIS')
|
|
36
|
+
```
|
|
37
|
+
↓ 変更
|
|
38
|
+
```python
|
|
39
|
+
t = np.arange(0, 5.5, 0.000005)
|
|
40
|
+
y = np.random.rand(len(t))*2-1
|
|
41
|
+
df = pd.DataFrame(np.vstack([t, y]).T)
|
|
42
|
+
```
|
|
43
|
+
|
2
短く
test
CHANGED
|
@@ -1,48 +1,3 @@
|
|
|
1
|
-
[Why it returns nan while using scipy.signal.filtfilt](https://stackoverflow.com/questions/50254625/why-it-returns-nan-while-using-scipy-signal-filtfilt)
|
|
2
|
-
と同じことが起きてるのだと思います
|
|
3
|
-
|
|
4
|
-
|
|
5
|
-
> def highpass(x, samplerate, fp, fs, gpass, gstop):
|
|
6
|
-
|
|
7
|
-
の関数内の、
|
|
8
|
-
|
|
9
|
-
> return y
|
|
10
|
-
|
|
11
|
-
のすぐ上に
|
|
12
|
-
|
|
13
|
-
```python
|
|
14
|
-
print(wp)
|
|
15
|
-
print(ws)
|
|
16
|
-
print(N)
|
|
17
|
-
```
|
|
18
|
-
|
|
19
|
-
を追加して実行して、結果表示を確認したら、フィルタの次数「N」が非常に大きくなってますよね
|
|
20
|
-
時間刻み「dt」で決まる「samplerate」に対する「fp」と「fs」の設定により、「wp」と「ws」の違いがわずかであり、非常に急峻なフィルタを設定してるためだと思います
|
|
21
|
-
|
|
22
|
-
|
|
23
|
-
「samplerate」を小さく(「dt」を大きく)するとか、「samplerate」に対する「fp」と「fs」の設定を見直す(「fp」と「fs」の差を大きくする)とかして、フィルタの次数「N」がもっと小さくなれば、「NaN」が出ずに計算できるようになると思います
|
|
24
|
-
たとえば、下記の変更のどちらかを行ってみてください
|
|
25
|
-
|
|
26
|
-
```python
|
|
27
|
-
fp_hp = 520
|
|
28
|
-
fs_hp = 500
|
|
29
|
-
```
|
|
30
|
-
↓ 変更
|
|
31
|
-
```python
|
|
32
|
-
fp_hp = 1000
|
|
33
|
-
fs_hp = 500
|
|
34
|
-
```
|
|
35
|
-
または
|
|
36
|
-
```python
|
|
37
|
-
fp_hp = 520
|
|
38
|
-
fs_hp = 250
|
|
39
|
-
```
|
|
40
|
-
|
|
41
|
-
数値は適当なので、上記変更でうまくいかなければ、「fp」と「fs」の差がもっと大きくなるように変えてみてください
|
|
42
|
-
|
|
43
|
-
|
|
44
|
-
そのあたりを変えたくないなら、下記の変更を試してみてください
|
|
45
|
-
|
|
46
1
|
```python
|
|
47
2
|
b, a = signal.butter(N, Wn, "high")
|
|
48
3
|
y = signal.filtfilt(b, a, x)
|
|
@@ -72,12 +27,3 @@
|
|
|
72
27
|
axs[1].set_xlim(0, 1500)
|
|
73
28
|
plt.show()
|
|
74
29
|
```
|
|
75
|
-
|
|
76
|
-
|
|
77
|
-
「df_filter」の周波数特性を見ると、回答の前半に書いた「fp」と「fs」の差を大きくする方法の結果よりも、後半に書いた下記に変更する方法の結果の方が、急峻なフィルタ特性が実現できることが分かります
|
|
78
|
-
(「fp」と「fs」を変えてないのだから当然ですが)
|
|
79
|
-
|
|
80
|
-
```python
|
|
81
|
-
sos = signal.butter(N, Wn, "high", output='sos')
|
|
82
|
-
y = signal.sosfiltfilt(sos, x)
|
|
83
|
-
```
|
1
追記
test
CHANGED
|
@@ -55,3 +55,29 @@
|
|
|
55
55
|
|
|
56
56
|
参考
|
|
57
57
|
[A.2P Pythonによるデジタルフィルタ](https://program4ptotat.livedoor.blog/archives/13572699.html)
|
|
58
|
+
|
|
59
|
+
|
|
60
|
+
【追記】
|
|
61
|
+
コードの最後に下記を追加したら、「df」と「df_filter」の周波数特性の比較ができます
|
|
62
|
+
|
|
63
|
+
```python
|
|
64
|
+
df_f = np.abs(np.fft.fft(df.values[:, 1], norm='ortho'))
|
|
65
|
+
df_filter_f = np.abs(np.fft.fft(df_filter.values[:, 1], norm='ortho'))
|
|
66
|
+
freq = np.fft.fftfreq(len(df_f), d=df.T.iloc[0,1])
|
|
67
|
+
|
|
68
|
+
fig, axs = plt.subplots(2, 1)
|
|
69
|
+
axs[0].plot(freq, df_f, '.')
|
|
70
|
+
axs[0].set_xlim(0, 1500)
|
|
71
|
+
axs[1].plot(freq, df_filter_f, '.')
|
|
72
|
+
axs[1].set_xlim(0, 1500)
|
|
73
|
+
plt.show()
|
|
74
|
+
```
|
|
75
|
+
|
|
76
|
+
|
|
77
|
+
「df_filter」の周波数特性を見ると、回答の前半に書いた「fp」と「fs」の差を大きくする方法の結果よりも、後半に書いた下記に変更する方法の結果の方が、急峻なフィルタ特性が実現できることが分かります
|
|
78
|
+
(「fp」と「fs」を変えてないのだから当然ですが)
|
|
79
|
+
|
|
80
|
+
```python
|
|
81
|
+
sos = signal.butter(N, Wn, "high", output='sos')
|
|
82
|
+
y = signal.sosfiltfilt(sos, x)
|
|
83
|
+
```
|