DFT IPの作成 13

どうやら他にも問題があるようで、入力信号の振幅が最大値の75%を超える値がある場合に波形が変になる。具体的にはスペクトル波形が残ってしまうという以前の症状と同じになる。
どうもおかしい。。。そもそも、方式的に無理があるのか？？

そこで、Cモデルを改良してwavデータを処理出来るようにした。以下がCモデルのDFT演算処理を行う部分で極めてシンプルだ。RTLも基本的にこれと同じ筈なんだがなぁ。。。


これで実機と同じデータを処理させて実機と同じ結果になるかを確認する。基本的に同じ処理を行っている訳だから、同じ結果が得られるはずだ。以下が使用法だ。


my_dft -i wavファイルへのパス  と打つと、指定ファイルのDFT演算を行い結果をファイルに出力する。 -o オプションで出力する項目を指定する。


-r オプションは出力レートで、0 (default)の場合は2048点演算する度にファイルを出力する。一方、1の場合は1点演算する度にファイルを出力する。

defaultの出力形式 (-o 0)で出力されたファイルをgnuplotでグラフ化するシェルスクリプトも作成した。


以下が操作例だ。


my_dft -i wavファイル名を実行すると、カレントディレクトリに出力ファイルが生成される。拡張子は.spcとした。この状態でpic_gen.cshを実行すると各出力ファイルのデータのグラフを作成しjpeg形式でファイルに出力する。これはjpgというサブディレクトリを作成してその下に出力される。このjpegファイルはavidemuxというソフトを使ってavi形式の動画ファイルに変換することも出来る。このavidemuxというソフトは優れもので、mp3やwav形式の音ファイルも結合させることが出来る。さらに、これは別のソフトだが、ffmpegを使えばこうやって生成したaviファイルをmpeg形式の圧縮ファイルにすることも出来る。

以下はその例
低域側の200ポイントの部分のみをグラフ化している。
※ 音源は「甘茶の音楽工房」というサイトのロイヤリティフリーの音源を使用させて頂いた。


この方法で、動作が異常になるwavデータを処理した結果が以下だ。

同じデータを実機で再生させると以下のようになる。


明らかに実機の方の動作が異常だが、ソフトウェアの方は正常に処理できているので、再帰的に演算を行うという本方式自体には問題なさそうだとも言える。では、一体何が間違っているのか？？？

シミュレーションで確認したいのだが、Xilinxの浮動小数演算器のverilogソースがゲートレベルになってしまっているために非常に時間がかかると言うことを先日のブログに書いた。DFT IPの単体シミュレーション環境でなら少しは短縮できるか？と期待してやってみたのだが、サンプリングレート48KHzの2048点分に相当する42.6msecの時間分のシミュレーションに43時間を要するという有様だ。(;_;)　波形ファイルはFST形式でダンプさせているが、サイズは25GBにもなった。まぁ、波形ファイルはダンプする対象範囲を狭めたり、$dumpon/offを使えば小さくできるかも知れない。問題はシミュレーションの実行時間だ。これを何とか高速化する必要がある。そこで、浮動小数演算を行うVPIタスクを作成しこれを使ってみることにした。演算器を置き換えるため、不具合の原因が演算器にある場合は不具合現象が再現できなくなるが、逆にこれで正常動作が確認できれば、原因は演算器か或いはもっと別の要因、例えばACタイミング的な事等と切分けできるかも知れない。

以下がVPIタスクのプログラムの一部で、加減乗除の2項演算用の関数だ。この他に、今調査したい現象とは関係ないが、sqrtとlogの単項演算用も作成した。


verilogの方は以下のように書き換えた。


で、実行時間の方だがすんげー早くなった。263msecをシミュレーションするのに90分程度だ。
元は16.5usecのシミュレーションをするのに1分を要する計算だが、VPIタスク化した場合はこれが2.9msec/分と177倍程も高速化できたことになる。後は、ダンプファイルが肥大化しないように工夫して、何とか不具合の原因が特定できるといいのだが。。。

DFT IPの作成 12

例の誤差問題だが、浮動小数演算部の丸め誤差の累積によるものだろうと考え、縮小係数を掛けることで妥協することにしたのだった。　前々回のブログで発生する誤差について式を作って検討した上でそうすることにしたのだが、実は釈然としない感じをずっと持っていた。
が、原因が判った。
原因は、xp - xp-N を計算するモジュール dft_dlyline.v内のポインタのズレだった。

