Skip to main content ↑ ↓

Notice

Please note that most of the software linked on this forum is likely to be safe to use. If you are unsure, feel free to ask in the relevant topics, or send a private message to an administrator or moderator. To help curb the problems of false positives, or in the event that you do find actual malware, you can contribute through the article linked here.
Topic: Strange FFT results (Read 3815 times) previous topic - next topic
0 Members and 1 Guest are viewing this topic.

Strange FFT results

I'm using an FFT algorithm I have found on dspguide.com, but converted in vb2008. Here it is:
Code: [Select]
    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: [Select]
        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?

Strange FFT results

Reply #1 –
Quote from: Totem on
I'm using an FFT algorithm I have found on dspguide.com, but converted in vb2008.

Quote from: Totem on
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: [Select]
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
 

Strange FFT results

Reply #2 –
Quote from: SebastianG on
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