前回の投稿後に再び@rithmety様のご指導ありコードを改良することができました。ありがとうございます。おかげで今回はFIRフィルタの計算に取り組めます。100次のローパスフィルタとな。結構強力っす。でもフィルタの実体はほとんど１行。Rustも強力にして簡潔。例によってアカラサマなmutなし。表示以外のforなし。
※『RustにいればRustに従え』関係記事 index はこちら
※動作確認は、Windows11のWSL2上にインストールしたUbuntu20.04LTS上のrustc 1.64.0 (a55dd71d5 2022-09-19) で行っています。
今回の計算例
今回のコードの前半は、前回作成の被テスト波形の生成部分です。@rithmety様のご指導により、ZIP使った方法に改めたもの。ただし、その後のフィルタ処理で簡潔だった@rithmety様コードに若干尾ひれがついてます。
被テスト波形は正弦波220Hzと正弦波783.689Hzを足し合わせたもの。そしてそこに
-
- カットオフ周波数 500Hz、ローパス
- 100次
のFIRフィルタを適用してみると。さすれば哀れ783.689Hzは抑制され、220Hzの低い信号だけが残る筈。
実験に使用したRustのソースコード
いまだRustのファイル入出力とか画像出力とかおぼつかないので、１００次のフィルタの係数をmain内でリテラルで定義しているという横暴さ、出力は標準出力へ数値表示という手抜きのコードであります。まあ、動けばよいのだ、と。
use std::f64::consts;
fn gen_sine(t: f64, fq: f64, phdeg: f64, mag: f64) -> f64 {
mag * (2.0 * consts::PI * (fq * t + (phdeg / 360.0))).sin()
}
fn main() {
let fs: f64 = 1.0e4; // Sampling Frequency [Hz]
let n_samples: usize = 4096; // Number of Samples
let ts = (0..n_samples).map(|i| i as f64 / fs);
let s220 = ts.clone().map(|t| gen_sine(t, 220.0, 0.0, 1.0));
let s784 = ts.clone().map(|t| gen_sine(t, 783.609, 0.0, 1.0));
let under_test = s220.zip(s784).map(|(s220, s784)| s220 + s784).collect::<Vec<f64>>();
// FIR LOW-PASS filter, Fs=10kHz, n=100, cut off=500Hz
let filtlen: usize = 100;
let firfilt: [f64; 100] = [
0.001006, 0.0029796, 0.0047385, 0.0060993, 0.0069097, 0.007065, 0.0065199, 0.005296, 0.0034822, 0.0012295,
-0.0012606,-0.0037535,-0.0060021,-0.0077703,-0.0088561,-0.0091128,-0.0084662,-0.0069255,-0.0045876,-0.0016326,
0.001688, 0.0050705, 0.0081847, 0.0107025, 0.0123291, 0.0128323, 0.0120688, 0.0100035, 0.0067214, 0.002429,
-0.0025536,-0.0078113,-0.0128617,-0.0171889,-0.0202833,-0.0216821,-0.0210086,-0.0180063,-0.0125661,-0.0047423,
0.0052415, 0.0170011, 0.0300105, 0.0436333, 0.057162, 0.0698647, 0.0810332, 0.0900316, 0.0963398, 0.0995893,
0.0995893, 0.0963398, 0.0900316, 0.0810332, 0.0698647, 0.057162, 0.0436333, 0.0300105, 0.0170011, 0.0052415,
-0.0047423,-0.0125661,-0.0180063,-0.0210086,-0.0216821,-0.0202833,-0.0171889,-0.0128617,-0.0078113,-0.0025536,
0.002429, 0.0067214, 0.0100035, 0.0120688, 0.0128323, 0.0123291, 0.0107025, 0.0081847, 0.0050705, 0.001688,
-0.0016326,-0.0045876,-0.0069255,-0.0084662,-0.0091128,-0.0088561,-0.0077703,-0.0060021,-0.0037535,-0.0012606,
0.0012295, 0.0034822, 0.005296, 0.0065199, 0.007065, 0.0069097, 0.0060993, 0.0047385, 0.0029796, 0.001006];
let result = (0..n_samples - filtlen + 1).map(|i| {
(0..filtlen).fold(0.0, |acc, x| under_test.get(i + x).map_or(acc, |fv| acc + fv*firfilt.get(x).unwrap()))
});
for item in result {
println!("{0:3.7}", item);
}
println!();
}
実験結果
under_testに格納されている被テスト波形、前回グラフ化したものを再掲します。こんな感じ（まだRustでグラフ描けないので表計算ソフトっす。）
２つの周波数が混じっているので、結構ガタガタな波形です。
そして以下がFIRフィルタを適用後の波形です。標準出力をファイルにリダイレクトして表計算ソフトにに取り込んでグラフにしてます。
ローパスなんで当然なんですけど、低い周波数だけ綺麗に取り出せてますなあ。100次ともなると切れ味鋭いっす。上のグラフから周波数を「だいたい」で計算すると222Hzとな。まあ、予定通りでないかい。
まあね、ホントはFFTかけてフィルタ前後の周波数を調べたいところなんだけれども、まだRustでFFTかけられんし。その前にファイル入出力とか当然なあたりやっておけよ、自分。