I'm using an FFT algorithm I have found on dspguide.com, but converted in vb2008. Here it is:
Public Shared Sub FFT2(ByRef ReX() As Single, ByRef ImX() As Single, ByVal Length As Int32)
Dim nm1 As Int32 = Length - 1
Dim nd2 As Int32 = Length / 2
Dim m As Int32 = Math.Log(Length) / Math.Log(2)
Dim j As Int32 = nd2
Dim tr, ti As Single
Dim k As Int32
For I As Int32 = 1 To Length - 2
If I < j Then
tr = ReX(j)
ti = ImX(j)
ReX(j) = ReX(I)
ImX(j) = ImX(I)
ReX(I) = tr
ImX(I) = ti
End If
k = nd2
Do While k <= j
j -= k
k /= 2
Loop
j += k
Next
Dim le, jm1 As Single
Dim le2, ip As Single
Dim ur, ui As Single
Dim sr, si As Single
For l As Int32 = 1 To m
le = 2 ^ l
le2 = le / 2
ur = 1
ui = 0
sr = Math.Cos(Math.PI / le2)
si = -Math.Sin(Math.PI / le2)
For j = 1 To le2
jm1 = j - 1
For i As Int32 = jm1 To nm1 Step le
ip = i + le2
tr = ReX(ip) * ur - ImX(ip) * ui
ti = ReX(ip) * ui - ImX(ip) * ur
ReX(ip) = ReX(i) - tr
ImX(ip) = ImX(i) - ti
ReX(i) += tr
ImX(i) += ti
Next
tr = ur
ur = tr * sr - ui * si
ui = tr * si - ui * sr
Next
Next
End Sub
It works fine if the input signal (which is passed as ReX parameter) is a single impulse: the fft produces a ReX that contains the same value in all its cells. Recombining ReX and ImX produces the input signal as it was (even if there is a small, unevitable, noise). But, if the input is not an impulse, never matter how simple it is, the fft outpus completely wrong results, and the recombined waveform is absolutely different from the initial one. The code I used is the following:
Dim ReX(), ImX() As Single
ReDim ReX(2 ^ 8 - 1)
ReDim ImX(2 ^ 8 - 1)
Dim I As Int32 = 0
For I = 0 To 255
ReX(I) = 0
ImX(I) = 0
Next
ReX(0) = 50 * Math.Cos(0)
ReX(1) = 50 * Math.Cos(Math.PI / 16)
ReX(2) = 50 * Math.Cos(Math.PI / 8)
Stop
FFT.FFT2(ReX, ImX, ReX.Length)
Dim x(ReX.Length - 1) As Single
For I = 0 To 4
x(I) = (ReX(0) / (x.Length)) * Math.Cos(0 * 2 * Math.PI * I / (x.Length - 1)) + _
(ImX(0) / (x.Length)) * Math.Sin(0 * 2 * Math.PI * I / (x.Length - 1))
For J As Int32 = 1 To x.Length / 2 - 1
x(I) += (ReX(J) / (x.Length / 2)) * Math.Cos(J * 2 * Math.PI * I / (x.Length - 1)) + _
(ImX(J) / (x.Length / 2)) * Math.Sin(J * 2 * Math.PI * I / (x.Length - 1))
Next
x(I) += (ReX(x.Length / 2) / (x.Length)) * Math.Cos(x.Length * Math.PI * I / (x.Length - 1)) + _
(ImX(x.Length / 2) / (x.Length)) * Math.Sin(x.Length * Math.PI * I / (x.Length - 1))
Next
Stop
The Stop statements are used for debugging. x should contains the same values containted in ReX before passing it to fft method (i.e. the first three points of a cosine wave of 1Hz and 50x amplitude). Could you tell me what's wrong?