MATLAB の正規化周波数について

MATLAB の Signal Processing Toolbox を使っていると、このページの例のように横軸が「Normalized Frequency (正規化周波数)」になっているグラフが頻出する。しかし、MATLAB の正規化周波数は信号処理の教科書や技術書に書かれているものとやや異なる。そのため、使いはじめたばかりの学生（かつての自分含む）が混乱することがある。

実は、MATLAB の正規化周波数は「ナイキスト周波数」が基準となっている。以下、公式の FAQ から引用する。

This toolbox uses the convention that unit frequency is the Nyquist frequency, defined as half the sampling frequency. The cutoff frequency parameter for all basic filter design functions is normalized by the Nyquist frequency. For a system with a 1000 Hz sampling frequency, for example, 300 Hz is 300/500 = 0.6. To convert normalized frequency to angular frequency around the unit circle, multiply by π. To convert normalized frequency back to hertz, multiply by half the sample frequency. Frequency Response - MATLAB & Simulink

抄訳: Signal Processing Toolbox の単位周波数はナイキスト周波数であり、標本化周波数の半分です。(略) 標本化周波数が 1000 Hz ならば、300 Hz は（正規化周波数では）300/500=0.6 となります。正規化周波数を単位円上の角周波数に変換するには π を乗算し、ヘルツに変換するには標本化周波数の半分を乗算してください。

例

たいていの教科書では、正規化周波数は F = f / f_s と書かれている。ここで f は周波数、f_s は標本化周波数である。同様に正規化角周波数は Ω = ω / f_s = 2πf / f_s = 2πF と書かれている。この定義に基づいた教科書的な IIR フィルタの例を見てみよう。IIR フィルタの伝達関数は

$\displaystyle H(z)=\frac{\sum^M_{m=0}a_mz^{-m}}{\sum^N_{n=0}b_nz^{-n}}$

で与えられる。ここで $z=e^{j\Omega}=e^{j2\pi F}$ として、ある低域通過フィルタの振幅特性を計算すると以下のようになる。

f:id:hark00:20140808193608p:plain:w250 f:id:hark00:20140808193613p:plain:w250

図の横軸は教科書的な正規化周波数である。左図は横軸を 0.5 まで、右図は 1.0 までとした。正規化角周波数に変換する場合は 2π を乗算すればよいので、0.5 は π、1.0 は 2π に変換される。

さて、標本化定理によると、標本化前の信号（原信号）の最大周波数が f のとき、原信号を標本化後の信号から正しく復元するためには、標本化周波数を 2f より大きくしなければならない。言いかえると、標本化周波数が f_s ならば f_s / 2 以上の周波数を持つ原信号は正しく復元できない。この標本化周波数の半分 f_s / 2 をナイキスト周波数と呼ぶ。

冒頭で引用したように、MATLAB の正規化周波数はナイキスト周波数を基準としている。つまり、f_s / 2 が単位周波数となるように設計されている。そのため上図では 0.5 であった横軸の点が 1.0 となる。下図のように。

f:id:hark00:20140808205525p:plain:w250

さきほどは正規化角周波数に変換するために 2π をかけたが、今回は π をかければよい。MATLAB の正規化周波数の単位が ×π rad/sample になっているのは、これが理由である。

なお、MATLAB に似た無料の数値計算ソフトである Scilab の正規化周波数は、自分の知るかぎり教科書通りである。そのため、使用に際して上記のような混乱はなかったと思う。