So-net無料ブログ作成
検索選択

パフォーマンス比較の続きとFFT実装の結論 [任意点数のFFT]

一通りやろうと思っていたFFTの実装を終わらせて、それぞれの効率を比較してみた。離散Fourier変換の計算はいろいろな規則性があってそれを上手く使う賢いアルゴリズムがいっぱいあって面白い。
今日は昨日の続きでNの大きな2のベキの場合の比較をする。

ひとつ、驚愕の結果が得られた。実は、いまだに信じられない。なにかの間違いじゃないか、という気がしてしょうがない。

2のベキの場合

Nが2のベキの場合を比較してみたのが図-6。
0912fig6.png
こっちには
  1. 定義式通り
  2. 混合基数
  3. Chirp z変換
  4. fftw
  5. 2の基数のFFT(独自実装)
  6. vDSPのvDSP_fft_zip()関数を使ったもの
の6つを並べた。定義式通りはNの大きな値では何桁も遅くなって測定が大変だった。せっかく測ったのにlog-logプロットでも離れ過ぎてグラフが見づらくなるので、外にはみ出させて描かないことにした。

どれもほぼ直線にのっているが、定義式通りの線だけが傾きが違う。これだけ傾きが2で、つまりN2に比例しているということ。それ以外はどれもN1.2前後でK-T型のアルゴリズムが効いていることがわかる。

独自実装の2の基数のFFTはあんがいがんばってfftwとそう遠くないスピードが出ている。といっても倍ほど違うけど。混合基数のクラスは2のベキの場合には2の基数のFFTと式の上では一致するけど、パフォーマンスは6、7倍違っている。かけ算回数は同じだけどメソッド呼び出しが増えているのでそこで最適化が効かなくなっていると思われる。

意外だったのはfftwよりもvDSPの方が断然速い!ということ。10進ヒト桁違う。これは驚き。なんか間違ったんじゃないか、と思って計算結果をfftwのものと比較してみたりした。

vDSPもfftwも今回使ったのはfloat幅のものだけなので計算精度はちょっとおちるけど、N=1024でどちらも十進で5桁の精度は出ていることを確認した(ちなみに僕のコードの精度も同程度だった)。

ほかにもvDSPがなにかズルしてるのではないかと思ってあら探ししたが、これ、というのは見つからなかった。

PowerPCの時代のベンチマークではfloat幅の場合にfftwよりvDSPの方が速いという結果はあったが、それでもここまで差は出ていない。

ひとつには、今回使ったfftw3.2.2はずいぶん昔にコンパイルしたものでgccが3.*時代だった。それ以外はgcc4.2でコンパイルしている。こないだ見たようにgcc4.0から自動ベクトル化ができるようになって、fftwのmakeでベクタユニット使用の指定をしたときとは違った使い方をする、と言うのことがあるかもしれない。しかしそれにしても一ケタ差がでるとは思えない。

僕の以前のメインマシンであるPowerBook G4(もちろん今でも現役)でコンパイルして比較するとfftwとvDSPはとんとんの結果を出すのでIntelのCPUに対して固有の最適化をしているんだろう。

屋根屋のふんどし泣き別れのメモリ配置がそんなに効くんだろうか?

アルゴリズムはまったく同じはずなので、メモリ配置とかベクタユニットの使いこなしとかで、ここまで差が出るとは思えないんだけどなあ。ひょっとしたらGPUとかまで使ってたりするのかなあ。どう考えても不思議だ。実はカーネル空間で実行してたりして(今回ユーザ時間で測っていた。あとで実時間で測っても同じだった)。

いや、目からウロコでした。vDSPはフェイドアウトかなんて言ったけど、なかなかどうして。この実力ならそれに見合ったドキュメントを整備してほしいもんだとは思うけど。

結論

任意の点数に対応するFFTのアルゴリズムを一通りおさらいしてそれを実装した。その結果目標とするパフォーマンスには十分届くものができた。

