パンクでアナーキーな（？）MATLAB
まえがき
MathWorks か IEEE かどっち NaN だい！？～ atan2 の定義域を巡って
で、「極まれに MathWorks さん、世間標準に刃向かうパンクなところがある」と書きましたが、「じゃあ他の例は？」ということで、MATLAB と他の環境を行き来する上で注意すべき点をいくつかご紹介したいと思います。
Look out ! 'cause here she comes
・コメント、不等式
これは言語仕様なので仕方ないですが、
コメント
C：/* ...（複数行） */ （C++：// ...、/* ... （複数行）*/）
MATLAB：% ...、%{ ...（複数行）}%
不等式（論理反転）
C：a != b
MATLAB：a ~= b
とバラバラです。まあ、他の言語でもまた違いますが。
不等式に関しては、最近のバージョンの MATLAB エディターでは警告が出て直してくれます。
PC 2 画面で片方で MATLAB、もう片方で DSP 用 C の開発を同時進行してたときはほんとうに混乱して、ちょくちょく打ち間違えてました。
個人的には、できればみんな C と合わせて欲しかったな・・。
・strcmp(strA, strB)
文字列同士の比較関数です。
これも言語仕様ではありますが、C では文字コードを数値として大小比較を行うため、
if strA > strB：正値 ex. strcmp("AB", "AA")
if strA == strB：0
if strA < strB：負値
（正値・負値の値は環境依存）
MATLAB では
if strA > strB：0 (logical)
if strA == strB：1 (logical)
if strA < strB：0 (logical)
です。
これの辺に関しては、C では元々「0 が false でそれ以外が true」と扱っていたり、それなのに strcmp() では一致で 0 が返ってきたりと分かりにくく MATLAB の方が好きですが、論理が逆になるので注意が必要です。
C++ や最近の MATLAB では string 型 （C++ではクラスですが）をサポートしているので、比較演算子等も直接使えとても便利になりましたね。
・biquad フィルター係数符号
biquad フィルターとは、大変よく使われる応用範囲の広い巡回型フィルター構成です。
同じ構成のまま係数を入れ替えることにより、遮断周波数やゲインだけでなく、ハイパス・ローパス・シェルフ・バンドパス・ノッチ・オールパス等、様々なタイプのフィルター特性が実現できます。
ちなみに、巡回型であっても伝達関数の分子が分母で割りきれる場合はインパルス応答が有限となるため、IIR フィルターは巡回型でないと作れませんが巡回型フィルター ＝ IIR フィルターとは限りません。
biquad フィルターの伝達関数$${H(z)}$$は以下の式で表されます。
$$
H(z)=\frac{Y(z)}{X(z)}=\frac{b_0+b_1z^{-1}+b_2z^{-2}}{1+a_1z^{-1}+a_2z^{-2}}
$$
差分方程式を求めるために$${Y(z)=}$$の形に変形すると
$$
\frac{Y(z)}{X(z)}=\frac{b_0+b_1z^{-1}+b_2z^{-2}}{1+a_1z^{-1}+a_2z^{-2}}\\\
\\
(1+a_1z^{-1}+a_2z^{-2})Y(z)=(b_0+b_1z^{-1}+b_2z^{-2})X(z)\\\
\\
Y(z)=(b_0+b_1z^{-1}+b_2z^{-2})X(z) - (a_1z^{-1}+a_2z^{-2})Y(z)
$$
したがって差分方程式は
$$
y[n]=b_0x[n]+b_1x[n-1]+b_2x[n-2]-a_1y[n-1]-a_2y[n-2]
$$
となり、フィードバック側の係数 a には負号が付くことになります。
構成例はディレイの位置によって様々な型がありますが、よく使われる型のひとつ、直接型 II 転置構成（direct form II transposed）で見てみましょう。
こちらは Wiki からの引用です。
この構成図でも、フィードバック側の係数 a は負号が付けられています。
符号付きの係数が乗算器にセットされ、全て加算でフィルター演算が行われるわけです。
DSP 等のライブラリでもこれが多いかと思います。特性によって係数の符号は ± 変わるので、わざわざ引き算を用いる必要性がないからかもしれません。（関数によって異なる場合もありますが）
一方 MATLAB
ん？
んん？
んんん！？
そう、フィードバック側の係数から負号が取られ、代わりにデフォルトで「引き算」になっています。つまり最初の Wiki の例とは、乗算器にセットする係数の符号が異なることになります。
確かにこちらの方が式に忠実とも言えますし見たことはあるのですが（もしかしたら MATLAB を参考に？）、一般的には Wiki のように「（差分方程式に出てくる）負号も係数に含める」場合が多いかと思います。
もちろん計算結果としては同じですが乗算器にセットする係数としては符号が逆なので、「構成」と「係数」が合っていないと大変なことになります。
例えば 1 kHz で +3dB のゲインを持つピークフィルターがあるとします。
sosPEQ = [1.0376 -1.7991 0.77991 1 -1.7991 0.8175;1 0 0 1 0 0];
[h,f] = freqz(sosPEQ,512,48e3);
semilogx(f,20*log10(abs(h)))
grid on
xlabel('Frequency (Hz)')
ylabel('Gain (dB)')
このフィードバック側係数の符号を反転させてみましょう。
sosPEQx = sosPEQ;
sosPEQx(1,5:6) = -sosPEQ(1,5:6);
[hx,f] = freqz(sosPEQx,512,48e3);
semilogx(f,20*log10(abs(hx)))
grid on
xlabel('Frequency (Hz)')
ylabel('Gain (dB)')
当然ながらこのように全く別の特性となり、場合によってはシステムが動かなくなってしまうでしょう。
前に一度、「外注した DSP コードのフィルター部分がおかしいようなので見てもらえないか？」と頼まれたことがありました。チェックしたところ、すぐに「フィードバック係数の符号が逆」であることに気が付きました。フィルター係数設計部分も別の外注で、 DSP コードと構成が合っていなかったのです。
特に大手メーカーでは、20～30年前からすでにメーカー技術者の仕事は「仕様書だけ書いて外注コントロール」になってしまったところが多くあります。
最初の頃は開発経験がある人達がそれをやっていたのでまだ良かったのですが、（当然の結果として）そのうち、入社してから外注コントロールばかりで、開発経験がほとんどない「技術者」が生まれます。
「知識」自体はあってそれっぽい仕様書は見よう見まねで書けても、「検証ができない技術者」が増えてきてしまっているのではないでしょうか・・。
それはともかく最後にもうひとつ。
・正規化周波数
まず、48 kHz サンプリングで、遮断周波数 fc = 10 kHz のバタワースフィルターを設計し、グラフ表示してみましょう。
（fc は、ゲイン = -3 dB となる周波数です）
fs = 48e3;
fc = 10e3;
[b,a] = butter(6,fc/(fs/2));
freqz(b,a,[],fs)
設計通りの特性になっていますね。
この時点で何かに気付かれた方はかなり鋭いです。(-_☆)
freqz() をデフォルト設定で使ってみましょう。
>> freqz(b,a)
x 軸に「正規化周波数 (×π rad/sample)」と書いてあり、右端、つまり fs/2 が "1" です。
通常、どのデジタル信号処理の教科書を見ても、正規化周波数 "1" は "2π" つまり fs であると書いてあります。
少なくとも私は今まで、正規化周波数 "1" を "π" つまり fs/2 とする記述は、MATLAB 以外では見たことがありません。
最初の butter() の二つ目の引数は (0,1) の正規化周波数で指定します。MATLAB の古典的なフィルター関数はほぼそうです。
通常であれば fc/fs で良いはずですが、MATLABではその 2 倍の fc/(fs/2) で指定する必要があります。
freqz() で返り値を使って plot したい場合、
[h, w] = freqz();
の w は (0, 1) ではなく (0, π) で返ってくるので、通常通り $${\omega=2\pi f}$$ の関係から f = w/(2*pi)*fs とすれば周波数に変換ができます。
まあややこしいので、正規化周波数ではなく普通の周波数で返ってくる freqz(b,a,[],fs) の形を使った方が良いかもしれません。
Audio Toolbox ではこの問題を回避するため（？）か、正規化周波数ではなくサンプリング周波数と遮断周波数で指定するようになっています。
これはどうしてこうなっているのかほんと謎ですね。
理由をご存じの方がいらっしゃったら教えてください。
あとがき
思いついたのはこの辺りですが、もし他にもあったら情報をお寄せください。(u_u)
書いていて私も頭がごちゃごちゃしてきたので、誤り等あればご指摘ください。
今年の記事は（おそらく）これで最後になります。
今回のも入れると、今年は MATLAB コードが載っている記事が、去年より 1 本多い 17 本になりました！
みなさま、今年も大変お世話になりました。
MATLAB でもいじりながらよい年をお迎えください。m(._.)m
タイトル画像モデル：夢美