回答編集履歴

2

修正

2017/11/17 10:14

投稿

mkgrei
mkgrei

スコア8560

test CHANGED
@@ -1,6 +1,8 @@
1
- ちゃんとバンドパスフィルタでしたよ?
1
+ デーが無いの、てきとーに生成ています。
2
2
 
3
- データ無いので、きとーに生成していますが
3
+ 実空間・周波数空間の変換ミスっていました
4
+
5
+ パラメータはダイレクトに入力して、付随して決まるものは全部計算によって出すと、ミスがすくなるかと。
4
6
 
5
7
  ![イメージ説明](880e1ea3c2fb43ff1f954a35a8f6f88d.png)
6
8
 

1

修正

2017/11/17 10:14

投稿

mkgrei
mkgrei

スコア8560

test CHANGED
@@ -2,7 +2,7 @@
2
2
 
3
3
  データが無いので、てきとーに生成していますが。
4
4
 
5
- ![イメージ説明](8ff11f034a473973d2169518a656cabf.png)
5
+ ![イメージ説明](880e1ea3c2fb43ff1f954a35a8f6f88d.png)
6
6
 
7
7
 
8
8
 
@@ -18,61 +18,69 @@
18
18
 
19
19
  from scipy import signal
20
20
 
21
- N = 10000
22
21
 
23
- fs = 250 # サンプリング周波数
24
22
 
25
- nyq = fs / 2 # ナイキスト周波数
23
+ N = 1000
26
24
 
27
- fc1 = 10 / fs # カットオフ周波数
25
+ fs = 250
28
26
 
29
- fc2 = 50 / fs # カットオフ周波数
27
+ nyq = fs / 2
30
28
 
31
- numtaps = 125 # フィルタ係数の数
29
+ fc1 = 4.
32
30
 
33
- dt = 0.04 # サンプリング間隔
31
+ fc2 = 30.
34
32
 
35
- t = np.arange(0, N*dt, dt) # 時間軸
33
+ numtaps = 255
36
34
 
35
+ dt = 1. / fs
36
+
37
+ t = dt*np.arange(N)
38
+
37
- freq = np.linspace(0, 1.0/dt, N) # 周波数軸
39
+ freq = np.fft.fftfreq(t.shape[-1])*fs
38
40
 
39
41
 
40
42
 
41
43
  # フィルタ係数
42
44
 
43
- b = signal.firwin(numtaps, [fc1, fc2], pass_zero=False)
45
+ b = signal.firwin(numtaps, [fc1, fc2], pass_zero=False, nyq=nyq)
44
46
 
45
47
  #b = signal.firwin(numtaps, [fc1, fc2])
46
48
 
47
49
 
48
50
 
49
- data_pca = np.sum((np.random.random()*np.sin(i*t/(N*dt)) for i in range(N)), axis=1)
51
+ x = np.sum((np.random.random()*np.sin(i*t/(N*dt)) for i in range(N)), axis=1)
52
+
53
+ x /= np.max(np.abs(x))
50
54
 
51
55
  # scipyによるFIRフィルタ
52
56
 
53
- data_pca = signal.lfilter(b, 1, data_pca)
57
+ kx = signal.lfilter(b, 1, x)
54
58
 
55
59
 
56
60
 
57
- F = np.fft.fft(data_pca)
61
+ F = np.fft.fft(kx)
62
+
63
+ print(F.shape)
58
64
 
59
65
  Amp = np.abs(F)
60
66
 
61
67
 
62
68
 
63
- skip=5
69
+ skip=1
64
70
 
65
71
  plt.subplot(211)
66
72
 
67
- plt.plot(t[::skip], F[::skip], marker='.', ls='none', label='t-F')
73
+ plt.plot(freq[::skip], F[::skip], marker='.', ls='none', label='t-F')
74
+
75
+ plt.xlim(0, 50)
68
76
 
69
77
  plt.legend()
70
-
71
- plt.xlim(0, 100)
72
78
 
73
79
  plt.subplot(212)
74
80
 
75
81
  plt.plot(freq[::skip], Amp[::skip], marker='.', ls='none', label='w-A')
82
+
83
+ plt.xlim(0, 50)
76
84
 
77
85
  plt.legend()
78
86