「FFT (高速Fourier変換)」の編集履歴（バックアップ）一覧はこちら

FFT (高速Fourier変換) - (2010/03/10 (水) 13:47:41) の１つ前との変更点

追加された行は緑色になります。

削除された行は赤色になります。

#asciiart(blockquote){
/***********************************************************
	fft.c -- FFT (高速Fourier変換)
***********************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.14159265358979323846
/*
  関数 {\tt fft() }の下請けとして三角関数表を作る.
*/
static void make_sintbl(int n, float sintbl[])
 {
	int i, n2, n4, n8;
	double c, s, dc, ds, t;

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

	if (n8 != 0) sintbl[n8] = sqrt(0.5);
	for (i = 0; i < n4; i++)
		sintbl[n2 - i] = sintbl[i];
	for (i = 0; i < n2 + n4; i++)
		sintbl[i + n2] = - sintbl[i];
 }
/*
  関数 {\tt fft() }の下請けとしてビット反転表を作る.
*/
static void make_bitrev(int n, int bitrev[])
 {
	int i, j, k, n2;

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

	/* 準備 */
	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) free(sintbl);
		if (bitrev != NULL) free(bitrev);
		if (n == 0) return 0;  /* 記憶領域を解放した */
		sintbl = malloc((n + n4) * sizeof(float));
		bitrev = malloc(n * sizeof(int));
		if (sintbl == NULL || bitrev == NULL)  {
			fprintf(stderr, "記憶領域不足\n");  return 1;
		 }
		make_sintbl(n, sintbl);
		make_bitrev(n, bitrev);
	 }
	for (i = 0; i < n; 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;
		 }
	 }
	for (k = 1; k < n; k = k2)  {    /* 変換 */
		h = 0;  k2 = k + k;  d = n / k2;
		for (j = 0; j < k; j++)  {
			c = sintbl[h + n4];
			if (inverse) s = - sintbl[h];
			else         s =   sintbl[h];
			for (i = j; i < n; i += k2)  {
				ik = i + k;
				dx = s * y[ik] + c * x[ik];
				dy = c * y[ik] - s * x[ik];
				x[ik] = x[i] - dx;  x[i] += dx;
				y[ik] = y[i] - dy;  y[i] += dy;
			 }
			h += d;
		 }
	 }
	if (! inverse)    /* 逆変換でないならnで割る */
		for (i = 0; i < n; i++)  {  x[i] /= n;  y[i] /= n;   }
	return 0;  /* 正常終了 */
 }



//#define SIZE 256        // 信号サイズ
//#define PI 3.14159265    // 円周率

// ビットを左右反転した配列を返す
int* BitScrollArray(int arraySize)
 {
    int i, j;
    int* reBitArray = (int*)calloc(arraySize, sizeof(int));
    int arraySizeHarf = arraySize >> 1;

    reBitArray[0] = 0;
    for (i = 1; i < arraySize; i <<= 1)
     {
        for (j = 0; j < i; j++)
            reBitArray[j + i] = reBitArray[j] + arraySizeHarf;
        arraySizeHarf >>= 1;
     }
    return reBitArray;
 }

// FFT
void fft2(double* inputRe,
         double* inputIm,
         double* outputRe,
         double* outputIm,
         int bitSize)
 {
    int i, j, stage, type;
    int dataSize = 1 << bitSize;
    int butterflyDistance;
    int numType;
    int butterflySize;
    int jp;
    int* reverseBitArray = BitScrollArray(dataSize);
    double wRe, wIm, uRe, uIm, tempRe, tempIm, tempWRe, tempWIm;

    // バタフライ演算のための置き換え
    for (i = 0; i < dataSize; i++)
     {
        outputRe[i] = inputRe[reverseBitArray[i]];
        outputIm[i] = inputIm[reverseBitArray[i]];
     }

    // バタフライ演算
    for (stage = 1; stage <= bitSize; stage++)
     {
        butterflyDistance = 1 << stage;
        numType = butterflyDistance >> 1;
        butterflySize = butterflyDistance >> 1;

        wRe = 1.0;
        wIm = 0.0;
        uRe = cos(PI / butterflySize);
        uIm = -sin(PI / butterflySize);

        for (type = 0; type < numType; type++)
         {
            for (j = type; j < dataSize; j += butterflyDistance)
             {
                jp = j + butterflySize;
                tempRe = outputRe[jp] * wRe - outputIm[jp] * wIm;
                tempIm = outputRe[jp] * wIm + outputIm[jp] * wRe;
                outputRe[jp] = outputRe[j] - tempRe;
                outputIm[jp] = outputIm[j] - tempIm;
                outputRe[j] += tempRe;
                outputIm[j] += tempIm;
             }
            tempWRe = wRe * uRe - wIm * uIm;
            tempWIm = wRe * uIm + wIm * uRe;
            wRe = tempWRe;
            wIm = tempWIm;
         }
     }
 }

#define N 64

