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