落胆がらくた街

Dr.Rootはサポート終了製品です。

デジタル信号処理チートシートと、BLITと殴り合ったメモ

デジタル信号処理でよく出てくる式とか変数とか関係をリストにしてみた
シンセとかプログラミングするのにこんがらがった時はどうぞ

  • int R:サンプルレート。殆どの場合、44100か48000のどっちか。ほぼ定数とみなして良い。
  • int F=R/2:ナイキストレート。やはり実質ほぼ定数で、22050または24000。
  • float f:鳴らしたい周波数[Hz]。↓の f_{min}を用いて、 f_{min} \leq f \leq Fになる。
  • float f_{min}:鳴らしたい周波数の下限。大抵10[Hz]とかそんなもん。
  • float T=1/f:周期。単位は秒。
  • float p:現在の位相。2πでループするんだから、 0 \leq p \leq 2\pi。
  • float {\omega}=2{\pi}f:fにおける角速度。つまりf[Hz]のサイン波は現実世界では sin({\omega})=sin(2{\pi}f)。
  • float d=2{\pi}f/R:1サンプル毎に進む位相差。つまり1サンプル処理する毎にp+=d; if(p>=2*PI) p-=2*PI;をするということ。 \omegaのデジタル世界版。fの条件から、 0 \leq d \leq \piになる。
  • int H=(int)F/f=(int){\pi}/d:ある鳴らしたい周波数fにおいて、倍音が何個出せるか(F以下に収まる倍音の数)。整数倍音は全部持つ（ノコギリみたいに）と仮定。
  • float P=R/f:周波数fの信号が1週するのにかかるサンプル数。ぴったり整数個でループするとは限らないので実数。
  • float sin(2{\pi}\frac{f}{R}n):サンプルレートR、周波数fのサイン波のn番目の値。

あとBLITのメモ
\displaystyle \begin{align} Sinc_{M}(x)=\frac{sin(\pi x)}{Msin(\frac{\pi x}{M})}, \\ y(n)=\frac{M}{P}Sinc_M(\frac{M}{P}n) \\ =\frac{1}{P}\frac{sin(M\frac{\pi}{P}n)}{sin(\frac{\pi}{P}n)} \end{align}
ここでMはPを超えない最大の奇数(int)。
こいつを日本語訳すると、 「y(n)はナイキストレートまでの周波数で帯域制限したf[Hz]のインパルス列のn番目の値」。
Pは1周するのにかかるサンプル数なので上記チートシートの通りR/f。つまりサンプルレートを鳴らしたい周波数で割ったもの。
この帯域制限されたインパルス列から帯域制限されたSAWやSQUAREを作れる。
詳しくはg200kg氏のページ参照（必読）。

BLITのお話 | g200kg Music & Software

ところで、おもむろに P=R/f \Rightarrow \frac{1}{P}=f/R をBLITの式に突っ込んでみよう。
\displaystyle \begin{align} y(n)=\frac{1}{P}\frac{sin(M\frac{\pi}{P}n)}{sin(\frac{\pi}{P}n)} \\ =\frac{1}{P}\frac{sin(M{\pi}\frac{f}{R}n)}{sin({\pi}\frac{f}{R}n)} \end{align}
sinの中身をよく見てみると、f[Hz]のサイン波が sin(2{\pi}\frac{f}{R}n)なんだから、 sin({\pi}\frac{f}{R}n)は鳴らしたい周波数の半分の速度のsinだ。つまりBLITとは「f/2[Hz]のサインでMf/2[Hz]のサインを割ったもの」と言える。ここから、
BLIT(p)=\frac{f}{R}\frac{sin(M\frac{p}{2})}{sin(\frac{p}{2})}
が得られる。
y(n)は「BLITのn番目の値」を引っ張ってくるわけだが、BLIT波形の1周が綺麗に整数個で収まるとは限らない（=Pが整数とは限らない）ので、BLIT(p)=「p[rad]におけるBLITの値」を算出できるようにしておいたほうが都合がいい場合もあると思う。

F[Hz]に帯域制限されたf[Hz]のインパルス列が、なんでこの式で得られるのか。
式の算出過程は理想インパルス列(デルタ関数の級数)を理想LPF(時間領域でsinc関数)と畳み込んだものなんだが、得られた式は何がどういう意味を持っているのか。
早速Maximaさんにグラフを描いてもらった。 f:id:suzumodoki:20180121023528j:plain
極限の扱いが上手く出来ておらず、分母が0の時なにやら変な感じになってしまっているが、理想的には滑らかなsincの列なので脳内補完して図を眺める。
まずBLITの式を簡略化すると、振幅を決める \frac{1}{P}はこの際重要じゃないので無視して、意味があるのは \frac{sin(Mx)}{sin(x)}の部分。
ここで x = \pi\frac{f}{R}n だが一旦おいといて、結局こいつは角速度 MxのSinと xのSinの逆数との積だと言ってるわけだ。

まず \frac{1}{sin(x)}はこんな感じ。 f:id:suzumodoki:20180121025414j:plain
赤線が1/sin(x)なんだが、ソフトに書かせると無限の部分を真面目に描いちゃうので、あえて手書きで概略図にしている。
で、こいつにsin(Mx)をかけると… f:id:suzumodoki:20180121030914j:plain
こんな感じ。イマイチ分かりにくいが、大雑把にまとめると

  • 1/sin(x)が発散するタイミングで、 \lim_{x \to \infty} \frac{sin(Ax)}{x}=AよりMになる（インパルスがぴこっと立つ）
  • よって1/sin(x)は「インパルス列そのものの周波数」に対して影響を及ぼしている。
  • よく見るとx=0,π,2πで1/sin(x)が∞になるので、インパルスが立つのはsin(x)の倍の周期。BLITがf/2のsinで割っているのはこういう理由。
  • 分子のsin(Mx)は帯域制限する周波数に対するパラメータ。事実上Mが帯域制限周波数になる*1。
  • もしMが奇数だと、インパルスが立つ瞬間の1/sin(x)とsin(Mx)は同じ符号になるので、必ずインパルスはプラス側に立つ。
  • 逆にMが偶数だと、インパルスがプラスとマイナス交互に立つ。これを「バイポーラBLIT」と言う。
  • 矩形波作りたくてバイポーラBLITをやると周波数が半分になる、というのはこの理屈から1/sin(x)だけ倍の周波数にすれば解決する。多分。

といったところ。

現実的な問題と備忘録。
Sinc_M(0)=1より BLIT(0)=M/Pであるべきだ。黙ってsin/sinをプログラミングすると、sinはともかく割り算が出るので一々0除算じゃないかチェックしなきゃならず辛い。
この間「sinをテーブル引きにすると速い」というネタを自分の環境で試したら速いどころか100倍以上遅くてゲンナリだったんだけど、Sincとか1/sinみたいなのをテーブル化するんならやぶさかでない気もする。というわけで今後は本格的にどんな実装が速いのかとか、数学よりもっとC++とかx64で効率的なやり方を考えようと思います。どっとはらい（あ、あと位相変数は2πをUINT32に正規化するとオーバーフローのお陰で周期性チェックいらなくて便利じゃねっての思ったのでこれも今度試します）

*1:うまいことやればPを整数に決め打ちしても任意のMで帯域制限できそうなんだけど勉強不足で未検証