BLIT法メモ

https://ccrma.stanford.edu/~stilti/papers/blit.pdf

この論文では、アナログシンセサイザーの波形（パルストレインや鋸歯波）をエイリアスのないデジタル波形として合成する技術を紹介しています。方法としては、バンドリミットされた基本波形の総和やテーブルルックアップ技術が使われます。バンドリミットされたパルス波形や三角波形は、位相がずれた2つのバンドリミットされたインパルストレインの差を積分することで作り出されます。これらのインパルストレインは、ウィンドウドシンク関数の重ね合わせによって生成されます。さらに、一般的なバンドリミット波形の合成方法も紹介されています。これらの方法は、音質、計算効率、そして制御のしやすさという観点から評価されています。

1 イントロダクション

パルストレインや鋸歯波のように波形に不連続性があるアナログ信号や、三角波のように波形の傾きに不連続性がある信号は、対応する離散時間信号を得るためにサンプリング前にサンプリングレートの半分未満にバンドリミットされる必要があります。これらの波形をデジタルで生成する単純な方法は、不連続時間を最も近い利用可能なサンプリングインスタントに丸める必要があるため、エイリアスが含まれます。ここで主に扱われる信号はインパルストレイン、矩形パルス、および鋸歯波です。後者の2つの信号は、インパルストレインから積分によって導き出すことができるため、インパルストレインのアルゴリズムだけが詳細に説明されています。

バンドリミットの必要性:

• 不連続性のある波形: パルストレイン（パルス列）や鋸歯波のような波形には不連続性があり、三角波は波形の傾きに不連続性があります。
• バンドリミット: これらの波形をサンプリングしてデジタル信号に変換する際には、サンプリングレートの半分未満に周波数を制限（バンドリミット）する必要があります。これは、ナイキスト周波数（サンプリングレートの半分）を超える周波数成分がエイリアシングを引き起こさないようにするためです。

2. 単純な生成方法の問題点:
• エイリアス: 単純に波形をデジタル生成すると、不連続点を最も近いサンプリングタイミングに丸める必要があるため、エイリアスが発生します。エイリアスは、本来存在しない高周波成分が低周波成分として誤って再現される現象です。

3. 主に扱う信号:
• 対象の信号: この論文では主に、インパルストレイン、矩形パルス、鋸歯波の3つの信号について議論しています。
• 詳細なアルゴリズムの説明: 矩形パルスや鋸歯波はインパルストレインから積分によって導き出せるため、詳細に説明されるのはインパルストレインの生成アルゴリズムです。

2.なぜ単純な離散時間パルストレインはエイリアスを引き起こすのか

「明らかな」方法として、インパルストレインの離散時間バージョンを生成するには、それをユニットサンプルパルストレインで近似することが考えられます。ユニットサンプルパルス $${(\delta(n))}$$ は以下のように定義されます。

$$
\delta(n) \triangleq \begin{cases}
1 & n = 0 \
0 & |n| = 1, 2, 3, \ldots
\end{cases}
$$

ユニットサンプルパルスは整数 (n) に対してのみ定義されるため、問題が生じます。例えば、所望のインパルストレインの周波数 (f_1 = 1/T_1) とした場合、サンプルの周期は (P = T_1/T_s = F_s/f_1) となります。ここで (F_s = 1/T_s) はサンプリングレートです。この周期 (P) はほとんどの場合整数ではありません。ピッチの認識が非常に正確であるため、(M) を最も近い整数に丸めることは、周波数が非常に低い場合（約100 Hz以下）を除いて機能しません。そのため、ピッチを正確にするためには、インパルスの到達時間を非常に正確に計算し、各到達時間を最も近いサンプルインスタントに丸める必要があります。このプロセスはピッチ周期ジッタとしてモデル化でき、信号にノイズを加えます。また、これは幅が1サンプルの矩形パルスの周期列を一様にサンプリングすることとしてモデル化できます（図1）。

$$
p(t) \triangleq \begin{cases}
1 & |t| \leq T_s/2 \
0 & |t| > T_s/2
\end{cases}
$$

このパルスのフーリエ変換はsinc関数です。