以下は、修正した回路による単一正弦波信号のスペクトルの絶対値表示だが、完全な線スペクトルになった。



また、音楽データを再生させてもスペクトルが残るような事もない。

やった。！！　これで理論どおりだ。
妙な誤魔化しの手法も要らない。

修正は多分合っていると思う。
時間が無くなったので、続きはまた明日以降にしよう。

DFT IPの作成 11

DFT IPの作成をもう少し続けることにした。

前回のブログにも書いた通り、基本演算部の固定小数化をしてみたのだが、bit幅が48bitと広すぎたこともありなかなかタイミングがMetしない。手動で合成条件を色々変えてみたり、ISEのSmartXplorerで１０数パターンで合成させてみたりして何とかMetさせた。が、残念ながらスペクトルは誤差の影響が見られた。現状の浮動小数演算方式の方がまだましといった感じだった。
この件は一旦保留とすることにした。

元々このIPはAC97コントローラに繋いで使う（実際、今そういう使い方をしているが）つもりだったので、チャネルの2ch化（ステレオ対応）をすることにした。現行AC97コントローラの最高のサンプリングレートは48KHzであるのに対してIPは100MHzで動作させているので、１サンプル当たり2083点分のDFT演算が出来るので2ch分の演算は可能である。実際、回路の変更も非常に簡単で、変数用RAMの容量は倍になるが演算部は変わらない。あとは、タイミング等を若干変更するだけだ。

バイクの通過音のデータを再生させてみた。
2ch分のスペクトル表示をするために、表示部の回路も若干変更した。具体的にはch1（R）を赤色、ch2（L）を緑色で表示することにし、また、これまでは棒グラフ式表示だったのを線グラフ式にした。

ピアノ演奏のデータを再生させてみた。


上記データのログ表示モードでの再生
データに強いピークが含まれていないと、ログ表示は面白みがない。

ということで、2ch化は容易にできた。
次はスペクトルデータをDRAMに書込む機能を追加したいと思う。

DFT IPの作成 10

DFT IPのsnapshot版を下記に置いた。


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


ブログの文章や画像だけでは、動いてもいないものを動いたとしたり、出来てもいないものを出来たとしたり等いくらでも誤魔化せるので、このブログで作成したものは、成功も失敗も上手いも下手もひっくるめて、エビデンスとして基本的に公開することにしている。

DFT IPの作成を後どれくらい続けるかは思案中だが、今現在は、基本演算部を固定小数演算方式で実装してみようとしている。 

タイミングがなかなかMetしないんだこれがまた。。。

DFT IPの作成 9

このところ、寝ても覚めても誤差のことで頭がオッパイだ。いや、イッパイだ。

今回のアルゴリズムの式は以下であった。