int main(void)
 {
	int i;
	static float x1[N], y1[N], x2[N], y2[N], x3[N], y3[N];

	for (i = 0; i < N; i++)  {
		x1[i] = x2[i] = 6 * cos( 6 * PI * i / N)
					  + 4 * sin(18 * PI * i / N);
		y1[i] = y2[i] = 0;
	 }
	//if (fft(N, x2, y2)) return EXIT_FAILURE;
	fft2();
	
	for (i = 0; i < N; i++)  {
		x3[i] = x2[i];  y3[i] = y2[i];
	 }
	if (fft(-N, x3, y3)) return EXIT_FAILURE;
	printf("      元のデータ    フーリエ変換  逆変換\n");
	for (i = 0; i < N; i++)
		printf("%4d | %6.3f %6.3f | %6.3f %6.3f | %6.3f %6.3f\n",
			i, x1[i], y1[i], x2[i], y2[i], x3[i], y3[i]);
	return EXIT_SUCCESS;
 }
}

#asciiart(blockquote){
Const PI =3.14159265358979323846
/*
  関数 {\tt fft() }の下請けとして三角関数表を作る.
*/

Sub make_sintbl(n As Long, sintbl As *Single)
	Dim i As Long, n2 As Long, n4 As Long, n8 As Long
	Dim c As Double, s As Double, dc As Double, ds As Double, t As Double

	n2 = n / 2 : n4 = n / 4 : n8 = n / 8
	t = Sin(PI / n)
	dc = 2 * t * t : ds = Sqr(dc * (2 - dc))
	t = 2 * dc : c = sintbl[n4] = 1 : sintbl[0] = 0 : s = sintbl[0]
	For i = 0 To n4-1
	   sintbl[n2 - i] =   sintbl[i]
	   sintbl[i + n2] = - sintbl[i]
	Next

	If n8 <> 0 Then sintbl[n8] = Sqr(0.5)
	For i = 0 To i< n4-1
		sintbl[n2 - i] = sintbl[i]
	Next
	For i = 0 To n2 + n4-1
		sintbl[i + n2] = - sintbl[i]
	Next
End Sub
		
/*
  関数 {\tt fft() }の下請けとしてビット反転表を作る.
*/
Sub make_bitrev(n As Long, bitrev  As *Long)
	Dim i As Long, j As Long, k As Long,n2 As Long

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

Function fft(n As Long,x As *Single, y As *Single) As Long
	Dim i As Long, j As Long, k As Long, ik As Long, h As Long, d As Long
	Dim k2 As Long, n4 As Long, inverse As Long
	DIm t As Single, s As Single, c As Single, dx As Single, dy As Single
	
	/* 準備 */
	If n < 0 Then
		n = -n : inverse = 1  /* 逆変換 */
	Else
		inverse = 0
	End If
	n4 = n / 4
	If n <> last_n Or n = 0 Then
		last_n = n
		If sintbl <> NULL Then free(sintbl)
		If bitrev <> NULL Then free(bitrev)
		If n = 0 Then Exit Function  /* 記憶領域を解放した */
		sintbl = malloc((n + n4) * SizeOf(Single))
		bitrev = malloc(n * SizeOf(Long))
		If sintbl = NULL Or bitrev = NULL Then
			Print "記憶領域不足"
			fft = 1
			Exit Function
		End If
		make_sintbl(n, sintbl)
		make_bitrev(n, bitrev)
	End If
	For i = 0 To n-1    /* ビット反転 */
		j = bitrev[i]
		If i < j Then
			t = x[i] : x[i] = x[j] : x[j] = t
			t = y[i] : y[i] = y[j] : y[j] = t
		 End If
	Next
	k = 1
	While k < n /* 変換 */
		h = 0 : k2 = k + k : d = n / k2
		For j = 0 To k-1
			c = sintbl[h + n4]
			If inverse <> 0 Then
				s = - sintbl[h]
			Else
				s =   sintbl[h]
			End If
			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
			Wend
			h = h + d
		Next
		k = k2    
	Wend
	If inverse = 0 Then    /* 逆変換でないならnで割る */
		For i = 0 To i = n-1
			x[i] = x[i] / n : y[i] = y[i] / n
	 	Next
	End If
	fft = 0/* 正常終了 */
End Function


Const N = 64
Declare Function sprintf CDECL Lib"msvcrt" (p As *Byte, fmt As *Byte, ...) As Long
Sub main()
	Dim i As Long
	Dim x1[ELM(N)] As Single, x2[ELM(N)] As Single, x3[ELM(N)] As Single
	Dim y1[ELM(N)] As Single, y2[ELM(N)] As Single, y3[ELM(N)] As Single

	
	For i = 0 To N-1
		x1[i] = 6 * Cos( 6 * PI * i / N) + 4 * Sin(18 * PI * i / N)
		x2[i] = x1[i]
		y1[i] = 0 : y2[i] = y1[i] 
	Next
	If fft(N, x2, y2) Then Exit Sub
	
	For i = 0 To N-1
		x3[i] = x2[i] : y3[i] = y2[i]
	Next
	If fft(-N, x3, y3) Then Exit Sub
	Dim mes[589] As Byte
	sprintf(mes, Ex"      元のデータ    フーリエ変換  逆変換")
	Print MakeStr(mes)
	For i = 0 To N-1
		sprintf(mes, Ex"%4d | %6.3f %6.3f | %6.3f %6.3f | %6.3f %6.3f", _
			i, x1[i] As Double, y1[i] As Double, x2[i] As Double, y2[i] As Double, x3[i] As Double, y3[i] As Double)
		Print MakeStr(mes)
	Next
End Sub

#N88BASIC
main()

}