バックナンバーはこちら。
https://www.simulationroom999.com/blog/stock-predict-matlabpython-backnumber/
はじめに
前回までで、MATLAB、Python(Numpy)によるFFT出力の周波数分布をローテーションをできることが確認できた状態。
今回から実際に超簡易バンドパスフィルタを実施していく。
登場人物
博識フクロウのフクさん
イラストACにて公開の「kino_k」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=iKciwKA9&area=1
エンジニア歴8年の太郎くん
イラストACにて公開の「しのみ」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=uCKphAW2&area=1
MATLABで超簡易バンドパスフィルタ
じゃ、ついにバンドパスフィルタをやるところだね。
そうそう。
前回までで、FFTと周波数分布のローテーションは終わってるから、
特定周波数の抽出とIFFTをやるだけだ。
まるで簡単なような言い方してるけど、あまり簡単だった試しが無いんだよなぁ。
そんなことは無いよ。
少なくともIFFTはすでにやってるものだし、
太郎くんから見て不明点は「特定周波数の抽出」の部分だけだと思うよ。
そうそう！そこそこ！
まぁこれもコードで見れば一行で終わるやつだから。
ということは、チョロいってことで良さそうか。
ちょっと通常の書き方とは違うかもいしれないけど。
(やっぱり罠がありそう・・・。)
MATALBコード
じゃ、作って来たコードだ。
N=1024;
L=pi;
x=linspace(-L,L,N)';
k=-N/2:N/2-1;
ft=sin(x)+sin(3*x)+sin(7*x);
% 変換前波形
subplot(4,1,1)
plot(x,ft)
title('f(t)');
grid();
% FFT後にローテーション
Fw=fft(ft);
Fw_tmp=circshift(Fw,N/2);
subplot(4,1,2);
mask = (-10<=k & k<=10);
plot(k(mask),abs(Fw_tmp(mask)));
title('F(\omega)');
grid();
% 特定周波数のみ抽出
Fw_Filter=Fw_tmp;
Hz = 3;
Low = Hz-0.1;
High = Hz+0.1;
Fw_Filter( k< -High | (-Low < k & k <= 0) | (0 <= k & k < Low) | High < k)=0;
subplot(4,1,3);
mask = (-10<=k & k<=10);
plot(k(mask),abs(Fw_Filter(mask)));
title('F(\omega) Filter');
grid();
% IFFT前にローテーション
Fx = (ifft(circshift(Fw_Filter,N/2)));
subplot(4,1,4);
plot(x,real(Fx));
title('f(x)');
grid();結果
そして、結果だ。
おー！
ちゃんと3Hzだけ取り出せてる！
というわけで、
無事、バンドパスフィルタができたというわけだ。
コードを見てると？
ちょっと疑問に思ったのだけど、
おそらく周波数分布で3[Hz]と-3[Hz]以外を0にしてる部分のコードが妙なことになってるように見える・・・。
Fw_Filter( k< -High | (-Low < k & k <= 0) | (0 <= k & k < Low) | High < k)=0;
あー、そこね。
論理インデックス検索って書き方だな。
なんだそれは？
やはりそこで引っかかったか。
じゃ、次回はそこは少し掘り下げよう。
(こうなることを知っててやりやがったな・・・。)
まとめ
まとめだよ。
- MATLABで超簡易バンドパスフィルタ実施。
- コード開示＆結果確認。
- 想定通り、3Hzだけ抽出で来た。
- 3[Hz]と-3[Hz]以外を0にしてる部分のコードが妙。
- 論理インデックス検索という手法を使っている。
バックナンバーはこちら。
コメント