この式で、簡略化のためにxp-xp-Nをxp、exp(j2πf/N)をkとし、またこの演算で発生する誤差をεp(複素数）として以下のように表す。

すると、この式は以下のように展開できる。

さらに、書き直すと、

となる。したがって、誤差成分は以下のようになる筈だ。

kはexp(j2πf/N)なので値域は-1.0～+1.0だ。　εは浮動小数演算の仮数部で発生する丸め誤差だが、大きさは指数部の値によって決まる。入力信号は16bit整数であり、これまでは、これをそのままの大きさで演算していた。この場合、指数値は2^-15から2^15の範囲となる。これを最大値(32768.0)で正規化すれば、exp()とスケール感が同じになるので誤差も小さくなるかも？と考えた。しかし、これは単に縮尺を変えているだけとも言えそうなのであまり違いはないのかも知れない。そこで、実際どうなるか、単一周波数のスペクトルで見てみた。

以下は正規化していない場合（ログ表示）


以下は正規化した場合
。

やはり、違いはなかった。

では、乗算や加減算で丸めをせずに切り捨てしたらどうなるかを見てみた。coregenで生成する演算器は丸めの有無を選択できないため、この確認の為には乗算器と加減算器を自作する必要がある。

乗算器は以下のようにした。


加算器は以下のようにした。


減算器(y=a-b)は、基本的に加算器のRTLを流用しb入力の符号を反転して使った。

結果としては、乗算器のみ丸め無し版に置き換えた場合は結果にほとんど違いがなかったが、加減算器の場合のはより悪くなり、スペクトル波形が激しく振動したような、無茶苦茶noisyな波形になってしまった。

ということで、誤差対策としては、縮小係数をかける方法と、一定回数毎にアキュムレーションRAMをクリアする方法のどちらかしか出来なさそうだ。他に変数の桁数を広げる（単精度小数を倍精度にする。）というのも考えられるが、回路規模がデカくなる、タイミング的(100MHz動作）に厳しそうだ、桁数を広げたとしても量子化された値での演算である以上は誤差は必ず発生し累積されてくる、等々考えて桁数を広げる案は止めにした。

では、2者の内どちらが妥協できるか？ということで、それぞれの場合の波形を比べてみる。

2048点毎にアキュムレーションRAMをクリアする方式の場合は以下のようであった。


見事に線スペクトルになっているように見えるが、上記はログ表示であり、これを絶対値表示（√r^2 + i^2）にすると以下のようであった。 （ちょっと変だ）


一方、縮小係数をかける方式の場合は以下のようであった。

ログ表示


絶対値表示


と、言うことで今回は縮小係数をかける方式で妥協しちゃおうかと考えている。

学生のときにもっと数学を勉強しとくんだった。 (;_;)

DFT IPの作成 8

RTL全体でシミュレーションしてみたところ、以下のような結果だった。


上記は開始から約10msecの時点での波形だが、右端に小さなピークが見える。
入力信号は上図の左側のスペクトルに対応する単一周波数の正弦波なので、右側のピークはおかしい。　さらに波形とRTLを追って見たところバグっている箇所が分かった。入力信号は1024点の演算実行期間、値を保持していなければならないのだが、それが抜けていた。
（単体シミュレーションではテストベンチ側で値を保持していた。）


上記波形でxpが入力信号であるが、値が保持されていない。
これを以下のように修正した。


xdはxp - xp-Nの結果で、16bit整数値で演算してその結果を32bit浮動小数に変換しているが、この変換器でstb信号を入れて変換結果を保持させるようにした。
この修正の結果、表示されるべきスペクトル（上図の場合は左側）のみが表示されるようになった。

但し、誤差問題は解決していない。
これまで、Cモデルやシミュレーションでは発生せず実機動作で発生すると考えていたのだが、これは誤りで、恐らく、Cモデルやシミュレーションでも誤差は発生している。これが見えなかったのは実行時間が短すぎたためと思われる。シミュレーションの場合、私はIcarus verilogでシミュレーションをしているが、全体シミュレーション環境で10msec程度をシミュレーションするのに１２時間位要している。シミュレーション時間はRTLのノード数・・・というか回路の規模に比例する。浮動小数乗算器や加減算器等はcoregenで生成したモデルを使用しているが、これはネットリスト形式になっているので、演算器の部分はゲートレベルシミュレーションしているようなもので、これもシミュレーション時間が長時間になる要因になっているだろう。しかし、誤差問題はこの演算部分に関わることなので、ビヘイビアモデルを作って置き換えるという訳にもいかない。何れにしてもシミュレーション時間が短すぎる為に、誤差問題を見落としていたのは間違いないだろう。

で、この誤差のキャンセル方法について、これまでに考えたのは次の2つだった。
1. 2048点演算毎にアキュムレーションRAMをクリアする。
2. 縮小係数をかける。

結局、1の方法で普通に2048点毎の演算としたほうがいいのかも知れない。

DFT IPの作成 7

乱れている～、乱れている～、検量線っがー、
おーかしいよとー、言われるけどー、しーかたがない～のー

っていうフレーズが何故か頭の中でリフレインしている今日この頃如何お過ごしでしょうか？
あれは、黒木真由美の感情線だっけか？、懐かしいなぁー、あの頃は俺も若かったぁ、
色んなとこブヒブヒ言わせてさぁー、、、、太ってたからねぇー
それが、今じゃ、こんなに痩せ細っ・・・、あれ？、細ってない。むしろ増えてる！！！　おっかしいなー

なーんて言うおバカキャラはこれぐらいにして、DFT IPのデバッグだ。
いやー、ストレスが溜まってるときはこの位バカして発散しないと、ねぇ?

元々考えていた方式でやると、入力がゼロになってもアキュムレータRAMに値が残り、時間経過と
共に増加していくという問題が出た。以下はその様子だ。

画像では、入力値がゼロになってからのスペクトルの成長が分かり辛いが、じっと眺めてると、
所々で少しずつ成長しているのが分かる。
入力値がゼロになっていく過程で、アキュムレータRAMの値も減少していく筈なので、値が残るという事自体おかしな事なんだが。

やっている事（回路で処理している事）は、


という演算で、入力値ゼロの状態が継続すると、Xp-Xp-Nは0-0になるから、S(f,p)は増えもしなけりゃ、減りもしない筈で、正に何も足さない・何も引かない状態になる筈だ。が、その絶対値であるスペクトルが徐々に成長するってことは、上記演算過程で何らかの要因でプラス方向の誤差が発生していると考えられる。また、S(f,p)は複素数であり、これに複素数の定数が掛かる訳だから、exp(jα) x exp(jβ) = exp(j(α+β))な訳で、複素平面で考えると、ベクトルを常に回転させていることになるが、これが誤差により外側に膨らんでいっていると考えられる。

そこで、S(f,p)に1.0より小さい縮小係数を乗じて見ることにした。但し、その為に態々乗算器を追加するのは回路の無駄なので、正弦波テーブル(exp(j2πf/N)のテーブル）の値を縮小することにした。具体的には以下のプログラムでテーブルを作り直した。

0.9999fが縮小係数だ。 上図の通り、これ以外の値も試してみたが0.9999より大きいと効果が弱かったり、まったく効かなかったりだった。

このテーブルを使って上記と同じデータを再生してみたのが以下の画像だ。

でも、良く良く考えてみると、何かおかしい。 Cモデルでの検証ではこういう現象は出てなかった。
32bit浮動小数点データは仮数部が24bit（省略されている最上位の1を含む）のみであり、位の違う数値同士を演算した場合、当然ながら、その有効桁の範囲外に値が出ることもある。この場合に丸め処理を行う訳であるが、これは「最近点への丸めモード」としてIEEE754で規定されている。っていうことが色んな文献に書いてあるが、CQ出版の「ディジタル数値演算回路の実用設計」という書籍の表6-2がとても解りやすい。　この丸め処理がcoregenで生成した浮動小数演算器と、Cの処理系とで違うのかな？と考えて検証してみたが同等だった。

それではということで、単一正弦波のデータを作って再生させて線スペクトルになるかを見てみた。


その結果が上記画像で、線スペクトルになっているが、右端にも想定外のスペクトルがいる。
DFTは2048点で演算していて、その内の半分の1024点を表示しているので、折り返しは表示されない筈だし、Cモデルによる検証でも、また、DFTモジュールの単体シミュレーションでも表示されていない。
おまえ・・・誰？ (-_-;
ってな感じだ。

どっかバグってんのかな?、全体でシミュレーションやってみるか。

statistics

今日のページビューの中で、これは日本からだが、「dft+面白い」という検索から、このブログに来た方がいるようだ。。。これはウケた。「dft+おばか」だったらもっとウケたのに。


このブログのページビューは今年の夏頃一時的に1000件/月を越え、その後、ブログの更新の停滞と共に減少していたのだが、ここ最近増加傾向にある。不思議なのが、最近はアメリカからのページビューが増えている。日本語で記述している上に、それでも読みづらい文章の筈なのに、不思議だ。インタネットを自動で検索して回るソフトウェアによるページビューかも知れないな。



ちなみに、ブログを開始してから全期間の統計は以下の通りだ。
bloggerのシステムでは国名は上位10迄しか表示されないようだ。
過去、ここにリストされていない国名を見た記憶があるので実際はもう少し多いと思う。


内容としてはYUV422について書いたのがもっとも多く参照されている。

YUV422に関しての情報は意外とネットに転がっていないのかも知れない。
私の文章も大して中身は無いので、その意味では申し訳ない。

そうそう、統計つながりで、dviencのダウンロード数がもうすぐ300に到達しそうだ。
去年の10月に200だったので、１年で100件程の増加だ。


国数は32ヶ国に増えた。




おそらく絶対に有り得ないとは思うが、全大陸制覇したら、俺、興奮しちゃうね。

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