$$
P(\omega) \triangleq \int_{-\infty}^{\infty} p(t) e^{-j\omega t} dt = T_s \text{sinc}(f T_s)
$$

$$
\text{sinc}(x) \triangleq \frac{\sin(\pi x)}{\pi x}
$$

矩形パルスの周期列は以下のように構成されます。

$$
x(t) \triangleq \sum_{l=-\infty}^{\infty} p(t + lT_1)
$$

そのフーリエ変換は、シフト定理およびデルタ関数の性質により、

$$
X(\omega) \triangleq \sum_{l=-\infty}^{\infty} e^{j\omega l T_1} P(\omega) \propto P(\omega) \sum_{l=-\infty}^{\infty} \delta(\omega - l\omega_1)
$$

矩形パルス列のフーリエ変換は、サンプリングレートの倍数でゼロとなるシンク関数の重み付けされた無限の高調波スペクトルを持ちます。サンプリング (x(t)) は、次のようにエイリアスを引き起こします。

$$
y(n) \triangleq x(nT_s) \propto \sum_{k=-\infty}^{\infty} X(\omega + k2\pi F_s) = \sum_{k=-\infty}^{\infty} P(\omega + k2\pi F_s) \sum_{l=-\infty}^{\infty} \delta(\omega + k2\pi F_s - l\omega_1)
$$

望ましいパルストレインスペクトルは (k = 0) 項のみです。各非ゼロの (k) 項は、全周波数帯域にわたる高調波のエイリアスの列を引き起こします。エイリアスの量は非常に大きく、聞こえます。唯一エイリアスが発生しない周波数はDCです。シンクスペクトルの包絡線がサンプリングレートの倍数でゼロを通過するためです。この理由から、サンプリングレートに対して非常に低い周波数ではエイリアスが減少します。

1. ユニットサンプルパルスの近似:
   • インパルストレインを生成するための「明らかな」方法として、ユニットサンプルパルス (\delta(n)) を使用します。しかし、これは整数のサンプル点にのみ適用されるため、非整数の周期 (P) を持つ信号には適用できません。
   2. 不正確な周波数の問題:
   • 所望の周波数 (f_1) に対する周期 P が整数でない場合、インパルスの到達時間を最も近いサンプルインスタントに丸める必要があります。これがピッチ周期ジッタとしてノイズを引き起こし、エイリアスの原因となります。
   3. 矩形パルスのフーリエ変換:
   • 矩形パルスのフーリエ変換はシンク関数で表されます。これにより、矩形パルスの周期列のスペクトルが無限の高調波を持ち、サンプリングによりエイリアスが発生します。
   4. エイリアスの可聴性:
   • エイリアスの量は非常に大きく、聞こえます。唯一エイリアスが発生しないのはDC（直流成分）で、これはシンク関数の包絡線がサンプリングレートの倍数でゼロになるためです。

3.3 バンドリミット補間（Bandlimited Interpolation）

一般的な離散時間信号 ( x(nT_s) ) がその最高周波数の2倍で一様にサンプリングされた場合、正確なバンドリミット補間（シンク補間と呼ばれる）を以下のように行います：

$$
x(t) = \sum_{n=-\infty}^{\infty} x(nT_s) h_s(t - nT_s)
$$

ここで、

$$
h_s(t) = \text{sinc}(F_s t) = \frac{\sin(\pi F_s t)}{\pi F_s t}
$$

これは、理想的なローパスフィルタのインパルス応答です。新しいサンプリングレート ( F_s’ = 1/T_s’ ) で ( x(t) ) を再サンプリングするには、整数倍の ( T_s’ ) で式 (2) を評価する必要があります。

式 (2) の総和は、理想的なローパスフィルタのインパルス応答が無限に伸びるため、実際には実装できません。インパルス応答を有限にするためには、ウィンドウ処理を行う必要があります。これがデジタルフィルタ設計のウィンドウ法の基礎です（RabinerとGold 1975を参照）。このフィルタインパルス応答は非常に長いため、理想的なバンドリミット補間として頻繁に使用されます。

3.4 正確なウェーブテーブル補間（Exact Wavetable Interpolation）

補間された周期的ウェーブテーブル合成で最高の結果を得るために、シンク補間が理想的です。信号が周期的であるため、式 (2) は一周期に収束し、有限和として実装できます（Schanze 1995）：