MacOS Xの場合、最終的にはやはりvDSPの2の基数のFFTを利用するのがパフォーマンス上は望ましい。そうするとChirp z変換の実装にも使うべきである。

また、混合基数の今回の実装では、与えられたNを小さな素因数から分解したけど、大きい方から分解して、残りが2のベキになる場合にはvDSPを使う、というのが賢いやり方だろう。

また今回全然タッチしなかったRaderのアルゴリズムを勉強してもいいかもしれない。整数の性質は不思議だ。

FFTも奥が深い。

nice!(0)  コメント(11)  トラックバック(0) 

nice! 0

コメント 11

zuruyasumineko

チャープZ変換の任意サイズFFTは、DFTと同じ結果がでるのでしょうか?
だとしたらすごいですね。もしDFTと結果が一致するならDCTで実装してみて、高速リサイズができるかもしれません。
by zuruyasumineko (2013-10-07 04:05) 

decafish

コメントありがとうございます。
Chirp-Zは2のベキのFFTなどと同じで式の上では離散Fourier変換と一致します。ですので原理的には誤差が含まれるということはありません。ただ、floatやdoubleの有限精度で計算すると丸めや桁落ちなどによって計算誤差がでます。Chirp-Zは計算量が比較的多い(最大で倍の長さ-1の2のベキのFFTを2回)のでその累積誤差も多くなる可能性があります。
「高速リサイズ」というのはどう言うお話なのでしょうか?面白そうなお話でしたらぜひ伺いたいのですが....
by decafish (2013-10-07 22:01) 

zuruyasumineko

ご返事ありがとうございます。6年ほど前から、画像サイズのリサイズを考えていまして、リサイズは通常FIRフィルターを使って空間領域でやるのですが、ナイキスト周波数以下に帯域制限されてしまう画像の拡大縮小において、理想的な遮断特性で、エイリアス成分を除去できるのが、DCTを使ったリサイズです。FFTを使うには、データを鏡像対称に2倍にしてFFTして、実数部だけとりだすという作業が追加になりますが、データ数が素数でも、チャープZ変換なら、高速化できると思いますので、やりがいがあると思います。私のホームページとブログをお知らせしておきます。今回のお返事をきいてチャープZ変換によるFFTで任意サイズDCTリサイズを実装したいと思いました。

ブログ http://blogs.yahoo.co.jp/zuruyasumineko2002

by zuruyasumineko (2013-10-08 00:54) 

zuruyasumineko

実は、チャープZ変換には前から気づいていたのですが、0埋めがあるので
DFTとは違う結果がでるのではないかと思い込んでいました。
例えば、100データをDFTするのに、データ数128のFFTを使って、28個0埋めするということだと、データ数100のDFTとは結果が違ってきます。
なので、そうではないとすると、チャープZ変換によるFFTが、一般によく知られるようになってもいいのになあと思いました。
by zuruyasumineko (2013-10-08 01:09) 

decafish

なるほど。FIRフィルタでは、拡大では本来存在しない周波数成分が現れたり、縮小ではなくなるはずの成分に結果が影響されたりするが、DFTを使えばそういうことがなく、エイリアスみたいな目立つ現象も起こらないということでしょうか。

Chirp-Zではたしかに0を埋めますが、そのままFFTするのではなく、畳み込みにする(2回繰り返す)ことで、うまいぐあいに埋めた0によって不要な成分が消えてくれる、と私は理解しています。詳細はもう忘れてしまいましたが、初めて勉強したとき「これを考えた人はなんて頭がいいんだ!」と感心したことは覚えています。
ぜひChirp-Zを試してみて下さい。

ところでDFTを使って画像を拡大すると、見た目にはわからないように別の情報が埋め込めて、元に戻したとき画像劣化無しに取り出せる、ということになるのでしょうか? まあ、できてもあまり面白くはないかも知れませんが。
by decafish (2013-10-08 12:52) 

