離散フーリエ変換と高速フーリエ変換

[みけCAT]

チャットでh2so5さんに離散フーリエ変換のソースコードを教えてもらいました。
これはそれをC言語に直し、少し調整したものです。
でもDXライブラリを使うのでC++としてコンパイルしています。

コード:

#include "main.h"

/*離散フーリエ変換*/
void calcSpectrum(const int *signal, int length, long long *spectrum) {
	int i,j;

	int N = length;

	for (i = 0; i < N; i++) {

		double a = 0, b = 0;

		for (j = 0; j < N; j++) {
			int m = j * i;

			a += cos( 2 * PHI * m / N ) * signal[j];
			b += -sin( 2 * PHI * m / N ) * signal[j];
		}

		spectrum[i] = (long long)(a * a + b * b);
	}
}

これをもとに、DXライブラリでwavファイルを再生しながら再生位置のスペクトラムを表示するソフトを作りました。
このソフトを実行すると、自分の環境ではだいたい8～12fps位で動作しました。

そこで、
http://www31.atwiki.jp/abwiki/pages/261.html
を参考に（というかほぼ丸パクリ）高速フーリエ変換の関数を作ろうと思いました。
書いたものがこれです。

コード:

#include "main.h"

/*
ほぼ 
http://www31.atwiki.jp/abwiki/pages/261.html
のパクリ
*/

/*
関数 fft()の下請けとして三角関数表を作る.
*/

static void make_sintbl(int n,double* sintbl) {
	int i=0, n2=0, n4=0, n8=0;
	double c=0, s=0, dc=0, ds=0, t=0;

	n2 = n / 2 ; n4 = n / 4 ; n8 = n / 8;
	t = sin(PHI / n);
	dc = 2 * t * t ; ds = sqrt(dc * (2 - dc));
	t = 2 * dc ; c = sintbl[n4] = 1 ; sintbl[0] = 0 ; s = sintbl[0];
	for(i = 0 ; i<=n4-1 ; i++) {
		sintbl[n2 - i] = sintbl[i];
		sintbl[i + n2] = - sintbl[i];
	}
	
	if(n8 != 0) sintbl[n8] = sqrt(0.5);
	/*
	For i = 0 To i< n4-1
	sintbl[n2 - i] = sintbl[i]
	Next
	*/
	for(i = 0 ; i<=n2 + n4-1 ; i++) {
		sintbl[i + n2] = - sintbl[i];
	}
}

/*
関数 fft()の下請けとしてビット反転表を作る.
*/
static void make_bitrev(int n,int* bitrev) {
	int i=0, j=0, k=0,n2=0;

	n2 = n / 2 ; i = 0 ; j = 0;
	while(1) {
		bitrev[i] = j;
		i++;
		if(i >= n) break;
		k = n2;
		while(k <= j) {
			j = j - k ; k = k / 2;
		}
		j = j + k;
	}
}
/*
高速Fourier変換 (Cooley--Tukeyのアルゴリズム).
標本点の数 n は2の整数乗に限る.
x[$k$] が実部, y[$k$] が虚部 ($k = 0$, $1$, $2$,
\ldots, $| n| - 1$).
結果は x[], y[] に上書きされる.
$ n = 0$ なら表のメモリを解放する.
$ n < 0$ なら逆変換を行う.
前回と異なる $| n|$ の値で呼び出すと,
三角関数とビット反転の表を作るために多少余分に時間がかかる.
この表のための記憶領域獲得に失敗すると1を返す (正常終了時
の戻り値は0).
これらの表の記憶領域を解放するには $ n = 0$ として
呼び出す (このときは x[], y[] の値は変わらない).
*/
static int last_n;  /* 前回呼出し時の n */
static int* bitrev; /* ビット反転表 */
static double* sintbl; /* 三角関数表 */

int fft(int n,double* x,double* y) {
	int i=0, j=0, k=0, ik=0, h=0, d=0;
	int k2=0, n4=0, inverse=0;
	double t=0, s=0, c=0, dx=0, dy=0;
	
	/* 準備 */
	if(n < 0) {
		n = -n ; inverse = 1; /* 逆変換 */
	} else {
		inverse = 0;
	}
	n4 = n / 4;
	if(n != last_n || n == 0) {
		last_n = n;
		if(sintbl != NULL) HeapFree(GetProcessHeap(),0,sintbl);
		if(bitrev != NULL) HeapFree(GetProcessHeap(),0,bitrev);
		if(n == 0)return 0; /* 記憶領域を解放した */
		sintbl = (double*)HeapAlloc(GetProcessHeap(),0,(n + n4) * sizeof(double));
		bitrev = (int*)HeapAlloc(GetProcessHeap(),0,n * sizeof(int));
		if(sintbl == NULL || bitrev == NULL) {
			return 1;
		}
		make_sintbl(n, sintbl);
		make_bitrev(n, bitrev);
	}
	for(i = 0 ; i<=n-1 ; i++) { /* ビット反転 */
		j = bitrev[i];
		if(i < j) {
			t = x[i] ; x[i] = x[j] ; x[j] = t;
			t = y[i] ; y[i] = y[j] ; y[j] = t;
		}
	}
	k = 1;
	while(k < n) { /* 変換 */
		h = 0 ; k2 = k + k ; d = n / k2;
		for(j = 0 ; j<=k-1 ; j++) {
			c = sintbl[h + n4];
			if(inverse != 0) {
				s = - sintbl[h];
			} else {
				s = sintbl[h];
			}
			i = j;
			while(i > n) {
			
				ik = i + k;
				dx = s * y[ik] + c * x[ik];
				dy = c * y[ik] - s * x[ik];
				x[ik] = x[i] - dx ; x[i] = x[i] + dx;
				y[ik] = y[i] - dy ; y[i] = y[i] + dy;
				i = i + k2;
			}
			h = h + d;
		}
		k = k2;
	}
	if(inverse == 0) { /* 逆変換でないならnで割る */
		//For i = 0 To i = n-1
		for(i=0;i<=n-1;i++) {
			x[i] = x[i] / n ; y[i] = y[i] / n;
		}
	}
	return 0;/* 正常終了 */
}

