DFT IPの作成 3

DFTで求まるのは実数部と虚数部からなる複素数であり、スペクトルとして表示するためにはこれの絶対値(sqrt(r^2 + i^2))に変換する必要がある。  また、通常、対数表示されるので対数変換も必要になる。  この対数変換器について検討した。
まず、対数は以下のように任意の数を指数形式で表した場合の指数部のことであり、a(底)が10の場合を常用対数、ネイピア数の場合を自然対数と言う。スペクトルは常用対数で表す。

今作成しているDFT IPは内部演算を32bit単精度浮動小数形式で行っているが、対数変換を考える場合この浮動小数形式は都合が良い。  32bit単精度浮動小数は数値を±1.M×2^e 形式で表す。構造は以下のようになっている。

sは符号で0が正数、1が負数、exponent(e)は指数部、Mantissa(M)は仮数部である。
指数部は本来は整数で+/-の値を取りえるが、この形式では+127のオフセットが加算され正数となるようになっている。仮数部は整数部（1または0)が省略され小数部23bitのみが含まれる。
また、指数部と仮数部の特定の値の組み合わせに対して以下のような状態が割り当てられている。

このDFT IPでは値が不正規化数や無限大、非数になることは有り得ないので、正規化数に限定して検討を進める。 また、符号は負になることはないので正に限定して考える。

であり、仮数部の対数値(m)とオフセットを引いた指数部の値を加算すれば対数値が求まる。
次にmについて考える。

であり、グラフにすると以下のようになる。

凸型の曲線になっているので、計算して値を求めようとすると多項式か何かでの近似計算となりそうであり、それよりは、テーブル参照方式で求めたほうが良さそうだ。その場合、Mの小数部の上位何bitかをインデックスとして用いることになる。4,6,7,8bitの場合について見てみた。

当然だがbit数が多い程誤差は少いが、今回の用途では6bitでも良さそうな気がする。
つまり、64x32bitのテーブルを作成することにする。  対数変換器への入力は二乗和の平方根だったが、sqrt(x)はx^1/2であるので、開平せず二乗のまま対数変換をしその結果を2で割るという方法も可能で、こうすればsqrt演算器を省略することができる。ここまでの処理で求めた値は2を底とする対数値なので、最後にこれを常用対数に変換する必要がある。これは(log210）^-1をかければいい。ここでもlogが出てきてしまうが、固定値なので定数化できる。
と、言うことで、ここまでの処理イメージをまとめると以下のようになる。