zuruyasumineko

chirp-zの方は、DFTとあってるかの検証含めて、気長にやってみます。

別の情報が埋め込めるかについては、
拡大の際に、現画像の最大周波数成分と拡大画像のナイキスト成分の間に、いれることは可能だと思います。みためにわからなくするために、ホワイトノイズ的に入れることになると思います。また、画素値は、0~255の値をとりますから、絶対値が1の信号をいれてやれば、あまりわからないかもしれません。
by zuruyasumineko (2013-10-08 17:33) 

decafish

よく考えればDFTは線形なのであたりまえですね。広がった周波数領域の部分に、ホワイトノイズに見えるように分布させるためのアルゴリズムのほうが難しそうです。一種のスペクトラム拡散みたいなのを考えればいい、ということでしょうけど。
by decafish (2013-10-08 21:08) 

zuruyasumineko

その任意サイズDCTリサイズをチャープZ変換のFFTで作りました。
FFTで畳み込みを計算することは、今回はじめてやりました。
0埋めしても、誤差がでないことがわかりました。
真っ正直なDFTでやると計算に2,3日かかるところが、数分でおわるようになりました。ありがとうございました。
by zuruyasumineko (2014-08-02 19:47) 

decafish

コメントありがとうございます。
お役に立てたのなら、うれしいです。それほど高速化されるということはかなり点数が多いのですね。Cooley-Tukey型のFFTがいかに強力か、という見本のようです。

Chirp-Z変換によるFFTのアイデアは膝ぽんものですごく面白いと思います。考えだした人はきっと頭の中でああでもないこうでもない、とこねくり回していて、あるときぱっ、と気がついたんだろうと思います。
僕もそういう体験をしてみたいです(もう遅いですけど)。
by decafish (2014-08-03 20:35) 

zuruyasumineko

>FIRフィルタでは拡大では本来存在しない周波数成分が現れたり、縮小ではなくなるはずの成分に結果が影響されたりするが、DFTを使えばそういうことがなく、エイリアスみたいな目立つ現象も起こらないということでしょうか。

DFTだと、直流分からサンプリング周波数の1/2まで、周波数特性を自由につくれるので、遮断特性が、急峻なものもつくれると思ってました。
縮小ではエイリアスがないものもつくれると。
それでDFTで画像を拡大縮小してもいいんですが、実偶関数に対するDFTの半分の結果を使うDCTというのが、DFTより良い結果が得られるようです。それでも、画像全体にDCTをかける拡大縮小法ですと、人口画像(アニメ画、CG画、文字画像など)ではギブス現象で、画像はみにくくなります。風景写真では、高周波成分が少ないので、きれいな画像になりますが。といったことが実験でわかっています。最近では、JPEG、MPEGで使われる、8x8マクロブロック的に、細かく分割する手法をちょっと工夫して、つなげたときに境界がでないようにしてやるという大き目のマクロブロックで変換して、まわりはすててつなげるという手法で、ギブス現象の尾引きの長さを縮めるというとこまでいってます。lanczos窓をかけたsinc関数による拡大縮小よりもシャープな画像が得られるまでになりました。
by zuruyasumineko (2017-08-05 23:57) 

decafish

お久しぶりです。
DCTではDFTと違って幅が倍になって必ず偶関数になるので高周波成分が抑制されるのだととりあえず理解しているのですが、線形な変換同士で結果が違うということが不思議で、やっぱり僕はちゃんとは理解できていないようです。
JPEGのブロックノイズは煩わしいのでそれが抑えられる手法があればそれは面白そうです。
lanczos窓よりシャープというのはすごいですね。lanczos窓は拡大した時Fourier領域で考えると高周波成分をごにょごにょする(無い成分を作る)ことになると思うので、絵のシャープさという意味では結構いいんだと思っていましたが。
by decafish (2017-08-07 08:51) 

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

トラックバック 0

メッセージを送る