/*nが2のn乗であるかを判定*/
static int n_check(int n) {
	if(n==0)return 1;
	while(n%2==0)n/=2;
	return n==1?1:0;
}

/*この関数を使用する*/
void calcSpectrum_fft(const int *signal, int length, long long *spectrum) {
	double* x;
	double* y;
	int i;
	if(length==0) {
		fft(0,NULL,NULL);
		return;
	}
	/*不正*/
	if(!n_check(length)) {
		for(i=0;i<length;i++)spectrum[i]=0;
		return;
	}
	x=(double*)HeapAlloc(GetProcessHeap(),0,length*sizeof(double));
	y=(double*)HeapAlloc(GetProcessHeap(),0,length*sizeof(double));
	if(x==NULL || y==NULL) {
		if(x)HeapFree(GetProcessHeap(),0,x);
		if(y)HeapFree(GetProcessHeap(),0,y);
		return;
	}
	for(i=0;i<length;i++) {
		x[i]=signal[i];
		y[i]=0;
	}
	fft(length,x,y);
	for(i=0;i<length;i++) {
		spectrum[i]=(long long)(x[i] * x[i] + y[i] * y[i]);
	}
	HeapFree(GetProcessHeap(),0,x);
	HeapFree(GetProcessHeap(),0,y);
}

この関数を使うと、先程の離散フーリエ変換とはかなり違ったスペクトラムが表示されました。
まず、スペクトラムの値が非常に小さいです。
これは、理想的には離散フーリエ変換と同じスペクトラムになるものなのですか？
それとも、まったく別物なのですか？
理想的には同じになるのでしたら、どこを直せばいいか教えていただければありがたいです。
よろしくお願いします。
プロジェクト一式をうpしておきます。

[softya（ソフト屋)]

サイン波を入れて見ましたが、離散フーリエ変換は２つの山ができます。これも変です。
高速フーリエ変換では山さえ出来ません。窓関数に問題有るような気がします。

この１００～８００Hzデータを使用。
「サイン波wavファイル | 芽我拒グループ web支店」
http://mkgweb.undo.jp/local/advan/signwav

[みけCAT]

＞サイン波を入れて見ましたが、離散フーリエ変換は２つの山ができます。これも変です。
h2so5さんは結果に対称性があるからいいみたいな事を言っていた気がするのですが…どうなのでしょうか？

＞高速フーリエ変換では山さえ出来ません。窓関数に問題有るような気がします。
だから質問しているのです。窓関数とは何ですか？

[softya（ソフト屋)]

>h2so5さんは結果に対称性があるからいいみたいな事を言っていた気がするのですが…どうなのでしょうか？
私も遊びでフーリエ変換してみたぐらいのことしか無いですし、高速フーリエ変換(ＦＦＴ)のみなので離散フーリエ変換（ＤＦＴ）は未経験です。
調べて見ましたが、離散フーリエ変換の半分は対称性があるので表示する必要がないということですね。
２５６個のサンプリングなら１２８までを表示すれば良いということです。そうすれば山は１つだけになります。

>＞高速フーリエ変換では山さえ出来ません。窓関数に問題有るような気がします。
>だから質問しているのです。窓関数とは何ですか？

失礼しました窓関数そのものが無いですね。
サンプリング範囲を決めるための関数のことです。

うまく変換・表示できない原因は調べきれてないので、また待ちほど報告します。

[a5ua]

参照先のコードが間違っている可能性はないのでしょうか？

たとえば、以下のページを参考にしてFFTを実装してみたところ、O(n²)のDFTとほぼ同じ結果が得られました。
http://www.kurims.kyoto-u.ac.jp/~ooura/ ... mn1_2.html

FFTのアルゴリズムと窓関数は、直接関係するわけではないので、必要になったら意識すればよいかと思います。
（暗黙のうちに方形窓を使っていることになります。）

私の実装を載せておきます。コードはC++によるものです。
上記参照ページの、再帰呼び出しによるFFTを参考にしています。

コード:

#include "main.h"

#include <complex>
#include <vector>
#include <bitset>

static bool n_check(int n)
{
	return std::bitset<32>(n).count() <= 1;
}

