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

「FFT (高速Fourier変換)」(2010/03/10 (水) 13:52:43) の最新版変更点

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

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

#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()

}

#asciiart(blockquote){
Const PI =3.14159265358979323846
/*
  関数 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
		
/*
  関数 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のアルゴリズム).
  標本点の数  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[] の値は変わらない).
*/
Dim last_n As Long /* 前回呼出し時の  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()

}