2017年4月
FPGAにズームFFTを実装
2017.04.28
FFTで周波数を測るとき、周波数分解能を上げたい場合、どうすればいいでしょうか？
周波数分解能は、fs÷サンプリングポイント数ですから、例えば5MHzサンプリングで65536ポイントのFFTを行うのであれば、分解能は約76Hzとなります。
では、さらに分解能を上げたい場合は、どうすればよいでしょうか？
簡単に思いつく方法は、一つはFFTのポイント数を増やすことです。
262144ポイントのFFTにすれば、周波数分解能は19Hzになります。
でも、FFTの計算量はN×log2(N)ですから、ポイント数を増やすと計算量が増えて効率的ではありません。
そこで、ズームFFTというものがあります。
原理は、下の図をご覧ください。
ズームFFTというのは次のような手順で行います。
- ノイズの多い入力信号がADCから入る
- BPFで目的の周波数付近を切り取る
- 正弦波を掛けて周波数をシフトさせる
- LPFで低いほうの周波数だけ残す
- デシメーションを行って、データ数を減らす
原理はこのような感じなのですが、実際に実装するにはBPFやLPFの設計、正弦波発生回路が肝になります。
まず、BPFですが、FPGAの中に入るフィルタではFIRを作るのがとても簡単です。FIRというと乗算器と加算器でたくさんの信号を足し合わせるのですが、まともなフィルタを作るには数十から100次程度のフィルタにしなければなりません。
たくさんの加算を行いますが、クロック速度に対してデータレートが遅い場合にはリソースを大幅に節約することができます。
125MHzクロックで、5MHzの信号を処理する101次フィルタの周波数特性は下の図のようになりました。
101次にもなると、たくさん48bitのDSPブロックを使うんじゃないの・・・？と思うかもしれませんが、実際に使用したDSPブロックは4個程度でした。
すっごーいお得ですね。
正弦波はDDSのような感じで位相信号を生成しますが、SINやCOSの波形を作るにはCORDIC法を使います。CORDIC法はXILINXのLogicoreを使えば難しくありません。
このコアは「位相」の信号を入れればSINとCOSが出てくるというものです。
入力する位相の情報をScaled Radianにすれば、1がπラジアン、-1が-πラジアンとして解釈されます。ただし、位相の情報は整数部が2bitある固定小数点数なので注意が必要です。
下の図のようなフォーマットになります。
出力データは18bitで取り出すようにしていますが、SINやCOSは±1の値も取りうるので、フォーマットは下の図のように整数部のある固定小数点となります。
LPFは、別の案件で使った42次のFIRフィルタを使いました。これは0.6MHz程度までを通します。
この原理でFPGAに実装して試してみました。
まず、ファンクションジェネレータで作った410.00kHzの正弦波をCosmo-Zに入れ、FFTで見てみます。
BPFに通したスペクトラムは次のようになります。直流成分や低周波のノイズなどが削除されて、目的の信号だけが残りましたね。BPFで10dBほど減りました。
ここで、ローカル発信器の周波数を416.67kHz (=5MHz/12) として、周波数シフトさせてみます。周波数シフトさせた結果が緑のスペクトラムです。
ローカル発信器と入力信号を乗算した6.67kHzの信号が原点付近にシフトされましたが、820ｋHzのほうにも高調波が出ています。
これをLPFで削ります。
すると、目的の信号だけが取り出されて残るというわけです。ただ、ちょっとLPFの効きが甘かったようで ローカル発信機の周波数あたりに残ってしまっています。
さて、これで目的の信号が周波数ゼロ付近にシフトしたのですが、まだズームはしていません。ズームFFTをするにはデシメーションという操作が必要です。デシメーションというのは、N回に1回だけのデータを残して捨てるという操作です。
デシメーションで10個に1個だけ残すようにしたら、周波数が10倍にズームされました。
410kHzの入力信号と、局発の416.67kHzの差が、ズームされて見えているのがわかります。
このようなズーム回路ですが、使用したリソースはDSPブロックが7個、SLICEが665、BRAMが0個と、やっている演算の割にとてもコンパクトでした。
N分の1にデシメーションをかけてからFFTを行うことで、もともとのN倍の長さのFFTを行っているのと同じことになり、周波数分解能をN倍に上げることができます。ズームFFTの本質は「データを捨てる」というデシメーション操作にあるので、BPFやLPFはノイズを減らすという目的のためのものにすぎません。
ズームFFTの欠点は、必要なデータ数を集めるのにかかる時間が長くなることです。10倍にズームするならば、FFT自体の計算量は変わりませんが、必要なデータを集めるのに10倍の時間がかかってしまいます。
FPGAによるデジタル信号処理って、本当に面白いですね。
Cosmo-Zのアルミケースに新型が登場
2017.04.26
Cosmo-Zのアルミケースが新しくなりました。
構造自体は変わらないのですが、フロントパネルとリアパネルにシルク印刷とアルマイト加工が施されるようになりました。
TOKUDENの字が少し上によってしまったのと、もう少し大きくすればよかったかなと思います。
シルクで印刷されているので、何のコネクタの穴なのかがすぐにわかります。
中身はこんな感じ。
汎用計測器としての高級感を醸し出せるように今後も改良を続けていきたいと思います。
ハードウェアFFTの精度向上2
2017.04.09
ハードウェアFFTを行うFPGAで、丸めの計算を単純な切り捨てから、四捨五入(0捨1入)に変更してみました。本当は最近接偶数丸めが良いのですが、0.5ちょうどになることはあまりないので、とりあえずは0捨1入で行います。
まず、従来の単純切捨てのFFTの場合です。
周波数の全域にわたってFixed Patternなノイズが見えています。これを積算して1LSB以下の構造を見てみると、
このようなノイズの構造が見えてきます。
さて、いよいよ丸めを実装した場合です。下の図のようにFixed Patternなノイズが大幅に減りました。
ほとんどの場所でスペクトラムの値は0になっているのですが、それでも4か所ほど、Fixed Pattern Noiseが出てしまっています。
積算すると、1LSB以下の構造が見えてきて・・・
やはり、ノイズが出やすい周波数があるようです。
単純切捨てに比べるとノイズは1万分の1くらいに減っているとは思いますが、それでもゼロにはなりませんでした。
ほとんどの周波数でFixed Pattern Noiseは1LSB以下に抑圧されましたが、積算すると見えてくるので、ゼロにはできなかったということです。
ハードウェアFFTの精度向上
2017.04.06
FPGAを使ってハードウェアでFFTする回路を作ったのですが、ノイズに悩まされています。
そのノイズというのは、決まった周波数にいつもノイズがいるというものです。単一の正弦波(10kHz)をFFTすると下の図のような結果になります。
１本立っているのが真のスペクトラムなのですが、全域にわたってフロアノイズが出ています。なお、16bitの演算精度なので、結果が1LSBなら-90dB、2LSBなら-84dB程度のところになります。
つまり、スペクトラムの全域にわたって±2LSB程度の誤差が出ているというわけです。その中でも特に1.2MHzと0.6MHzの値に特徴的な盛り上がりがあります。
このFFT結果を何度も積算して平均化してみると、ノイズの構造が見えてきます。
このように周期的なFixed Patternなノイズが見えてきます。
ノイズが生じる原因は、FFTの演算精度を16bitの固定小数点(つまり整数)で行っていて、バタフライ演算の後で生じる除算をした際にLSBを単純に切り捨ててしまっているから、その誤差が積み重なったのだと考えられます。
そこで、FPGAではなく、パソコンのソフトウェアを作ってFFT計算誤差を調べてみました。
最初は1024ポイントのFFTで、浮動小数点で計算します。入力データは正弦波なのですが、周波数が整数ではない場所に設定されていて、サイドローブが広がるようにしています。
f=20あたりにピークがあって、あとは滑らかに減衰するだけです。
16bitの固定小数点で、誤差を単純に切り捨てた場合の結果を以下に示します。まさにFPGAでの計算のときのようなギザギザが生じてしまっています。このような形のギザギザは演算精度の問題であることがわかります。
そこで、最近接偶数丸めというのを使ってみました。
浮動小数点計算ほど滑らかではありませんが、特徴的なギザギザはなくなったようです。
FPGAのハードウェアでFFTを行う際にも、誤差は最近接丸めで処理すれば、もっと良い結果が出るかもしれないという知見が得られました。
Cosmo-Z 18bit版の性能評価
2017.04.05
アルバイトさんが、Cosmo-Zの18bi版拡張ADCボードの特性を測ってくれました。
Cosmo-Zの18bit拡張ボードというのは、下の写真にあるようなボードで、±1Vpp入力で5MHzサンプリングのボードです。
AD変換器にはAnalogDevicesのAD7960を使い、初段のプリアンプ（可変ゲインアンプ）にはAD8253を、２段目のプリアンプ（ADCドライブ用）にはTIのTHS4521を使っています。
ノーマルのCosmo-Zに比べると速度は遅いですが、ADCの性能やチャネル間のアイソレーションは抜群に良くなっています。また、入力にはx1 x10 x100 x1000の可変ゲインアンプが入っています。ゲイン1の設定の場合の分解能は7.6μVですが、ゲイン1000にすると1LSB=7.6nVとなります。
さて、性能を見てみましょう。まずはゲイン1の場合のノイズのヒストグラムです。
このADCは±4.096Vの範囲を18bitで変換するので、1LSB=4.096V*2/262144=31.25μVとなります。上のヒストグラムは半値全幅が7LSBなので約220μVに相当するのですが、可変ゲインアンプの後にゲインは4のプリアンプを乗せているので、入力に換算するとノイズの幅は55μVとなります。
次はゲイン10の場合です。
ノイズによる広がりは22LSBになりました。入力換算ノイズは17μVとなります。
ゲイン100の場合、広がりは96LSB程度になりましたが、プリアンプのゲインがトータルで400あるので入力に換算すると7.5μVとなります。
最後はゲイン1000の場合。広がりは320LSBほどで、入力に換算したノイズは2.5μVとなります。
ゲインを10倍、100倍、1000倍と変えていっても、ノイズの増え方は10倍、100倍、1000倍とはならず、それよりゆっくり増加するのは、ノイズが初段の可変ゲインアンプだけで生じているからではなく2段目のプリアンプやADC自身で生じているためです。そのため、初段のアンプでゲインを稼ぐことは有意義といえます。
また、オーディオアナライザを使って綺麗な正弦波を作り、このボードに入力し、FFTを行ってスペクトラムを見てみるとひずみ率は-80dB程度でした。
周波数特性は約1MHzで3dBダウンとなりました。
この周波数特性を決めているのは、２段目のプリアンプで、現在は帰還抵抗が5kΩ//33pFという構成なので妥当な結果と言えます。
なお、コンデンサを減らすと周波数特性は良くなりますが、熱雑音によるノイズも増えるし、そもそも2.5MHz以上はエイリアシングだから無い方が良いので、難しいところです。
この結果が妥当なものか検討してみます。
AD8253はゲインによってノイズの大きさが大きく変わります。
ゲイン1の場合、AD8253の出すノイズは45nV/√Hzです。5MHzくらいまで見ているとすると、AD8253で225μVのノイズが発生していることになります。また2段目のTHS4521は4.6nV/√Hzなので23μV、帰還抵抗は5kΩなので熱雑音は20μV。二乗して全部足し合わせると、227μVとなります。
大ざっぱな計算ですが、最初のヒストグラムで示した220μVというのとほぼ一致しますので、ゲイン1のノイズの99%の発生原はAD8253であると推定されます。
ゲイン10の場合、AD8253の出すノイズは60μVに減りますが、これがゲインで10倍されるので600μVのノイズが2段目のTHS4521に加わることになります。THS4521の出すノイズと熱雑音は変わらないので、観測されるトータルで600μVとなります。ヒストグラムから読み取ったノイズの幅は22LSB=687μVなので、これも結果とよく一致します。
ゲイン100の場合、AD8253の出すノイズは55μVに減りますが、ゲインで100倍されるので5500μVのノイズが2段目のTHS4521に加わることになります。ヒストグラムの広がりは約100LSB=3.1mVなので、実際は計算値ほど悪くないという結果でした。
ゲイン1000の場合、AD8253の出すノイズは50μVに減りますが、ゲインで1000倍されるので5mVのノイズが2段目のTHS4521に加わることになります。ヒストグラムの広がりは約300LSB=10mVなので、実際は計算値ほど悪くないという結果でした。
ゲイン100以上で、ノイズが想定より少なくなるのはAD8253の周波数特性がG=100,G=1000の場合には早い段階で下がり始めるからでしょう。
つまり、ノイズ発生源であるAD8253は、G=100以上で周波数特性が悪くなるため、高い周波数成分のノイズも減ったと考えられます。
固定小数点回路と浮動小数点回路
2017.04.01
XILINXのFPGAである種の信号処理回路を作っているのですが、固定小数点回路と浮動小数点回路で結果が一致しました。
固定小数点では、16bit表現のとき、
符号.有効数字15桁
のようにして-1～1を表現する場合、結果を32768で割ることになります。
結果がパワーのような「二乗」で表現されている場合、16bitでは表現しきれないので32bit（二乗の和を積和する場合は48bit程度に拡張して）計算して、最終結果を32768*32768で割ることになります。
このような桁の換算が結構混乱しますが、演算自体は高速です。また、適当なビット([34:19]とか)をスライスして切り出せば、結果を目ですぐに見ることができ便利です。
一方、浮動小数点の場合は桁の換算は不要なのですが、XILINXのIPコアはIEEE754という形式で数を扱うので、これについて理解しておかなければなりません。詳しくはWikipediaを参照してください。
私は下記の形式にします。
IEEE754で表現された数値は、16進表現をぱっと見て人間が理解できるものでもないので、デバッグはやりにくいです。
C言語のfloatと互換性があるのかどうかはよくわかりませんが、IEEE754の32bit形式で格納された値をfloat型に変換するには、
*(float *)(&val)
で、要するに無理やりキャストでうまくいくようでした。
固定小数点と、浮動小数点の両方で計算したある種の計算の結果が、先のグラフに示したとおりで、まったく同じになりました。