$$
x(t) = \frac{\sin(\pi t)}{2N} \sum_{n=-L}^{M-1} x(n) (\pm 1)^n \left[ \cot\left(\frac{\pi (t-n)}{2N}\right) + (\pm 1)^{N+1} \tan\left(\frac{\pi (t-n)}{2N}\right) \right]
$$

ここで、サンプリング間隔は ( T = 1 ) に正規化されており、周期は ( N = L + M ) サンプルです。これは、 ( x(t) ) の一周期にわたる時間エイリアスとして解釈できます。したがって、次のように示されます：

$$
x(t) = \sum_{\text{period}} x(nT_s) \text{Sinc}_N[F_s(t \pm nT_s)] C_N(n)
$$

ここで、

$$
C_N(n) = \begin{cases}
1, & \text{N奇数の場合} \
\cos(\pi F_s(t \pm n)/N), & \text{N偶数の場合}
\end{cases}
$$

3.5 DSF合成（DSF Synthesis）

DSF（離散和公式）は、バンドリミットされた周期信号を生成するために使用されます（Moorer 1975、Dodge and Jerse 1985、Roads 1996）。次の等式に基づいています：

$$
\sum_{k=0}^{N-1} a^k \sin(\theta + k\beta) = \frac{1}{1 - 2a \cos(\beta) + a^2} \left[ \sin(\theta) - a \sin(\theta - \beta) - a^N \sin(\theta + N\beta) + a^{N+1} \sin(\theta + (N-1)\beta) \right]
$$

この閉形式表現は、左辺に等式 $${( 2j \sin(x) = e^{jx} - e^{-jx} ) }$$を適用し、幾何級数の閉形式表現を適用することで導出されます：

$$
\sum_{k=0}^{N-1} z^k = \frac{1 - z^N}{1 - z}
$$

この式を用いて、広範なバンドリミット波形を生成します。バンドリミットは、生成される最高周波数が ( F_s/2 ) 未満になるように N を制御することで達成されます。例えば、 \theta = 0 および \beta = f_1 T_s N とすると、バンドリミット周波数に簡単に一般化されます。図4に示すように、DSFテーブルを用いて合成することができます。

3.6 バンドリミット・インパルス・トレイン（BLIT）合成

前述の技術のいくつかはBLITを生成するために使用できますが、他の波形も生成できます。ここでは、他の波形の拡張を考慮せずに、BLITを生成するために特別に設計された技術について説明します。特定の波形に限定することで、これらの方法は他の方法よりも効率的になります。また、生成する信号がインパルストレインであるため、数学的に他の波形よりもはるかに単純です。

3.7 シンクの和とSincM

サンプリング前の標準的な操作はアンチエイリアスフィルターを適用することです。理想的なアンチエイリアスフィルターの連続時間インパルス応答は、ゼロ交差間隔が1サンプルのシンク関数です。

$$
h_s(t) \triangleq \text{sinc}(F_s t) \triangleq \frac{\sin(\pi F_s t)}{\pi F_s t}
$$

理想的な単位振幅のインパルストレインの周期が (T_1) 秒であるとすると、

$$
x(t) = \sum_{l=-\infty}^{\infty} \delta(t + l T_1)
$$

この信号にアンチエイリアスフィルター (h_s) を適用すると

$$
x_f(t) \triangleq (x \ast h_s)(t) = \sum_{l=-\infty}^{\infty} h_s(t + l T_1) = \sum_{l=-\infty}^{\infty} \text{sinc}\left(\frac{t}{T_s} + l P\right)
$$

ここで (P = T_1 / T_s) はサンプル周期です（整数ではありません）。 (x_f) は周波数帯域 ((-F_s/2, F_s/2)) にバンドリミットされているため、次のようにサンプリングできます。

$$
y(n) \triangleq x_f(n T_s) = \sum_{l=-\infty}^{\infty} \text{sinc}\left(n + l P\right)
$$

この式は時間エイリアスを (P) サンプルの間隔でシンク関数に適用したものとして解釈できます。時間エイリアス化されたシンクは次のようになります。

