|
ありがとうございます。
void fft(int n, double theta, double ar[], double ai[],
double tmpr[], double tmpi[])
{
int radix, n_radix, j, m, r;
double xr, xi, wr, wi;
if (n <= 1) return;
/* ---- factorization ---- */
for (radix = 2; radix * radix <= n; radix++) {
if (n % radix == 0) break;
}
if (n % radix != 0) radix = n;
n_radix = n / radix;
/* ---- butterflies ---- */
for (j = 0; j < n_radix; j++) {
for (m = 0; m < radix; m++) {
xr = ar[j];
xi = ai[j];
for (r = n_radix; r < n; r += n_radix) {
wr = cos(theta * m * r);
wi = sin(theta * m * r);
xr += wr * ar[r + j] - wi * ai[r + j];
xi += wr * ai[r + j] + wi * ar[r + j];
}
wr = cos(theta * m * j);
wi = sin(theta * m * j);
tmpr[m * n_radix + j] = xr * wr - xi * wi;
tmpi[m * n_radix + j] = xi * wr + xr * wi;
}
}
for (r = 0; r < n; r += n_radix) {
fft(n_radix, theta * radix, &tmpr[r], &tmpi[r], ar, ai);
}
for (j = 0; j < n_radix; j++) {
for (m = 0; m < radix; m++) {
ar[radix * j + m] = tmpr[n_radix * m + j];
ai[radix * j + m] = tmpi[n_radix * m + j];
}
}
}
これを以下のように翻訳してみました。
Public Sub one_D_FFT(ByRef xr#(), ByRef xi#())
Dim n As Integer = xr.GetLength(0) - 1
Dim theta = 2 * Math.PI / n
Dim tmpr(n) As Double, tmpi(n) As Double
Call fft(n, theta, xr, xi, tmpr, tmpi, 0)
End Sub
Private Sub fft(ByVal n%, ByVal theta#, ByRef ar#(), ByRef ai#(), ByRef tmpr#(), ByRef tmpi#(), ByVal r_st%)
If n <= 1 Then Exit Sub
Dim radix As Integer = 2
Do While radix * radix <= n
If n Mod radix = 0 Then Exit Do
radix += 1
Loop
If n Mod radix <> 0 Then radix = n
Dim n_radix As Integer = n \ radix
For j As Integer = 0 To n_radix - 1
For m As Integer = 0 To radix - 1
Dim xr As Double = ar(r_st + j)
Dim xi As Double = ai(r_st + j)
For r As Integer = n_radix To n - 1 Step n_radix
Dim wr0 As Double = Math.Cos(theta * m * r)
Dim wi0 As Double = Math.Sin(theta * m * r)
xr += (wr0 * ar(r_st + r + j) - wi0 * ai(r_st + r + j))
xi += (wr0 * ai(r_st + r + j) + wi0 * ar(r_st + r + j))
Next r
Dim wr As Double = Math.Cos(theta * m * j)
Dim wi As Double = Math.Sin(theta * m * j)
tmpr(m * n_radix + j) = xr * wr - xi * wi
tmpi(m * n_radix + j) = xi * wr + xr * wi
Next m
Next j
For r As Integer = 0 To n - 1 Step n_radix
fft(n_radix, theta * radix, tmpr, tmpi, ar, ai, r)
Next r
For j As Integer = 0 To n_radix - 1
For m As Integer = 0 To radix - 1
ar(r_st + radix * j + m) = tmpr(n_radix * m + j)
ai(r_st + radix * j + m) = tmpi(n_radix * m + j)
Next m
Next j
End Sub
しかし、むちゃくちゃな結果しか得られないのですが・・・・。
どこが間違っていますでしょうか?
|