polynomial_oscillator

\[ P(x) = \sum_{n=0}^{\infty} a_n x^n = a_0 + a_1 x + a_2 x^2 + a_3 x^3 + \dots \]

以下は多項式オシレータのデモです。◎で示された制御点をドラッグすると波形が変わります。

多項式から波形を取り出す範囲を \([0, 1)\) とします。これで位相のスケーリングを行わずに済みます。詳細は「その他 -> 位相の計算」を参照してください。

波形の始点は 0 とします。これで多項式の a_0 が常に 0 となり、計算量が減ります。

多項式はユーザが指定した \(k\) 個の制御点を通過します。言い換えると、制御点によってユーザが波形を変更できるようにします。

この節では、指定した制御点を通るように \(a_n\) を求めます。

\(x\) が小さい制御点から順にインデックスを与えて \(x_0 < x_1 < \dots < x_{k+1}\) となるようにします。

インデックス \(0\) とインデックス \(k + 1\) の制御点は波形の両端です。インデックス \(0\) の制御点の値を \((x_0, y_0) = (0, 0)\) 、インデックス \(k + 1\) の制御点の値を \((x_{k+1}, y_{k+1}) = (1, 0)\) と固定します。 \(x_0\) と \(x_{k+1}\) の値は多項式から波形を取り出す範囲についての仕様に基づきます。 \(y_0\) は波形の始点を 0 とする仕様から 0 になります。 \(y_{k+1}\) はユーザが設定できるようにもできますが、簡略化のために 0 に固定しています。

ここまでの定義を使えば \(\mathbf{A} \mathbf{x} = \mathbf{b}\) の形にして \(a_n\) を求めることができます。

多項式 \(\hat{P}\) はベクトルを使って以下のように書き直すことができます。

\[ \hat{P}(x) = \sum_{n=1}^{k} a_n x^n = \begin{bmatrix} x & x^2 & x^3 & \cdots & x^k\end{bmatrix} \begin{bmatrix} a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_k\end{bmatrix} \]

ここで \(k\) 個の制御点の組 \((x_n, y_n)\) が分かっているので、以下の行列式を組み立てられます。

\[ \begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_k\end{bmatrix} = \begin{bmatrix} x_1 & x_1^2 & x_1^3 & \cdots & x_1^k \\ x_2 & x_2^2 & x_2^3 & & x_2^k \\ x_3 & x_3^2 & x_3^3 & & x_3^k \\ \vdots & & & \ddots & \vdots \\ x_k & x_k^2 & x_k^3 & \cdots & x_k^k \\ \end{bmatrix} \begin{bmatrix} a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_k\end{bmatrix} \]

上の式は \(\mathbf{b} = \mathbf{A} \mathbf{x}\) の形になっているので scipy.linalg.solve のようなソルバで解けます。

今回の応用では \(\mathbf{A}\) がそれほど大きくならないので、ソルバを自前で実装するときはガウスの消去法を使えば間に合います。また Wikipedia の LU decomposition の記事にもソルバのコードが掲載されていますが、対角成分が 0 となるときのピボッティングなどが省略されているので注意してください。以下は IntegerArpeggio で使用した、JavaScript で書かれたソルバのコードへのリンクです。

以下の条件を満たすゲイン \(g\) を求めて最大振幅を 1 に正規化します。

\[ \begin{aligned} & \text{Find} && g, \\ & \text{such that} && \max(|g P(x)|) = 1, \quad P(x) = \sum_{n=1}^{k} a_n x^n, \\ & \text{where} && x \in [0, 1). \end{aligned} \]

この問題は \(x \in [0, 1)\) の範囲における \(P(x)\) のピークを探して、最大のピークの絶対値の逆数を \(g\) とすれば解けます。