$$
y(n) = \left(\frac{M}{P}\right) \text{Sinc}_M\left(\frac{M}{P} n\right)
$$

ここで

$$
\text{Sinc}_M(x) \triangleq \frac{\sin(\pi x)}{M \sin(\pi x / M)}
$$

この関数はサンプル化されたバンドリミット・インパルス・トレイン（BLIT）の閉形式表現を提供し、DSFに似た方法で直接合成に使用できます。 (P) はサンプル周期、(M) は高調波の数であり、常に奇数です。なぜなら、インパルストレインにはDC成分が1つ、他に偶数個の非ゼロ高調波が存在するためです（サンプリングレートの半分に正確に一致する高調波は許可しません）。 (M / P) は常に1に近く、(P) が奇数の場合、 (P = M) であり、 (y(n)) は単に (\text{Sinc}_M(n)) となります。 (P) が (M) から離れると、式(3) は時間スケーリングと振幅スケーリングを実装します。

3.8 窓付きシンクの和（BLIT-SWS）

デジタルインパルストレインを合成するより効率的な方法は、一般的なバンドリミット補間のための窓付きシンク法に基づいているかもしれません。この技術は概念的には、インパルストレインのバンドリミット周期的ウェーブテーブル合成と同等です。バンドリミット補間を使用して、サンプルレートを整数に割り切れるピッチから所望のピッチに変換します。このレート変換により、各ユニットサンプルパルス (\delta(n)) がウィンドウ付きシンク関数 (w(t)h_s(t)) に置き換えられます。

高調波（エイリアス）の減衰速度とゼロ交差数

窓付けにより高調波の減衰速度が制約されるため、エイリアスは避けられません。ウィンドウの選択によってこれを制御することができます。さらに、オーバーサンプリング係数を持つことは、可聴範囲の上限とサンプリングレートの半分の間に十分なガードバンドを確保するために役立ちます。これにより、必要なウィンドウ長が短くなります。

高調波数とシンクの重なり

周波数が非常に高い場合、バンドリミットされた高調波数は非常に少なくなります。実際、最高オクターブでは基本周波数のみがバンド内にあります。このため、多くの周波数範囲で、BLIT（または他の高調波波形）をサイン波の単純な総和によって生成する方が効率的である可能性が高いです。低周波数では、BLIT-SWS法がより効率的であるため、再度それを使用することができます。

最適なトレードオフのポイントは、アルゴリズムを実装するシステムによって異なります。極端な例では、FFT-1システムでは、トレードオフ周波数が (F_s/(2N)) まで下がります。これはシステムがサイン波の総和を非常に効率的に実装するためです。逆に、サイン波生成がかなり高価なシステムでは、トレードオフ周波数が F_s/8 程度になることがあります。

実際には、BLIT-SWS法は他のより正確なバンドリミットBLITアルゴリズムに対しても利点があります。周波数がスイープしている場合（ヴィブラートやポルタメントなど）、信号の高調波が消えるか現れます（周波数スイープの方向によります）。正確にバンドリミットされたシステムでは、最高調波が1サンプルの期間で完全な振幅とゼロ振幅の間を移動します。これは特に低サンプリングレートでは聞こえるトランジェントを引き起こします。BLIT-SWS法はウィンドウ処理のため、これらのトランジェントを自然に回避できます。

4.エイリアスフリーの鋸波と方形波の生成

次に扱う主要なアナログ波形のクラスは、方形波と鋸波です。これらをBLIT（バンドリミット・インパルス・トレイン）を使用して積分することで、エイリアスフリーな形で生成する方法を紹介します。

4.1 BLITの逐次積分

鋸波

鋸波関数は以下のように生成できます：

$$
\text{Saw}_{\text{CTS}}(t) = \int_0^t \text{CIT}(\tau) d\tau
$$

ここで、$${\text{CIT}(\tau) }$$は連続時間のインパルストレイン、$${\text{C}_1 = \int_0^T \text{CIT}(\tau) d\tau }$$はインパルストレインのDC成分です。これにより、連続時間の信号を直接、ディスクリート時間に変換できます（インパルス不変変換を介して）。

$$
\text{Saw}(n) = \sum_{k=0}^n \text{BLIT}(k) \leftrightarrow \text{C}_2
$$

