今製作中のエンベロープジェネレータは、FPU のないSTM32G0 を使って計算によりエンベロープを生成しています。40 kHz サンプリングで2ボイス分の生成を行うので演算速度への要求はかなり厳しく固定小数点演算が必須です。Embassy での固定小数点演算について調べたことを簡単にまとめてみました。
どうやって実装するか？
基本の ADSR 曲線を生成する場合、以前検討した通りコンデンサへの充放電をシミュレーションする方法では、基本の計算式は以下のように非常に単純になります
loop {
current_value += (target - current_value) * k;
}ここで、 target は充放電で最終的に行きたいレベル、current_value は現在の値、k は充放電の速さをあらわす定数です。充放電は収束しないといけないので k の範囲は必ず 0 <= k < 1 でなくてはいけません。
これを固定小数点で演算してゆくわけですが、Embassy で固定小数点演算を実装する方法は大まかに以下の二択です。
- fixed クレイトを使う
- 整数とビットシフトを使って自分で実装する (Rust で書く)
アセンブリを使って実装する方法もあるかなと思うのですがそちらはまだ勉強していないので保留です。詳細は後述しますが、固定小数点演算の精度は以下の二つを検討しました。
- Q32.32 （符号付き 64 bit 整数、32 bit が小数部分）
- UQ0.32 （符号なし 32 bit 整数、32 bit が小数部分) – 負数は処理の場合分けで避ける
その性能は？
いつもそうですけど詳細を書くと大変煩雑になるので先に結果から。各条件で同じエンベロープ計算プログラムを走らせた場合の実行時間が以下です。実際に使ったプログラムは github の issue に掲載しています。計算回数は 40 万回で、40 kHz サンプリングのエンベロープを 10 秒生成する分量です。
| 実装 | 負数回避チェックなし | 負数回避チェックあり |
|---|---|---|
| I32F32 (fixed) | 4037 ms | n/a |
| i64 Q32.32 (自前) | 3836 ms | n/a |
| U0F32 (fixed) | 589 ms | 657 ms |
| u32 UQ0.32 (自前) | 588 ms | 635 ms |
結果を見て言えるのは
- Q32.32 の演算は全般非常に遅い
- fixed クレイトと自前実装の実行速度はほぼ同じ。単純な整数演算の自前実装のほうがほんの少し速い
Q32.32 の遅さは致命的で、10 秒の信号を 4 秒かけて計算するのを 2 ボイス分やっているのでおそらく間に合っていません。実は Q32.32 の実装に入れ替えたところで「あれ？音が変になったかも」と思ったのでこの検討を始めたのです。これまた後述しますが、この遅さは Q32.32 で必要な 128bit 演算からきていると思われます。
結果を見てのさしあたっての実装の選択
結果を見て、またその他の観察から実装の選択に関していえることは
- Q32.32 の固定少数点演算はやってはいけない。全力で避ける (fixed を使う使わないは無関係)
- fixed クレイトを使っても自分で実装しても実行速度には大差ない
- プログラムコードの読みやすさの観点では fixed の方がやや勝る
- ビルドしたイメージは fixed を使わない方が小さい。10 KB ぐらい差が出る
ということで、ただでさえイメージが大きくなって困る Embassy で今回のプロジェクトでは fixed 導入は良いことよりも悪いことの方が多いと判断し外しました。
検討の詳細
詳細を書くとどうしても煩雑になってしまうので読んでくれる人いるかなあと不安に思いつつ、結構自分で後で読んだりもするので書きます。検討内容はざっと github エンベロープジェネレータプロジェクトの issue に書いたんですが、検討内容の背景などはそちらでもあまり書かなかったのでこちらに書きます。
固定小数点の精度の決定
固定小数点演算をする場合、整数部と小数部にどれぐらいのビット数を割くか決めることが精度良い計算をするため重要になります。
製作中のエンベロープジェネレータの出力は 12 bit で、それに見合う精度の計算をするには 32 bit あれば十分なようです。これは実験を繰り返して試行錯誤で決めました。また、内部ではプログラムをわかりやすくするために、値の範囲を 0.0 から 1.0 の間に限定しています。その範囲のうち、0.5 の値が DAC の最大出力となるように設定しています。こうすると、固定小数点の表現はほとんどの場合 UQ0.32 で事足ります。一つ問題なのは、計算の途中で負の値が出てくる場合が時々あること。これを条件分岐をせずに処理するには整数部と小数部が同じ長さであると手っ取り早く、そういうわけで符号付きの Q32.32 の精度が出てきました。
実装は？
Q32.32 の符号付き演算はの実装は、fixed クレイトを使うと以下のようになります。
let k = I32F32::from_num(0.00001f32);
let mut current_value = I32F32::from_bits(0x7fffffff);
let mut target = I32F32::from_num(0.2f32);
let n = 40_000 * 10;
for _ in 0..n {
let delta = (target - current_value) * k;
current_value += delta;
}
let output = current_value.to_bits();ライブラリは非常に巧妙に実装してあって、固定小数点演算しているようには外からは見えません。一方、整数を使って自分で実装する場合以下の通りです。
//// Multiplies two Q32.32 numbers. Returns Q32.32 result.
#[inline(always)]
fn mul_q32_32(a: i64, b: i64) -> i64 {
((a as i128 * b as i128) >> 32) as i64
}
{
let k: i64 = k.to_bits(); // derived from the fixed I32F32 value above
let mut current_value: i64 = 0x7fffffff;
let target: i64 = target.to_bits();
for _ in 0..n {
let delta = mul_q32_32(target - current, k);
current_value += delta;
}
let output = current_value;
}mul_q32_32() 関数で固定小数点の掛け算をしますが、整数部と小数部が同じ長さだと負数の処理も単純なビットシフトで済んでしまいます。ところが、プログラムの字面はすっきりしていますが中身はどうもそうでもないようです。問題は 128 bit 整数を使っていることのようで、レジスタのやりくりが大変ならしく、この関数の実行には大変な時間がかかります。fixed クレイトも同様のしょりをしているらしくほぼ同じ時間がかかります。
もう一つ検討した固定小数点の精度が符号なし小数点部のみの UQ0.32 です。整数部がないというのはちょっと極端に見えますが、前述のように、エンベロープジェネレータの内部値はこれで問題ないように設定されています。一つ問題なのが符号なしなので負の値の処理ができないこと。演算時には注意深く負の値を避ける必要があります。
UQ0.32 の実装は以下の通り。fixed を使ったものと自前実装したものです。
#[inline(always)]
fn mul_uq0_32(a: u32, b: u32) -> u32 {
((a as u64 * b as u64) >> 32) as u32
}
{
// using crate fixed
let k = U0F32::from_bits(k.to_bits() as u32);
let mut current_value = U0F32::from_bits(0x7fffffff);
let target = U0F32::from_bits(target.to_bits() as u32);
for _ in 0..n {
if current_value > target {
let delta = (current_value - target) * k;
current_value -= delta;
} else {
let delta = (target - current_value) * decay_ratio;
current_value += delta;
}
}
let output = current_value.to_bits();
}
{
let k: u32 = k.to_bits() as u32;
let mut current_value: u32 = 0x7fffffff;
let target: u32 = target.to_bits() as u32;
for _ in 0..n {
if current_value > target {
let delta = mul_uq0_32(current_value - target, decay_ratio);
current_value -= delta;
} else {
let delta = mul_uq0_32(target - current_value, decay_ratio);
current_value += delta;
}
}
let output = current_value;
}上記のプログラムでは if 分岐するまでもなく current_value > target ですが、実際にはノブを回すことによって target が動いたりするので current_value < target となる場合も起こります。しかしそういうわけなので、条件分岐は大きく片側に偏り分岐予測が効果的に効くようで、条件分岐による実行時間のペナルティはあまり大きくありません。
