nabeの雑記帳

長年理解するのを拒否していたFFTとここ何週間か格闘しています。楽々とアルゴリズムを導出・実装できる人はいいのでしょうが、FFTのようなアルゴリズムをきちんと理解して実装するのは思うほど容易ではありません。既に導出されたアルゴリズム（や雛形サンプル）を何も考えず実装すれば簡単ですが、そういうことが生理的にできない場合、改良とかを考え出して簡単には実装できなくなるんですよね。

前置きはともかく、今日やっと実数専用のFFT/IFFTが作成できました。実数専用にすることで理論的には通常のFFTの半分の時間で処理できます。

• rfft.cpp

見ての通り、実数FFTの対称性を利用したものすごく面倒くさく複雑な作りになっているんですが、有名なFFTの概略と設計法のソースと速度比較をしたら倍近く遅い……。64K点、double型、VIA C3 500MHzで動作させて、自作が100ms、リンク先のソースが60ms。

いかに細かい最適化をしても実装アルゴリズムの時点で差があると非常に大きいですね……。あとでアルゴリズムの差を検証予定。

ソース解読メモ

大浦氏のFFTソース(fft4g.c)解読メモ。つまりリバースエンジニアリングです。

サブルーチンの機能

IFFTをFFTで代用する方法

なんだけども、複素共役を考えると

つまり前処理と後処理としてデータの複素共役を取ってあげれば、同じFFTルーチンを使い回すことができる（効率を考えて、初段部分と終段部分だけ個別特別に実装すれば、中段部分は使い回せる）。

その他メモ

• 一般的なFFT(DFT)と比べ虚数部の符号が逆になっている。*1

周波数間引きFFTである

やっとこのプログラムの要が理解出来ました。

先にデータをスクランブルして（並べ替えて）いるため「時間間引きFFT」に見えますが、実際には周波数間引きFFTです。sin/cosテーブルもわざわざスクランブル位置に格納しています。

こうすることでFFT実行時のデータアクセスをシーケンシャルに行え*2、それが功を奏して実行速度が飛躍的に速くなっています。またsin/cosテーブル、偶数番地に cos、奇数番地に sin を格納し、全体としても π/2 しか用意しないことで、データアクセス量を減らしキャッシュが効きやすくなっています。

256点の複素FFTを実行させたときの、cos/sinテーブル参照位置およびデータ参照位置は次のようになります。

 call cftmdl(512, 8)
  k1=1 (0.923880,0.382683), k2=2 (0.707107,0.707107)
    32 36 40 44
    33 37 41 45
    34 38 42 46
    35 39 43 47
  k1=2 (0.980785,0.195090), k2=4 (0.923880,0.382683)
    64 68 72 76
    65 69 73 77
    66 70 74 78
    67 71 75 79
  k1=3 (0.831470,0.555570), k2=6 (0.382683,0.923880)
    96 100 104 108
    97 101 105 109
    98 102 106 110
    99 103 107 111
  k1=4 (0.995185,0.098017), k2=8 (0.980785,0.195090)
    128 132 136 140
    129 133 137 141
    130 134 138 142
    131 135 139 143
  k1=5 (0.881921,0.471397), k2=10 (0.555570,0.831470)
    160 164 168 172
    161 165 169 173
    162 166 170 174
    163 167 171 175
  k1=6 (0.956940,0.290285), k2=12 (0.831470,0.555570)
    192 196 200 204
    193 197 201 205
    194 198 202 206
    195 199 203 207
  k1=7 (0.773010,0.634393), k2=14 (0.195090,0.980785)
    224 228 232 236
    225 229 233 237
    226 230 234 238
    227 231 235 239
 call cftmdl(512, 32)
  k1=1 (0.923880,0.382683), k2=2 (0.707107,0.707107)
    128 144 160 176
    129 145 161 177
（中略）
    142 158 174 190
    143 159 175 191

ほんと、この実装はすごいなぁ。

電池動作機器を設計するのに、電池の情報がなかなか纏まっておらず集めてみました。実際に測定したわけではなく、ネット上の情報を適当にかき集めただけですので、あくまで参考程度でお願いします。

より正しい値や、他の情報がありましたらお知らせください。

int型は整数型ですが、int や long や long long などは環境によってサイズが違ったりします。これらの型は、元もとサイズ（バイト長）を固定する目的で作られたものではないからです。

プラットホーム汎用でプログラムを書く際、固定長のデータを扱う時（バイナリデータ列など）は特に注意する必要があります。C99という規格で固定長整数型として int32t などが規定されましたが、すべての環境で使えるわけではありません。また Windows プラットホームならば、__int32 などが利用出来ますが、Windows以外では利用出来ません。

自分は、だいたいの環境でうまく動くマクロを作って、これを使っています。

#if defined(__C99__) || (defined(__GNUC__) && __GNUC__ >= 3)
#	include <inttypes.h>
#	include <stdint.h>
#else
#	if defined(__GNUC__)
		typedef short			 int16_t;
		typedef unsigned short		uint16_t;
		typedef int			 int32_t;
		typedef unsigned int		uint32_t;
		typedef long long		 int64_t;
		typedef unsigned long long	uint64_t;
#	elif (_MSC_VER || __BORLANDC__)
		typedef __int16			 int16_t;
		typedef unsigned __int16	uint16_t;
		typedef __int32			 int32_t;
		typedef unsigned __int32	uint32_t;
		typedef __int64			 int64_t;
		typedef unsigned __int64	uint64_t;
#	endif
#endif

ハミング符号(15,11)の生成行列。

生成多項式

生成行列の元データ

0001 : crc=03h 0011
0002 : crc=06h 0110
0004 : crc=0Ch 1100
0008 : crc=0Bh 1011
0010 : crc=05h 0101
0020 : crc=0Ah 1010
0040 : crc=07h 0111
0080 : crc=0Eh 1110
0100 : crc=0Fh 1111
0200 : crc=0Dh 1101
0400 : crc=09h 1001

生成行列（先頭11行省略、左がbit0、右がbit10）

p1 = [10011010111]
p2 = [11010111100]
p3 = [01101011110]
p4 = [00110101111]

後学のために買いたい書籍

誤り訂正符号入門