$$
\leftrightarrow \text{Saw}(z) = \frac{z}{z-1}(\text{BLIT}(z) \leftrightarrow Z(\text{C}_2))
$$

これは単純な総和と一極のデジタルフィルターを用いた場合に自明に実装できます。\text{C}_2 はBLITの平均値であり、信号を無限大に拡大させないために引き算する必要があります（または飽和する）。

Rectangle (矩形波)

矩形波は以下のように計算できます：

$$
R_{\text{ECTS}}(t) = \int_0^t \text{CIT}(\tau) d\tau \leftrightarrow \text{C}_3
$$

これを離散化すると：

$$
\text{Rect}(n) = \sum_{k=0}^n \text{BLIT}(k) \leftrightarrow \text{C}_4
$$

$$
\leftrightarrow \text{Saw}(z) = \frac{z}{z-1}(\text{BP-BLIT}_k(z) \leftrightarrow Z(\text{C}_2))
$$

ここで、BP-BLITは「バイポーラBLIT」であり、パルスが交互に符号を持つインパルストレインです。BP-BLITの効率的な生成方法についてはセクション5で説明します。バイポーラインパルストレインはDC成分がゼロであるため、(\text{C}_3 = 0) となります。矩形波は (k_0) で制御され、パルス幅変調（PWM）を与えることができます。これらの方程式の範囲は ([0, \text{Period}]) です。BLITの実装に応じて、PWM制御は ([0,1]) の範囲にもなります。積分器の初期条件に応じて、出力にはDCオフセットが存在します。

Triangle (三角波)
三角波は以下のように生成できます：

$$
\text{Tri}{\text{CTS}}(t) = \int_0^t R{\text{ECTS}}(\tau) d\tau \leftrightarrow \text{C}_5
$$

ここで、(\text{C}_5) は矩形波のDC成分です。これを離散化すると：

$$
\text{Tri}(n) = \sum_{k=0}^n \text{Rect}(k) \leftrightarrow \text{C}_6
$$

$$
\leftrightarrow \text{Tri}(z) = \frac{z}{z-1}(\text{Rect}(z) \leftrightarrow Z(\text{C}_6))
$$

$${(\text{C}_3)}$$ と $${(\text{C}_6)}$$ は矩形波のデューティサイクルと積分から生じるDCオフセットの関数です。三角波の適切な振幅を得るためには、矩形波の$${ (\text{BLIT})}$$ と同じ振幅にするために、周波数とデューティサイクルに依存したスケーリングを行います。

$$
\text{Tri}(n) = \sum_{k=0}^n g(f,d)(\text{Rect}(k) \leftrightarrow \text{C}_6)
$$

$$
g(f,d) = \frac{2f}{d(1 - e^{-d})}
$$

ここで、(f) は周波数（サンプル単位）、(d) はデューティサイクル（0から1の範囲）です。

4.2 DSFでそれが可能か？

DSF（離散和公式）を使用して、矩形波や三角波を直接生成できるか疑問に思うかもしれません。答えはノーです。なぜなら、積分から示されるように、正方形波の高調波の振幅は全て (1/f) であり、三角波の高調波は全て (1/f^2) です。DSFの高調波は (a^k) であり、これはDSFが正確な矩形波や三角波を生成しないことを意味します。ただし、状況によっては、信号が制御信号として使用されない場合など、(1/n) または 1/n^2 に近似するフィットが十分に近い場合もあります（例えば知覚的に）。

4.3 適切なスケーリング/オフセットと非定常状態の修正

積分器が出力を無限大（あるいは少なくとも数値システムやDACの能力を超える範囲）にしないようにするためには、積分器への入力にDCオフセットを避ける必要があります。既に述べたように、BP-BLITはDCオフセットがないため、矩形波の積分には特別なオフセットは必要ありません。しかし、最初の積分器の初期条件は、出力にDCオフセットを生じさせる可能性があり、これは2回目の積分前にキャンセルする必要があります。このオフセットの値は信号のデューティサイクルに依存します。そのため、正しい初期条件は以下に基づいて変更されます：（1）所望の位相、および（2）所望のデューティサイクル。三角波のスロープが周波数に比例するため、2回目の積分には周波数依存のスケーリングも必要です。パラメータが変更されたときに、古い状態（変更の直前）が積分器の「新しい初期条件」として機能することを覚えておくことが重要です。

