「C++ Composer XE 2011」ベンチマーク記事の隙が大きすぎて、アーキテクチャ批評が書けない

Sourceforge Magazineのデベロッパー向け記事にアセンブラコードで見るC++ Composer XEの強力な最適化機能がリストされています。さて、その記事、SSEの性能評価をしてやろうかと思ったらそれ以前の問題で膝を折りたい気分です。
一例は2ページ目

シンプルなベクトルの内積計算を行わせて、その生成コードをIntelコンパイラとVS2008のコンパイラで比較しようとしています。
プログラムはこんな感じ（リスト1より引用）

int i;
float a[128];
float b[128];
float d;
  ：
  ：
d = 0.0;
for (i = 0; i < 128; i++) {
  d += a[i] * b[i];
}
printf("%f\n", d);
  ：
  ：

このコードからIntelコンパイラは4並列浮動小数点演算を行うSSE命令をびっちりと出力してその実力を見せつけています。DSP屋としてはSSEに並列ロードストア命令がないのが不満ですがIntelやAMDのx86はスーパースケーラーの力業で命令を並列化しますから、そこはあまり問題ではありません（そういう事にしておきます）。
一方のVS2008コンパイラはFPUを使ったスタックコードをこれもびっちりとはき出します。
勝負あった、と思いきや、100万回実行した結果はこんな感じです（表4より引用）。

コンパイラ 実行時間
インテル C++ Composer XE 2011 9.325秒
Visual Studio 2008 9.722秒

ええええ？って、1MACに170サイクル以上かけていますから、明らかにprintf()の時間まで含んでいます。「ここで使用しているコードはシンプルであるため軽微な違いしか見られないものの、SSEを利用することによるパフォーマンス向上が確認できる」と記事にはありますが、少々すっとぼけた過ぎていないでしょうか。何のためにアセンブリ出力まで引っ張り出して比較しているのかと問いたださずにはいられません。
うがった見方をすれば、フォン・ノイマン・ボトルネックのために思ったほど差が出ない事を隠蔽するためにprintf()を取り除かなかったとも考えられます。
MAC演算は信号処理の基本演算であり、SSEがどの程度の性能向上をもたらすか興味があるのですが、残念ながらそれをえぐり出す鋭さはこの記事にはありませんでした。
Intelコンパイラはお試し版があるそうですが、わざわざダウンロードして調べるほど興味があるわけじゃないんですよね。

Real FFT

FFT

FFTは複素信号を処理します。入力は複素信号であり、出力も複素信号です。ところが、実数しか必要としないアプリケーションはたくさんあります。そのような場合、虚部を0としてFFTを行うわけですが、もったいない点は否めません。何とかならないかと思うのが人情です。そしてうまい具合に何とかする手だてがあります。Real FFTです。

実数データの複素FFT結果

はじめに、実数データr[]を複素FFTにかけた場合の結果について考えます。

R=FFT(r)

ここでFFT結果の実部Re(R)と虚部Im(R)を調べてみると、おもしろいことがわかります。

  • Re(R[])は左右対称になる（偶関数）
  • Im(R[])は点対称になる（奇関数）

ちょっとScilabで適当なデータを作ってFFTにかけた結果をお見せしましょう。入力データは実数配列です。
最初にお見せするのは実部 Re(R)。左右対称になっています。
[f:id:suikan:20101219195452p:image]
次が虚部 Im(R)。点対称になっています。

この結果から次のことが言えます。
まず、実部のデータですが、中点を軸にして紙を折るようにして対応する両側のデータを足しあわせると*1、結果はそれぞれのデータが２倍になります。なぜならデータは線対称になっており、同じデータが左右に広がっているからです。
一方、実部のデータですが、同じように中点を軸にして両側のデータを足しあわせると、対応するおのおののデータは絶対値が同じで符号が逆なのですべて0になってしまいます。

虚数データの複素FFT結果

つぎに、実部のない虚数だけのデータをFFTにかけるとどうなるか考えてみましょう。これはFFTの性質から簡単に導き出せます。
FFTは線形な変換です。つまり、

FFT(Ar) = A FFT(r)

が成立すると言うことです。ここから直ちに

FFT(ir) = iFFT(r) = iR[]

だということがわかります。言い換えると、純粋虚数データのFFT結果は、実数データのFFT結果にiをかけたものになります。
したがって、irのFFT結果であるiRは次のようになります。

  • Re(iR)= -Im(R)は点対称になる（奇関数）
  • Im(iR)= Re(R)は左右対称になる（偶関数）

実データの虚データのFFT結果の比較

ここまでの結果をまとめるとこうなります。

データ処理 結果の実部 結果の虚部
実データのFFT 偶関数 奇関数
虚データのFFT 奇関数 偶関数

これを利用すると、複素データのFFT結果から、実データのFFT結果と虚データのFFT結果を分離することができます。具体的には中心を軸として対応するデータの和と差をとることで、偶関数データと奇関数データを取り出すことができます。これは非常におもしろい性質です。

Real FFTの実装方法

複素データのFFT結果から実データのFFT結果と複素データのFFT結果を取り出せることから、実データFFTの計算量を以下のようにして削減できます。

  1. N点の実データを前半と後半に二分する。
  2. 後半のデータにiをかけて前半のデータと足す
  3. 合成した複素データをN/2点FFTにかける
  4. 結果から実部のFFT結果と虚部のFFT結果を取り出す
  5. 虚部のFFT結果に-iをかけて実部のFFT結果の後ろにつなげる
  6. N点バタフライをかける

これによって、N点実データのReal FFTをO(N/2 log(N/2))で実行できます。Nが十分大きければ、N点複素FFTの半分の時間で実FFTを実行できます

