前回投稿後にTwitterで@rithmety様のご指導あり、気になっていたところが解決。ありがとうございます。これでFIRフィルタの計算ができると思ったら、いけません。定数PI（円周率）はどこにいるの？だいたいSIN関数はどこ？頭の固いRust素人（老人）は戸惑うことばかり。今回は被テスト波形を計算するところまで。
※『RustにいればRustに従え』関係記事 index はこちら
※動作確認は、Windows11のWSL2上にインストールしたUbuntu20.04LTS上のrustc 1.64.0 (a55dd71d5 2022-09-19) で行っています。
今回途上の目論見
前回、移動平均が計算できたので、その余勢を駆ってFIRローパスフィルタ（１００次予定）によって、遮断周波数を上回る波形をビシビシと抑制し、通過域の信号を取り出すべし、との目論見であります。別シリーズ投稿にてScilabを使用して既にやってみているのでお楽。それに次数は大きいけれどもやることは移動平均と似たようなもんだと。手ごろな例題か？
しかし、実験に際しては何か被テスト用に信号波形が必要っす。ノイズ波形に踏み込むとそれはそれで大変そうなので、まあフツーにSIN波を適当に２周波数分作って混ぜてやればいいんでないの、と思ったのです。が、RustでSIN求めるにはどうするの、円周率はどうなの、という体たらくです。なお、別シリーズでScilab使った以下の記事で波形つくってます。
手習ひデジタル信号処理(59) Scilab、入力信号定義用の関数を手作り
ScilabのSIN波形作成関数のコードが以下に。ぶっちゃけ今回は以下のコードのRust版ができればOKと。
// Generate Sine wave
// Fs: sampling frequency [Hz]
// nSample: number of samples
// fq: frequency [Hz]
// phdeg: phase [deg]
// mag: manitude
function [t,y]=genSin(Fs, nSample, fq, phdeg, mag)
t = ((1:nSample) - 1) / Fs;
y = mag*sin(2*%pi*(fq*t+(phdeg/360)))
endfunction
今回作成してみたRustプログラム
例によってアカラサマな mut 無、計算時にfor無のコードが以下です（最後の結果印字用にfor使ってますが。）
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)).collect::<Vec<f64>>(); // A3
let s784 = ts.clone().map(|t| gen_sine(t, 783.609, 0.0, 1.0)).collect::<Vec<f64>>(); // G5
let under_test = (0..n_samples).map(|i| s220[i] + s784[i]);
for item in under_test {
println!("{0:3.7}", item);
}
}
π（円周率）なんですが、定数を格納したモジュールがあり、そこにありました。
また単なる実験なので周波数など、FIRフィルタのボード線図みながら適当な数字を選べば十分なのでありますが、ことさらに、ソの音とかラの音とかを選んでみました。各音の周波数については以下のサイトを使わせていただきました。あざーす。
プチモンテ様 各音の周波数一覧 [Web Audio API]
実行結果
上記のRustプログラムを cargo buildした実行ファイルの出力をファイルにリダイレクト。それを表計算ソフトで読み込んでグラフにしてみたものが以下です（まだRustでGraph描く方法知らんし。）
まあ、計算できてるんでないかい？確かめてないけど、だいだいな感じ。いい加減な。