リーキー積分器

純粋な積分器の使用は、数値誤差（ラウンドオフなど）の存在下で問題を引き起こすことがあります。これらの誤差は積分器に蓄積し、信号に不要で予測不可能なオフセットを生じさせ、特に2回目の積分で望ましい波形を生成する能力を損なう可能性があります。したがって、積分を実装するフィルタの極を単位円からわずかに内側に移動させます。これらの「リーキー」積分器は、悪い初期条件や数値誤差を徐々に忘れるため、永遠に積み上げることはありません。リーキー積分器のインパルス応答は、ゆっくりとゼロに減衰する指数関数です。これには次の2つの効果があります：

1. 減衰率は周波数に有効な下限を設定します（特に信号が制御信号として使用される場合）。周期が減衰時間のオーダーになると、矩形波などの「保持」出力はもはやそのサイクルの部分に対して一定ではなくなります。オーディオ信号では、これは問題ではありません（最低周波数の高調波を減衰させるだけで、通常は聞こえません）が、制御信号では、信号の正確な形状が重要な場合、大きな問題となる可能性があります。
2. 定常状態では、リーキー積分器の出力にはDC成分がありません（BP-BLITにはDC成分がないため）、初期条件に関係なく、リーキー積分器は最終的にそれらを忘れます。したがって、一時的なDCオフセットがあっても（リーク率で減衰します）、リーキー積分器の存在だけで全てのオフセット問題に対処できます。

それでもなお、DCオフセットの不在（過渡的であれ定常的であれ）が必要な場合、リーキー積分器の存在は適切なオフセットの計算をはるかに難しくします。

振幅の定義

MooreはDSF論文で振幅補償について議論しています。同様の補償がBLIT生成にも必要です。使用する補償は、信号の使用目的に依存します。信号がオーディオ信号として使用される場合、信号パワーや心理音響的なラウドネス測定が適していますが、信号が制御信号として使用される場合、最大値（チェビシェフ）測定がより適しています。

5.DSFによるバイポーラBLITの生成

前述のセクションでは、BLITから矩形波や三角波を生成する方法を説明しましたが、BLITの生成方法については詳しく述べていません。これは、どの技法を使用しても（数値精度の違いはありますが）問題ないためです。ただし、矩形波を生成するために必要なバイポーラBLITを生成するための興味深い方法がいくつかあります。

BLITの差分

標準的（直接的）な方法でバイポーラBLITを生成するには、2つのBLITの差分を取ります。1つはもう一方と比較してシフトされています。これはBLIT-SWSで使用される技術です。

DSFメソッド

DSFが効率的に受け入れられる場合（効率性の観点から）、アルゴリズムにいくつかの改良を加えることで、バイポーラBLITを直接生成するために良い結果が得られます。

まず、DSFの式でsinをcosに置き換えることで、BLITを生成できることに注意します（これにより、元のDSF方程式の分子のsinをcosに置き換え、分母は同じままになります）。図14を参照してください。

まず、DSFの式において負のaを使用することで、正のa信号からちょうど半周期シフトされた信号が生成されることが示されます（図15）。これは、シフトされたBLITをtをオフセットすることなく生成するための効率的（またはエレガント）な方法を提供します。これにより次のことが示されます：

$$
y = \sum_{k=1}^N (0.99^k \cos(0 + 3 k t)), N=20
$$

同じことがcosの和に対しても適用され、図16に示されるように、50%デューティサイクルのバイポーラBLITを生成します。

PWM

他のデューティサイクルの場合、DSFに関するもう一つのバリエーションがあります。複素数aを取り、DSFの実部または虚部を取り、調和波にsin(k/la)（またはcosine）振幅エンベロープを課すことで、櫛フィルタリングに相当し、DSFをシフトされたバージョン（符号反転を伴う場合もある）と合計することと同等です。これらはすべて数学的に示すことができます。図17を参照してください。この方法は、2つの実際のDSF BLITの差を評価するのとほぼ同じ複雑さですが、少しエレガントです。