回答編集履歴

3

周波数特性グラフ追加

2022/10/13 02:15

投稿

jbpb0
jbpb0

スコア7651

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
+ ![「df」と「df_filter」の周波数特性の比較](https://ddjkaamml8q8x.cloudfront.net/questions/2022-10-13/8560c68c-0fb1-46c2-bb63-ca6f0a87401b.png)

2

短く

2022/10/13 02:04

投稿

jbpb0
jbpb0

スコア7651

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

追記

2022/10/13 01:58

投稿

jbpb0
jbpb0

スコア7651

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
+ ```