\(P(x)\) のピークは、以下のように \(P(x)\) を微分した式 \(P'(x)\) が 0 になる点と定義します。

\[ \mathtt{peaks} = \left\{ P(x) \mid \forall x \enspace \text{such that} \enspace P'(x) = 0 \right\}, \quad P'(x) = \frac{d P(x)}{dx}. \]

一般的には root finding アルゴリズムに \(P'(x)\) を渡せばピークが求められます。ただし root finding アルゴリズムは多数あり、選定も実装も手間です。そこで何とかならないかと試行錯誤していたところ、今回の応用には制御点があり、隣り合う制御点の間には最大でも 1 つのピークしか現れないことに気が付きました。以降ではこの観察に基づいて話を進めていきます。ただし、制御点の間のピークの数はあくまでも観察に基づく仮説であって、証明された事実ではないので例外があるかもしれません。

制御点の間のピークが 1 つ以下であれば、二分探索で安定して \(P'(x)\) の解を探すことができます。以下に \(g\) を求める手続きを書き下します。途中で出てくる \(\mathrm{sgn}\) は符号関数です。

以下は上の手続きを JavaScript で実装したコードへのリンクです。

関数 \(P'(x)\) はごく普通の多項式であり、特に精度が必要な応用でもないので、二分探索よりもニュートン法のほうが適しているかもしれません。

オシレータの位相の計算は様々なバリエーションがあるのでまとめておきます。コードは C++ です。

おすすめの方法です。 floor を使う利点は phase が負の値になっても [0, 1) の範囲に丸められることです。

phase += frequencyHz / sampleRateHz;
phase -= std::floor(phase);

std::fmod の利点は位相の範囲の上限を任意に設定できることです。欠点は割られる数 (1 つ目の引数) が負の値のときに出力も負の値となることです。したがって FM などによって負の周波数が出てくるときは対応が必要です。

phase = std::fmod(phase + frequencyHz / sampleRateHz, upperBound);
if (frequencyHz < 0) phase = -phase;

// あるいは

phase = std::fmod(phase + frequencyHz / sampleRateHz, upperBound);
phase *= std::copysign(float(1), frequencyHz);

また、 sin や cos の位相を回すときは fmod を使って upperBound = 2 * pi と設定するとよさそうに思えますが、私が試した範囲では floor を使ったほうが速いです。

constexpr auto twoPi = float(2) * std::numbers::pi_v<float>;

// fmod
phase = std::fmod(phase + twoPi * frequencyHz / sampleRateHz, twoPi);
output = std::sin(phase);

// floor
phase += frequencyHz / sampleRateHz;
phase -= std::floor(phase);
output = std::sin(twoPi * phase);

位相の値が上限を超えたときに分岐を設ける方法があります。

phase += std::clamp(frequencyHz / sampleRateHz, float(0), float(0.5));
if (phase >= float(1)) phase -= float(1);

// あるいは

phase += std::max(frequencyHz / sampleRateHz, float(0));
while (phase >= float(1)) phase -= float(1);

while を使う方法は frequencyHz の大きさに比例して計算量が増えるので、使わないでください。似たようなコードを見たことがあるので一応紹介しています。

std::remainder は引数が両方とも正の値のとき、出力が負の値となるので位相の回転には不適です。

符号なし整数を使えば丸めの処理はコンパイラが何とかしてくれます。 DSP のコードがすべて整数演算なら高速です。浮動小数点数を使っているときはキャストやスケーリングが入るので速度面での利点は微妙です。

constexpr auto scaler = float(std::numeric_limits<unsigned>::max());
unsigned frequencyInt = unsigned(scaler * frequencyHz / sampleRateHz);
phase += frequencyInt;

// 使用例。
output = std::sin(twoPi * float(phase) / scaler);

標準ライブラリの <cinttypes> に様々な大きさの整数が定義されています。

波形をルックアップテーブルから読み出すときは符号なし整数とビットマスクを組み合わせる実装があります。ただし、ルックアップテーブルの長さは 2^n に制限されます。

constexpr auto scaler = unsigned(0xfffff); // 任意の 2^n - 1 の定数。
unsigned frequencyInt = unsigned(float(scaler) * frequencyHz / sampleRateHz);
phase += frequencyInt;
phase &= scaler;

// 使用例。
output = wavetable[phase];

以下は前述の floor による位相の回転を固定小数点数で実装した例です。

// 定義。
constexpr auto fractionBits = int32_t(16);
constexpr auto fractionMult = int32_t(1) << fractionBits;
constexpr auto fractionMask = fractionMult - int32_t(1);
int32_t phase = 0;
int32_t sampleRateFixed = int32_t(sampleRateHz) << fractionBits;
int32_t frequencyFixed = int32_t(frequencyHz) << fractionBits;

auto div = [](int32_t x, int32_t y) -> int32_t {
  return int32_t((int64_t(x) * fractionMult) / y);
};

// 位相の回転。
phase += div(frequencyFixed, sampleRateFixed);
phase &= fractionMask; // phase -= floor(phase) と同じ意味。

上記のコードの固定小数点数のフォーマットは符号 1 ビット、整数部 15 ビット、小数部 16 ビットで、 Q15.16 あるいは Q16.16 と表記されます。 Qm.n という表記は Q フォーマットと呼ばれ、小数点より上の桁数が m ビット、小数点より下の桁数が n ビットの固定小数点数という意味です。

Q フォーマットは Texas Instruments の資料 (Appendix A.2, p.A-3) と ARM の資料 (4.7.9 Q-format, p.4-24) に掲載されていますが、 m に符号ビットを含むかどうかで違いがあります。 Texas Instruments では符号ビットを含まない、 ARM では符号ビットを含む、としています。

固定小数点数を自前で実装すると数学関数の用意が大変なので、特に理由がなければ fpm のようなライブラリを使うことを推奨します。