static void fft_recursive(std::vector<std::complex<double>> &a)
{
	int n = a.size();
	double phase = 2 * M_PI / n;

	if (n <= 1) {
		return;
	}

	n /= 2;

	std::vector<std::complex<double>> t_l(n), t_r(n);
	
	for (int j = 0; j < n; ++j) {
		int k = n + j;
		t_l[j] = (a[j] + a[k]);
		t_r[j] = (a[j] - a[k]) * std::polar(1.0, phase * j);
	}

	fft_recursive(t_l);
	fft_recursive(t_r);

	for (int j = 0; j < n; ++j) {
		a[2 * j] = t_l[j];
		a[2 * j + 1] = t_r[j];
	}
}

static void fft(std::vector<std::complex<double>> &a)
{
	fft_recursive(a);
}


void calcSpectrum_fft(const int *signal, int length, long long *spectrum)
{
	if (!n_check(length)) {
		return;
	}
	std::vector<std::complex<double>> a(length);
	for (int i = 0; i < length; ++i) {
		a[i].real(signal[i]);
	}
	fft(a);
	for (int i = 0; i < length; ++i) {
		spectrum[i] = (long long)std::norm(a[i]);
	}
}

[みけCAT]

>h2so5さんは結果に対称性があるからいいみたいな事を言っていた気がするのですが…どうなのでしょうか？
私も遊びでフーリエ変換してみたぐらいのことしか無いですし、高速フーリエ変換(ＦＦＴ)のみなので離散フーリエ変換（ＤＦＴ）は未経験です。
調べて見ましたが、離散フーリエ変換の半分は対称性があるので表示する必要がないということですね。
２５６個のサンプリングなら１２８までを表示すれば良いということです。そうすれば山は１つだけになります。

プロジェクトにこの機能を実装しました。

参照先のコードが間違っている可能性はないのでしょうか？

たとえば、以下のページを参考にしてFFTを実装してみたところ、O(n2)のDFTとほぼ同じ結果が得られました。

http://www.kurims.kyoto-u.ac.jp/~ooura/ ... mn1_2.html
の「リスト1.2.1-3. 基数 2 の周波数間引き FFT (修正版)」
をコピペしたら離散フーリエ変換とほぼ同じ、肉眼では違いがわからない程度のスペクトラムが生成されました。
画面上部に表示している最大値が若干違うようです。

コード:

#include "main.h"

/*
http://www.kurims.kyoto-u.ac.jp/~ooura/fftman/ftmn1_2.html
リスト1.2.1-3. 基数 2 の周波数間引き FFT (修正版)
*/

void fft(int n, double theta, double ar[], double ai[])
{
	int m, mh, i, j, k, irev;
	double wr, wi, xr, xi;

	/* ---- scrambler ---- */
	i = 0;
	for (j = 1; j < n - 1; j++) {
		for (k = n >> 1; k > (i ^= k); k >>= 1);
		if (j < i) {
			xr = ar[j];
			xi = ai[j];
			ar[j] = ar[i];
			ai[j] = ai[i];
			ar[i] = xr;
			ai[i] = xi;
		}
	}
	for (mh = 1; (m = mh << 1) <= n; mh = m) {
		irev = 0;
		for (i = 0; i < n; i += m) {
			wr = cos(theta * irev);
			wi = sin(theta * irev);
			for (k = n >> 2; k > (irev ^= k); k >>= 1);
			for (j = i; j < mh + i; j++) {
				k = j + mh;
				xr = ar[j] - ar[k];
				xi = ai[j] - ai[k];
				ar[j] += ar[k];
				ai[j] += ai[k];
				ar[k] = wr * xr - wi * xi;
				ai[k] = wr * xi + wi * xr;
			}
		}
	}
}

/*nが2のn乗であるかを判定*/
static int n_check(int n) {
	if(n==0)return 1;
	while(n%2==0)n/=2;
	return n==1?1:0;
}

/*この関数を使用する*/
void calcSpectrum_fft(const int *signal, int length, long long *spectrum) {
	double* x;
	double* y;
	int i;
	/*不正*/
	if(!n_check(length)) {
		for(i=0;i<length;i++)spectrum[i]=0;
		return;
	}
	x=(double*)HeapAlloc(GetProcessHeap(),0,length*sizeof(double));
	y=(double*)HeapAlloc(GetProcessHeap(),0,length*sizeof(double));
	if(x==NULL || y==NULL) {
		if(x)HeapFree(GetProcessHeap(),0,x);
		if(y)HeapFree(GetProcessHeap(),0,y);
		return;
	}
	for(i=0;i<length;i++) {
		x[i]=signal[i];
		y[i]=0;
	}
	fft(length,2*PHI/length,x,y);
	for(i=0;i<length;i++) {
		spectrum[i]=(long long)(x[i] * x[i] + y[i] * y[i]);
	}
	HeapFree(GetProcessHeap(),0,x);
	HeapFree(GetProcessHeap(),0,y);
}

回答をくださりありがとうございました。

プロジェクト一式は日記の方にうpしました。
http://dixq.net/forum/blog.php?u=536&b=2388