DFT IPの作成 6

実機で動かしてみた。

現状のSystem Architectureは以下の通り。

AC97コントローラの再生のRchをDFTに入力し、DFTの出力をDRAMを介さず、直接VIF(フレームバッファ)に入れている。  VGAの表示サイズがXGA (1024x768)で横方向がDFTの出力データ数と一致するので、VIFのsyncgenのxadr(水平方向カウンタ）の値をDFTのpramのアドレス入力として与え、その出力をsyngenのyadr（垂直方向カウンタ）の値と比較して、hitした場合に赤、hitしない場合に黒を出力するようにして、スペクトル表示を行うようにした。
以下はDFTの当該部だ。


また、以下はVIFの当該部だ。


この構造で音楽データを再生させながらスペクトル表示させて見たところ、再生終了（音量がゼロ）後もスペクトルが残ったり、また、そのまま放置しておいても、スペクトルが徐々に増加したりといった現象が見られた。どうも、演算誤差の累積が原因のようだ。　そこで、とりあえず、データ入力2048点毎にアキュムレータの値をクリアするようにして、動作確認を進めてみることにした。これは2048点毎にDFTを実行することと同じことになるが、この回路の設計思想からは外れるので、誤魔化しではある。dft_coreの当該部を以下に示す。 165行目のscntでデータをカウントし、2048個に到達するとsclrが1になり、アキュムレータRAMからのリード値を0に置き換えている。(180,181行目)


回路規模だが、現状 dft_core部は以下の通りだ。


dft_core部の中身はDFT演算のみであり、絶対値演算と平方根、対数変換はdft_powerという別のモジュールで行う。この、dft_power部の回路規模は以下のようになった。


現状、浮動小数演算器はハードウェア乗算器(DSP48)を使わない構成で生成しているため、LUTの使用量が大きくなっている。
dft_top (dft_core + dft_power)部としての規模は以下のようになった。


これまでの数値は、論理合成結果なので、マッピング後の結果とは若干違ってくるかも知れない。

また、何とかこれをSpartan3E 250Eにも入れられないかとやってみたのだが、駄目だった。
ムググッ、おのれ250Eめ、今度固定小数方式でリベンジしてやる。
（その前に、誤差問題を対応する必要があるが。）

Spartan6 xc6slx45への、全体の合成・マッピングは出来て、タイミングもMetした。


で、データを再生して表示してみた。

これは、絶対値を平方根した値を表示している。
これの対数変換による表示は以下だが、ちょっとパッとしない。

サイレン音を再生・表示してみた。サイレンの音程の変化に合わせて、スペクトルが移動する様が、ちょっと面白い。

次に、自動車の通過音をやってみた。

音楽をやってみた。

きゃ～～、>_< アレサー、、おしっこチビりそう。

今度はテンポの早いやつをやってみた。

いやー、めでたい。ウチナーンチュ最高！ (俺、ウチナンチュでっす。)

これで、誤差問題が解決すればもっとおめでたい。
DFT演算ではexp(2πf/N)を乗じているので、ベクトルを常時回転させていることに相当するが、入力がゼロでもスペクトル値が増加すると言うことは、演算誤差により回転時に外側に膨らんでいっていると考えられる。だから、内側への補正が出来れば解決出来そうな気がするのだが、、、要検討だ。

げっ！、もうこんな時間だ。寝なきゃ。

DFT IPの作成 5

この土日は色々と野暮用があって作業が捗らなかった。
Nを2048にし、演算はその半分の1024のみを行うようにし、また、対数変換器を作成し
RTLに組み込んでシミュレーション実行までは出来た。
以下、上段は絶対値のスペクトルで、下段が対数でのスペクトルだ。


あ、それと、各種浮動小数演算器のレイテンシーも現実的な値で生成し直した。
具体的には乗算器・加減算器ともに4にした。
もう少しで、実機での動作確認に移れると思っており、AC97コントローラに組み込んで動作させてみるつもりだ。　が、ここで、対数値を20倍すべきか10倍すべきかが心許なくなってきた。。。
表示しようとしているのはパワーか？電圧か？で変わってくるはずだが。
こま～～～～けぇ～～～こたぁ～～～、いいから一度実機で動かしてみようか。
落ち着いて考えよう。

DFT IPの作成 4

今回は位相角は表示するつもりはないのだが、面白いのでちょっと検討してみた。
位相角はDFT出力のsin,cosから以下のように求める

ここで、sinθ,cosθがそれぞれ以下のような指数形式で表せるとすると、

となり、除算がなくせる。
そこで、tan^-1を例によってテーブル方式とする場合、インデックスはα^βになる訳だが、
α^βはβから一意に定まるので、インデックスとしてβを使うようにしても良いはずだ。
即ち、

という訳で、前回検討した対数変換器でsinθとcosθの対数を求め、その差分をインデックスとしてテーブル参照すればθも求められるのではないかと思う。
但し、0°や90°等の特異点は別に判定する必要があり、処理イメージとしては以下のような感じか?

この方式はテーブルの持ち方が肝になるような気がするのと、精度や回路規模、スループット、レイテンシー等の面でどれほど実用性があるかは疑問（ちゃんと考えていない）が、sinθ/cosθの除算をしなくて済むというのは面白い。

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が出てきてしまうが、固定値なので定数化できる。
と、言うことで、ここまでの処理イメージをまとめると以下のようになる。

AC'97 Codecを制御してみる 5

AC97コントローラのsnapshot版を以下に公開した。
（今作成中のDFT IPの前にやっていたやつです。）

http://www.hi-ho.ne.jp/bravo-fpga/

DFT IPの作成 2

RTLだが、まずは1024点版を作ってみている。 演算部は当初、整数部16bit、小数部16bitの固定小数方式でやってみたのだが、累積誤差が大きくなるので結局、浮動小数方式でやることにした。

現状のRTLは以下の通りだ。ただし、まだ最終型ではない、シミュレーションでの動作確認用と言ったところだ。

dft_core.v ... ipコアのトップモジュール

5行目はステートマシンで、rst解除後の内部変数（テーブル類）の初期化と、演算回数の処理を行う。　mul_fp32, add_fp32, sub_fp32, は浮動小数演算器で、ISEのCoregenで生成した物を使っている。今回は動作確認の段階なので、レイテンシーは1にしているが、100MHzで動作させる場合は、おそらく、3～4位にする必要があるだろう。従って最終的にはパイプラインの段数が増加する。

dft_dlyline.v ... xp, xp-N用の遅延モジュール


dft_arctbl.v ... sin,cos値のテーブル

1/4周期分のみromに持ち、残りの部分の値はこれを写像して生成している。
romの中身は32bit float形式だが、このデータは以下のようにプログラムを書いて作成した。


dft_spctbl.v ... スペクトル値保持用モジュール


実態は只のRAMだ。

RTLとしては浮動小数演算器を除いては以上で、結構シンプルに出来ていると思う。
但しスペクトル値にする為には、さらに乗算器2個と√演算器、log演算器等が必要になる。

以下はシミュレーションの様子だ。  基本周期の64倍の正弦波を入れている。


上記の拡大


現状1024点としているがスペクトルの折り返しがあるので使えるのは半分になる。
従って、有意な1024点を得るためには計算は2048点で行う必要があるが、DFTの場合、必要な分のみを計算できるので、Nが2048点でありながら計算は1024点のみを行うようにも出来る。
また、スペクトル波形を画面表示することを考えるとこの位の点数で十分な気がしている。


いや～～っっっ、 FPGAって本当に面白いですね !!

DFT IPの作成

AC97 Codec用のコントローラは大体出来たので、次にオーディオ帯域用のスペクトルアナライザを作って見たくなった。スペクトルアナライザと言えば、おそらくFFT方式が常套手段だと思うが、今回はDFT (Discrete Fourier Transformation)方式でやってみることにした。

DFTの式と言えば、教科書等には大体以下のような式が記述されていると思う。

この式は一つの周波数に関しての式になっており、全体的なスペクトルを求める形式にすると以下のようになると思う。

ここで、信号Xnについて、信号は-∞から+∞へ連続している信号列と定義し直す。

つまり、

そして、(1)式を、任意のp点からp+(N-1)までの区間の信号列に対して行うように書き直す。

また、今回はスペクトルアナライザを作成するので、p～p+(N-1)の区間の次は、p+1～p+(N-1)+1、…と連続して変換していきたい。

この(4)式はS(f,p) を使うと以下のように書ける。

ここで、

なので、結局、

となる。。。。 うほっ、、、これは面白い。！！
この式は、直前のスペクトル値を用いて回帰的に計算すれば、exp(2πf/N)の乗算だけで済むことを表している。乗算・加減算の回数も(1)式に比べてかなり少なくなる。

ここまでの式では、p点からp+(N-1)までの信号を用いて、p点のスペクトルを計算するという表現になっていたが、これだと未来の信号を用いて現在の値を求めるような形になり、少々気持ち悪い。

そこで、過去のN点から現在の値を求める形式に書き直す。

これを、実数部と虚数部に展開してまとめると以下のようになる。

(1)式の場合、乗算は2N回、加算が2(N-1)回だが、上記式の場合、乗算4回、加減算が6回だ。

念のため、C言語のプログラムを作成してシミュレーションしてみた。
プログラムは以下のとおり。

32～41行で入力信号を作成している。以下に波形を示す。

0～1024の区間はx1、x128倍、x256倍の周波数の合成波で、2048～3062の区間はx64倍の周波数の信号だ。従って、処理の初めのほうでは、スペクトルはx1,x128,x256の３箇所で観測され、次にx64で観測され、最終的にはスペクトルは観測されなくなる筈だ。
上記プログラムの出力から、各サンプル毎のスペクトルをgnuplotで描画・出力し、avidemuxでこれらを合成して、動画を作成した。このアルゴリズムをFPGAに実装して、上記のような信号列を入れれば以下のようなスペクトル波形が得られるはずだ。

うほほっ、、、面白い。 ちゃんと周波数領域に変換できている。
この方式だと、コンパクトな回路が出来そうだ。
使用周波数帯は、オーディオ帯域用を考えているが、これに対してFPGA内部が100MHzでの動作となるので、例えばオーディオのサンプリングレートが48KHzの場合速度比は2083倍になるので、１サンプル毎に2048点程度のDFTは楽勝で出来そうだ。

スペクトル全体を求める式は(2)式風に書くと、以下のように書ける。

ここで、正弦波形の対称性を考慮すると、

なので、(9)式は以下のようにも出来る。


つまり、exp(2πf/N)部は1/4象限分のみ持てばよくて、かつ、演算部を4並列することで、１サンプル以内の点数を4倍に増やせる。したがって、回路規模は若干大きくなるがこの方法でやれば8192点程度は出来そうだ。

それでは、RTL作成に進もう。

AC'97 Codecを制御してみる 4

DRAMコントローラとクライアントモジュールとのI/FはMPIFという独自I/Fになっており、これは通常のREAD,WRITEコマンドの他にFILLとSWAPコマンドがある。


SWAPは同一アドレスへのREADとWRITEをアトミックに実行できるコマンドで、マルチCPU構成時のMutex制御（排他制御）等での使用を想定して設けたコマンドだが、それ以外にも使用できそうだ。

今作成中のDMACはREADを使っているが、SWAPを使うようにすればDMAバッファをリードした後にゼロクリアしたり、あるいは、マイクやLINE IN等の録音データで置き換えることができる。
バッファをcirculationモード、即ちリングバッファとして使用すれば、録音データの遅延再生機能を実現できる。LM4550は下図に示すとおり入力データをDACデータと合成して出力する機能もあるので、この機能と合わせるとリバーブレータ（残響効果器）のようなことができそうだ。



レジスタ仕様は以下のように変更した。


以下はシミュレーションの様子 (clear_to_zero 選択時）


replace_to_ADC_data モードで、遅延が250uS, 10mS, 100mS, 500mSになるようにして、動かしてみた。

250uS

10mS


100mS


500mS

AC'97 Codecを制御してみる 3

AC'97 Codec用の制御回路を作成して、信号を出したり音楽データを再生させることはできた。
元々はCodecをDACやADCとして使ってみようと考えていたのだが、意外と出力のノイズが目立つので方針を変更して、本来のAC'97用コントローラの形にまとめてみることにした。

以下にブロック図を示す。


このコントローラは3つのブロックで構成される。 CodecへのPCMデータの出力はDMAによる方法とCPUによるレジスタへの書き込みによる方法を選べる。
DMAの場合、DRAM上にバッファ領域を確保しそのアドレスをレジスタに設定する。


DMACにはバッファの使い方で2つのモードを持たせた。１つはBEGINからENDまで転送したら動作完了となる non-circulation モードで、もう1つはENDに達したらBEGINから再び転送を行う circulation モードである。 前回の実験では音楽データ1曲分をDRAMに展開して再生させたが、これだと非常に大きなメモリ空間を占有してしまうことになる。　音楽データ全体はフラッシュメモリや外部記憶上にあり、そこから少しずつDRAM上のバッファに転送して出力するというやり方が一般的だろうと思う。  ただし、この場合はソフトウェアの介在が不可欠で、DMACのバッファからの読出しに同期してバッファのデータを更新する必要がある。このコントローラの場合は、DMACのポインタがTHRESHOLDレジスタで示すアドレスを通過したタイミングと、ENDレジスタで示すアドレスを通過したタイミングで割込みを発生できる。THRESHOLD通過割込みが発生した場合、BEGINからTHRESHOLD間のデータは転送済みなので、ソフトウェアはその領域を次のデータで更新する。同様にして、END通過割込みではTHRESHOLDからEND間のデータを更新する。もちろんこの方法の場合はDMACをcirculationモードで使用する。

全体的な外部仕様（レジスタ仕様）は以下のとおりだ。