Fast Fourier transform

From ScienceZero
Jump to: navigation, search
'Complex radix 2 forward Fast Fourier Transform
'Array size is restricted to 2m and the number of operations is O*n*log2(n).

  n = 2 ^ m

  j = 0
  For i = 0 To n - 2
     If (i < j) Then
        tmpR = dataR(i)
        tmpI = dataI(i)
        dataR(i) = dataR(j)
        dataI(i) = dataI(j)
        dataR(j) = tmpR
        dataI(j) = tmpI
     End If
     k = n / 2
     While (k <= j)
        j = j - k
        k = k / 2
     Wend
     j = j + k
  Next i
   
  c1 = -1
  c2 = 0
  l2 = 1
  For l = 0 To m - 1
     l1 = l2
     l2 = l2 + l2
     u1 = 1
     u2 = 0
     For j = 0 To l1 - 1
       For i = j To n - 1 Step l2
           i1 = i + l1
           t1 = u1 * dataR(i1) - u2 * dataI(i1)
           t2 = u1 * dataI(i1) + u2 * dataR(i1)
           dataR(i1) = dataR(i) - t1
           dataI(i1) = dataI(i) - t2
           dataR(i) = dataR(i) + t1
           dataI(i) = dataI(i) + t2
       Next i
       z = u1 * c1 - u2 * c2
       u2 = u1 * c2 + u2 * c1
       u1 = z
    Next j
    c2 = -Sqr((1 - c1) / 2)
    c1 = Sqr((1 + c1) / 2)
   Next l

  For i = 0 To n - 1
     dataR(i) = dataR(i) / n
     dataI(i) = dataI(i) / n
  Next i


'Complex radix 2 reverse Fast Fourier Transform
'Array size is restricted to 2m and the number of operations is O*n*log2(n).

  n = 2 ^ m

  j = 0
  For i = 0 To n - 2
     If (i < j) Then
        tmpR = dataR(i)
        tmpI = dataI(i)
        dataR(i) = dataR(j)
        dataI(i) = dataI(j)
        dataR(j) = tmpR
        dataI(j) = tmpI
     End If
     k = n / 2
     While (k <= j)
        j = j - k
        k = k / 2
     Wend
     j = j + k
  Next i
   
  c1 = -1
  c2 = 0
  l2 = 1
  For l = 0 To m - 1
     l1 = l2
     l2 = l2 + l2
     u1 = 1
     u2 = 0
     For j = 0 To l1 - 1
       For i = j To n - 1 Step l2
           i1 = i + l1
           t1 = u1 * dataR(i1) - u2 * dataI(i1)
           t2 = u1 * dataI(i1) + u2 * dataR(i1)
           dataR(i1) = dataR(i) - t1
           dataI(i1) = dataI(i) - t2
           dataR(i) = dataR(i) + t1
           dataI(i) = dataI(i) + t2
       Next i
       z = u1 * c1 - u2 * c2
       u2 = u1 * c2 + u2 * c1
       u1 = z
    Next j
    c2 = Sqr((1 - c1) / 2)
    c1 = Sqr((1 + c1) / 2)
   Next l