まとめ

FFTは幾分魔術的なアルゴリズムですが、根底にある線形性を活用するとこのように劇的に実データに対する演算量を減らすことができます。同じ性質を使ってステレオデータの左右の信号を一挙にFFTにかけるといったこともできます。

*1:ロールシャッハテストのように

Codesourcery G++ LiteでCortex-M4Fプログラムをコンパイルする

CORTEX-M4

2010年秋版のCodesourcery G++ LiteはCORTEX-M4Fに対応しています。試行錯誤した結果、以下のようにそれらしいコンパイル結果を得ることができました。

takemasa@takemasa-desktop:~/foo$ cat main.c
volatile float a[20], b[20]; 

int main(void)
{
	float sum=0.0;
	int i;

	for( i=0; i<20; i++ )
		sum += a[i]*b[i];
	return (int) sum;
}

takemasa@takemasa-desktop:~/foo$ 
takemasa@takemasa-desktop:~/foo$ arm-none-eabi-gcc -O2 -mfpu=fpv4-sp-d16 -mfloat-abi=hard  -mthumb -mcpu=cortex-m4 -S main.c
takemasa@takemasa-desktop:~/foo$ cat main.s
	.syntax unified
	.cpu cortex-m4
	.eabi_attribute 27, 3
	.eabi_attribute 28, 1
	.fpu fpv4-sp-d16
	.eabi_attribute 20, 1
	.eabi_attribute 21, 1
	.eabi_attribute 23, 3
	.eabi_attribute 24, 1
	.eabi_attribute 25, 1
	.eabi_attribute 26, 1
	.eabi_attribute 30, 2
	.eabi_attribute 18, 4
	.thumb
	.file	"main.c"
	.text
	.align	2
	.global	main
	.thumb
	.thumb_func
	.type	main, %function
main:
	@ args = 0, pretend = 0, frame = 0
	@ frame_needed = 0, uses_anonymous_args = 0
	@ link register save eliminated.
	push	{r4}
	movw	r1, #:lower16:a
	movw	r2, #:lower16:b
	flds	s15, .L6
	movs	r3, #0
	movt	r1, #:upper16:a
	movt	r2, #:upper16:b
.L2:
	ldr	r4, [r1, r3, lsl #2]
	ldr	r0, [r2, r3, lsl #2]
	adds	r3, r3, #1
	fmsr	s0, r4
	fmsr	s1, r0
	ldr	r4, [r1, r3, lsl #2]
	ldr	r0, [r2, r3, lsl #2]
	fmacs	s15, s0, s1
	fmsr	s13, r4
	fmsr	s14, r0
	adds	r3, r3, #1
	cmp	r3, #20
	fmacs	s15, s13, s14
	bne	.L2
	ftosizs	s15, s15
	fmrs	r0, s15	@ int
	pop	{r4}
	bx	lr
.L7:
	.align	2
.L6:
	.word	0
	.size	main, .-main
	.comm	a,80,4
	.comm	b,80,4
	.ident	"GCC: (Sourcery G++ Lite 2010.09-51) 4.5.1"

実機で確認していませんのでこれで決まりとは言えませんが、それらしい結果が出ています。苦労したのは-mfpu=fpv4-sp-d16 -mfloat-abi=hardを指定する点です。それでも後者はドキュメントを検索すればわかりましたが、前者の-fpu=fpv4-sp-d16には難儀しました。ほとんど隠しパラメータ扱いです。
ARMの浮動小数点コプロセッサはNEONとVFP(Vectored Floating Point)と呼ばれるコプロセッサが主流ですが、VFPは倍精度命令を持っています。CORTEX-M4はこのVFPの単精度派生品のFPv4-spで、GCCはそれをfpv4-sp-d16と呼んでいます。ああめんどくさい。
性能のほうは、ベクトルのスカラ積を実行したときに、1MACに9サイクルです。事前のおおざっぱな見積もりがほぼ的中しています。

FFT性能を推測する

CORTEX-M4

CORTEX-M3関係の資料を洗い直していたところ、NXPのニュースリリースに「[http://www.nxp.com/news/content/file_1679.html:title=120MHzのCORTEX-M3マイコンで、1024pt 16bit FFTが
0.89mS以下]」という記述を見つけました。
ここから少し推論してみます。CORTEX-M4コアのLPC4300は150Mhzで動作しますので、同じプログラムを0.71mS以下で実行できます。っここで荒っぽい推論ですが、CORTEX-M4Fコアが、32bit浮動小数点演算を16bit固定小数点演算と同じサイクル数で実行できるとします。すると、LPC4300は32bit 1024pt浮動小数点FFTを0.71mSで実行できることになります。
1024pt FFTは512ptの複素FIRフィルタの実装に使われます。48KHzサンプルの時、512サンプルはおおよそ11mSですので、CORTEX-M4Fは、512pt複素FIRを1/7強のMIPS消費量で実行することになります*1。
これはすこぶるいい値に思えます。いよいよLPC4300が楽しみです。

*1:FFTとIFFTのペアで実行するので

プロジェクトほったらかしだった orz

ちょっと人から聞かれまして、Kobanzame版のTalkthroughを引っ張り出そうとしたのですが…ありませんでしたっ！
なんてこと、記憶の引き出しを探っても出てきません。自分のPCを検索してようやくPizzaFactory版がTOPPERS/JSP for Blackfinプロジェクトにアップロードされていることを突き止めました。
イニシャライザでクロック周波数などをいじっていますので、少しプログラムを変更した後、コマンドラインからビルドできることを週末に確認する予定です。