Hydrogenaudio Forums > Strange FFT results

Totem

Mar 6 2009, 21:00

I'm using an FFT algorithm I have found on dspguide.com, but converted in vb2008. Here it is:

CODE

    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:

CODE

        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?

SebastianG

Mar 16 2009, 10:15

QUOTE (Totem @ Mar 6 2009, 21:00)

I'm using an FFT algorithm I have found on dspguide.com, but converted in vb2008.

QUOTE (Totem @ Mar 6 2009, 21:00)

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.

It looks like you are not aware of the fact that this FFT algorithm returns a permutation of what you're interested in. If you want to reorder your elements you need to swap some elements: For every pair of indices (i,j) with i<j and binary(i)==reverse(binary(j)) swap elements at index i and j. Here, "binary" returns a binary string of length log_2(n) and "reverse" reverses the string.

Example: Elements to swap for n = 16 are
(1,8), (2,4), (3,12), (5,10), (7,14), (11,13)

you would code this permutation like this:

CODE

for (int i=0, j=0; i<n; ++i) {
  if (i<j) {
    swap(data[i],data[j]);
  }
  for (int k=n/2; k>0;) {
    j ^= k;
    if ((j & k)==0) {
      k = k >> 1;
    } else {
      break;
    }
  }
}

In Basic you would probably write "j = j xor k" instead of "j ^= k;"

Note: Indices start from 0 (ZERO) ;-)
Note2: n needs to be a power of two

Cheers!
SG

SebastianG

Mar 16 2009, 13:46

QUOTE (SebastianG @ Mar 16 2009, 10:15)

you would code this permutation like this:

Sorry, for the trouble. I just saw that this permutation is part of your FFT code.

Still, the indexing looks weird. Could it be that you're off by one?

Cheers!
SG