# VBForums CodeBank > CodeBank - Visual Basic 6 and earlier >  Simple FFT module

## Ben321

I wrote this FFT that you can put in a module file and use anywhere in your program.


```
Public Const Pi As Double = 3.14159265358979

Public Type Complex
    Re As Double
    Im As Double
End Type


Dim FFTSize As Long
Dim FFTSize_Log2 As Long
Dim BOR_LUT() As Long
Dim Pow2_LUT() As Long


Public Sub SetupFFT(ByVal Size As Long)
    Dim n As Long
    Dim n1 As Long
    Dim n2 As Long
    Dim i As Long
    
    FFTSize = Size
    FFTSize_Log2 = Log(FFTSize) / Log(2)
    
    ReDim BOR_LUT(FFTSize - 1)
    
    For n = 0 To FFTSize - 1
        n1 = n
        n2 = 0
        For i = 0 To FFTSize_Log2 - 1
            n2 = (n2 * 2) Or (n1 And 1)
            n1 = n1 \ 2
        Next i
        BOR_LUT(n) = n2
    Next n
    
    ReDim Pow2_LUT(FFTSize_Log2 - 1)
    For n = 0 To FFTSize_Log2 - 1
        Pow2_LUT(n) = 2 ^ n
    Next n
End Sub


Public Sub FFT(ByRef Data1() As Complex, ByRef Data2() As Complex, ByVal Sign As Long)
    Dim Layer As Long
    Dim BlockCount As Long
    Dim BlockSize As Long
    Dim HalfBlockSize As Long
    Dim BlockNumber As Long
    Dim ButterflyNumber As Long
    Dim n0 As Long
    Dim n1 As Long
    Dim n2 As Long
    Dim DataTemp() As Complex
    Dim Twiddles() As Complex
    Dim n As Long
    
    ReDim Twiddles(FFTSize - 1)
    For n = 0 To FFTSize - 1
        Twiddles(n) = GetTwiddle(n, Sign)
    Next n
    
    
    DataTemp() = Data1()
    
    For Layer = 0 To FFTSize_Log2 - 1
        BlockCount = Pow2_LUT(Layer)
        BlockSize = FFTSize \ BlockCount
        HalfBlockSize = BlockSize \ 2
        For BlockNumber = 0 To BlockCount - 1
            n0 = BlockNumber * BlockSize
            For ButterflyNumber = 0 To HalfBlockSize - 1
                n1 = n0 + ButterflyNumber
                n2 = n1 + HalfBlockSize
                
                'do butterfly calculation
                Data2(n1) = ComplexAdd(DataTemp(n1), DataTemp(n2))
                Data2(n2) = ComplexMult(ComplexSub(DataTemp(n1), DataTemp(n2)), Twiddles(ButterflyNumber * BlockCount))
                
            Next ButterflyNumber
        Next BlockNumber
        DataTemp() = Data2()
    Next Layer
    
    For n = 0 To FFTSize - 1
        Data2(n) = DataTemp(BOR_LUT(n))
    Next n
End Sub


Private Function ComplexAdd(ByRef Value1 As Complex, ByRef Value2 As Complex) As Complex
    ComplexAdd.Re = Value1.Re + Value2.Re
    ComplexAdd.Im = Value1.Im + Value2.Im
End Function

Private Function ComplexSub(ByRef Value1 As Complex, ByRef Value2 As Complex) As Complex
    ComplexSub.Re = Value1.Re - Value2.Re
    ComplexSub.Im = Value1.Im - Value2.Im
End Function

Private Function ComplexMult(ByRef Value1 As Complex, ByRef Value2 As Complex) As Complex
    ComplexMult.Re = Value1.Re * Value2.Re - Value1.Im * Value2.Im
    ComplexMult.Im = Value1.Re * Value2.Im + Value1.Im * Value2.Re
End Function


Private Function GetTwiddle(ByVal k As Long, ByVal Sign As Long) As Complex
    With GetTwiddle
        .Re = Cos(2 * Pi * k / FFTSize)
        .Im = Sin(2 * Pi * k / FFTSize) * Sign
    End With
End Function
```


You just call SetupFFT before using it, to set the size of the FFT (must be a power-of-2), and then you perform an FFT or IFFT on your data by calling the FFT sub. You also need to call SetupFFT again, any time you want to change the size of your FFT. I moved as much stuff out of the FFT sub to the SetupFFT sub, so that this setup stuff doesn't get run with every single you run the FFT sub (which would slow it down). To do an FFT, call it with the Sign parameter set to -1. To do an IFFT, call the FFT sub with the Sign parameter set to 1. The FFT parameters Data1 and Data2 are arrays of type Complex (a user-defined type). Data1 is the input, and Data2 is the output. Data1 is time domain, and Data2 is frequency domain, when performing an FFT. Data1 is frequency domain, and Data2 is time domain, when performing an IFFT.

----------


## Ben321

I've made a FixedPoint version of it, using a userdefined type called FP1000, which is a fixed point using 1000 as the denominator, and the complex type is now called FPComplex. All the fixed point and fixed point complex functions and userdefined types have been moved into their own separate module now.


```
Public Type FP1000 'fixed point with denominator of 1000
    Numerator As Long
End Type

Public Type FPComplex
    Re As FP1000
    Im As FP1000
End Type


Public Function DoubleToFP(ByVal Value As Double) As FP1000
    DoubleToFP.Numerator = Value * 1000
End Function

Public Function FPtoDouble(ByRef Value As FP1000) As Double
    FPtoDouble = Value.Numerator / 1000
End Function



Public Function LongToFP(ByVal Value As Long) As FP1000
    LongToFP.Numerator = Value * 1000
End Function

Public Function FPtoLong(ByRef Value As FP1000) As Long
    FPtoLong = Value.Numerator \ 1000
End Function



Public Function FPAdd(ByRef Val1 As FP1000, ByRef Val2 As FP1000) As FP1000
    FPAdd.Numerator = Val1.Numerator + Val2.Numerator
End Function

Public Function FPSub(ByRef Val1 As FP1000, ByRef Val2 As FP1000) As FP1000
    FPSub.Numerator = Val1.Numerator - Val2.Numerator
End Function

Public Function FPMult(ByRef Val1 As FP1000, ByRef Val2 As FP1000) As FP1000
    FPMult.Numerator = (((Val1.Numerator) \ 10) * ((Val2.Numerator) \ 10)) \ 10
End Function

Public Function FPDiv(ByRef Val1 As FP1000, ByRef Val2 As FP1000) As FP1000
    FPDiv.Numerator = (Val1.Numerator * 1000) \ Val2.Numerator
End Function




Public Function FPCAdd(ByRef Val1 As FPComplex, ByRef Val2 As FPComplex) As FPComplex
    With FPCAdd
        .Re = FPAdd(Val1.Re, Val2.Re)
        .Im = FPAdd(Val1.Im, Val2.Im)
    End With
End Function

Public Function FPCSub(ByRef Val1 As FPComplex, ByRef Val2 As FPComplex) As FPComplex
    With FPCSub
        .Re = FPSub(Val1.Re, Val2.Re)
        .Im = FPSub(Val1.Im, Val2.Im)
    End With
End Function

Public Function FPCMult(ByRef Val1 As FPComplex, ByRef Val2 As FPComplex) As FPComplex
    With FPCMult
        .Re = FPSub(FPMult(Val1.Re, Val2.Re), FPMult(Val1.Im, Val2.Im))
        .Im = FPAdd(FPMult(Val1.Re, Val2.Im), FPMult(Val1.Im, Val2.Re))
    End With
End Function
```

And the FFT code has been changed to use fixed point math in the butterfly calculations specifically, as it's the most computationally intensive, and will benefit most from using the CPU's faster integer math instructions (instead of the slower floating point math instructions). Here's the code for the new FFT module that uses this fixed point math.


```
Public Const Pi As Double = 3.14159265358979


Dim FFTSize As Long
Dim FFTSize_Log2 As Long
Dim BOR_LUT() As Long
Dim Pow2_LUT() As Long


Public Sub SetupFFT(ByVal Size As Long)
    Dim n As Long
    Dim n1 As Long
    Dim n2 As Long
    Dim i As Long
    
    FFTSize = Size
    FFTSize_Log2 = Log(FFTSize) / Log(2)
    
    ReDim BOR_LUT(FFTSize - 1)
    
    For n = 0 To FFTSize - 1
        n1 = n
        n2 = 0
        For i = 0 To FFTSize_Log2 - 1
            n2 = (n2 * 2) Or (n1 And 1)
            n1 = n1 \ 2
        Next i
        BOR_LUT(n) = n2
    Next n
    
    ReDim Pow2_LUT(FFTSize_Log2 - 1)
    For n = 0 To FFTSize_Log2 - 1
        Pow2_LUT(n) = 2 ^ n
    Next n
End Sub


Public Sub FFT(ByRef Data1() As FPComplex, ByRef Data2() As FPComplex, ByVal Sign As Long)
    Dim Layer As Long
    Dim BlockCount As Long
    Dim BlockSize As Long
    Dim HalfBlockSize As Long
    Dim BlockNumber As Long
    Dim ButterflyNumber As Long
    Dim n0 As Long
    Dim n1 As Long
    Dim n2 As Long
    Dim DataTemp() As FPComplex
    Dim Twiddles() As FPComplex
    Dim n As Long
    
    ReDim Twiddles(FFTSize - 1)
    For n = 0 To FFTSize - 1
        Twiddles(n) = GetTwiddle(n, Sign)
    Next n
    
    
    DataTemp() = Data1()
    
    For Layer = 0 To FFTSize_Log2 - 1
        BlockCount = Pow2_LUT(Layer)
        BlockSize = FFTSize \ BlockCount
        HalfBlockSize = BlockSize \ 2
        For BlockNumber = 0 To BlockCount - 1
            n0 = BlockNumber * BlockSize
            For ButterflyNumber = 0 To HalfBlockSize - 1
                n1 = n0 + ButterflyNumber
                n2 = n1 + HalfBlockSize
                
                'do butterfly calculation
                Data2(n1) = FPCAdd(DataTemp(n1), DataTemp(n2))
                Data2(n2) = FPCMult(FPCSub(DataTemp(n1), DataTemp(n2)), Twiddles(ButterflyNumber * BlockCount))
                
            Next ButterflyNumber
        Next BlockNumber
        DataTemp() = Data2()
    Next Layer
    
    For n = 0 To FFTSize - 1
        Data2(n) = DataTemp(BOR_LUT(n))
    Next n
End Sub


Private Function GetTwiddle(ByVal k As Long, ByVal Sign As Long) As FPComplex
    With GetTwiddle
        .Re = DoubleToFP(Cos(2 * Pi * k / FFTSize))
        .Im = DoubleToFP(Sin(2 * Pi * k / FFTSize) * Sign)
    End With
End Function
```

----------


## wqweto

Is Data2 the output and Data1 the input? Can you name these accordingly? Are +1 and -1 only valid values for Sign? Can you make it a boolean (named Negative or Positive) so that incorrect values are impossible to pass?

The "speed up" in the second implementation approach is dubious provided that you constantly (like in the innermost loop) call GetTwiddle which calls both Sin and Cos which are orders of magnitude slower that f-point multiplication (or division).

You have to concentrate on optimizing part of the code where speed gains are possible. This most likely includes precomputing LUTs for the GetTwiddle and stop calling this function in the innermost loop.

cheers,
</wqw>

----------


## The trick

*CFFT* class:


```
Option Explicit

Private Declare Function memcpy Lib "kernel32" _
                         Alias "RtlMoveMemory" ( _
                         ByRef Destination As Any, _
                         ByRef Source As Any, _
                         ByVal Length As Long) As Long
                         
Private Const PI            As Double = 3.14159265358979
Private Const MAX_FFT_LOG   As Long = 20
Private Const MAX_FFT       As Long = 2 ^ MAX_FFT_LOG

Private m_dRotCoef()    As Double
Private m_lFFTSize      As Long
Private m_lFFTLog2      As Long
Private m_fBuffer()     As Single

Public Property Get FFTSize() As Long
    FFTSize = m_lFFTSize
End Property

Public Property Let FFTSize( _
                    ByVal lValue As Long)
    Dim lIndex  As Long
    Dim lLog2   As Long
    
    If lValue <= 1 Or lValue > MAX_FFT Then
        Err.Raise 5
    End If
    
    lLog2 = 2
    
    For lIndex = 1 To MAX_FFT_LOG
            
        If lLog2 >= lValue Then
            
            m_lFFTSize = lLog2
            m_lFFTLog2 = lIndex
            Exit For
            
        End If
        
        lLog2 = lLog2 * 2
            
    Next
    
    ReDim m_fBuffer(1, m_lFFTSize - 1)
    
End Property

Public Sub FFT( _
           ByRef fSignal() As Single)
    Dim lCount  As Long
    
    lCount = UBound(fSignal, 2) + 1
    
    If lCount > m_lFFTSize Or UBound(fSignal, 1) <> 1 Then
        Err.Raise 5
    End If

    FFT_ fSignal, False

End Sub

Public Sub IFFT( _
           ByRef fSignal() As Single)
    Dim lCount  As Long
    
    lCount = UBound(fSignal, 2) + 1
    
    If lCount > m_lFFTSize Or UBound(fSignal, 1) <> 1 Then
        Err.Raise 5
    End If

    FFT_ fSignal, True

End Sub

Public Sub FFTR( _
           ByRef fSignal() As Single, _
           ByRef fResult() As Single)
    Dim lIndex  As Long
    Dim lCount  As Long
    
    lCount = UBound(fSignal) + 1
    
    If lCount > m_lFFTSize Then
        Err.Raise 5
    End If
    
    For lIndex = 0 To lCount - 1
        m_fBuffer(0, lIndex) = fSignal(lIndex)
        m_fBuffer(1, lIndex) = 0
    Next
    
    FFT_ m_fBuffer, False
    
    fResult = m_fBuffer
    
End Sub

Public Sub IFFTR( _
           ByRef fSignal() As Single, _
           ByRef fResult() As Single)
    Dim lIndex  As Long
    Dim lCount  As Long
    
    lCount = UBound(fSignal, 2) + 1
    
    If lCount > m_lFFTSize Or UBound(fSignal, 1) <> 1 Then
        Err.Raise 5
    End If
    
    memcpy m_fBuffer(0, 0), fSignal(0, 0), lCount * 8
    
    FFT_ m_fBuffer, True
    
    ReDim fResult(lCount - 1)
    
    For lIndex = 0 To lCount - 1
        fResult(lIndex) = m_fBuffer(0, lIndex)
    Next

End Sub

Private Sub FFT_( _
            ByRef fData() As Single, _
            ByVal bIsInverse As Boolean)
    Dim i   As Long:    Dim j   As Long
    Dim n   As Long:    Dim K   As Long
    Dim io  As Long:    Dim ie  As Long
    Dim in_ As Long:    Dim nn  As Long
    Dim ur  As Single:  Dim ui  As Single
    Dim tpr As Single:  Dim tpi As Single
    Dim tqr As Single:  Dim tqi As Single
    Dim wr  As Single:  Dim wi  As Single
    Dim sr  As Single:  Dim ti  As Long
    Dim tr  As Long

    nn = m_lFFTSize \ 2: ie = m_lFFTSize
    
    For n = 1 To m_lFFTLog2
    
        wr = m_dRotCoef(0, m_lFFTLog2 - n): wi = m_dRotCoef(1, m_lFFTLog2 - n)
        
        If bIsInverse Then
            wi = -wi
        End If
        
        in_ = ie \ 2: ur = 1: ui = 0
        
        For j = 0 To in_ - 1
        
            For i = j To m_lFFTSize - 1 Step ie
            
                io = i + in_
                tpr = fData(0, i) + fData(0, io): tpi = fData(1, i) + fData(1, io)
                tqr = fData(0, i) - fData(0, io): tqi = fData(1, i) - fData(1, io)
                fData(0, io) = tqr * ur - tqi * ui: fData(1, io) = tqi * ur + tqr * ui
                fData(0, i) = tpr: fData(1, i) = tpi
                
            Next
            
            sr = ur: ur = ur * wr - ui * wi: ui = ui * wr + sr * wi
            
        Next
        
        ie = ie \ 2
        
    Next

    j = 1
    
    For i = 1 To m_lFFTSize - 1
    
        If i < j Then
        
            io = i - 1: in_ = j - 1: tpr = fData(0, in_): tpi = fData(1, in_)
            fData(0, in_) = fData(0, io): fData(1, in_) = fData(1, io)
            fData(0, io) = tpr: fData(1, io) = tpi
            
        End If
        
        K = nn
        
        Do While K < j
            j = j - K: K = K \ 2
        Loop
        
        j = j + K
        
    Next
    
    If bIsInverse Then
        Exit Sub
    End If
    
    wr = 1 / m_lFFTSize
    
    For i = 0 To m_lFFTSize - 1
        fData(0, i) = fData(0, i) * wr: fData(1, i) = fData(1, i) * wr
    Next

End Sub

Private Sub Class_Initialize()
    Dim lIndex  As Long
    Dim dStep   As Double
    
    ReDim m_dRotCoef(1, MAX_FFT_LOG - 1)
    
    dStep = 1
    
    For lIndex = 0 To MAX_FFT_LOG - 1
        
        m_dRotCoef(0, lIndex) = Cos(dStep * PI)
        m_dRotCoef(1, lIndex) = -Sin(dStep * PI)
        
        dStep = dStep / 2

    Next
    
    FFTSize = 1024
    
End Sub
```

----------


## Ben321

> Is Data2 the output and Data1 the input? Can you name these accordingly? Are +1 and -1 only valid values for Sign? Can you make it a boolean (named Negative or Positive) so that incorrect values are impossible to pass?


Data1 and Data2 were names I picked back when I first designed the function, as just temporary names. Obviously DataIn and DataOut would be better, but Data1 and Data2 were shorter and quicker to type so I used it during rapid development of the function. I just haven't taken the time to change them after that.

Sign values of 1 and -1 are the only valid values. I intentionally avoided boolean, because with boolean I would need to use an if-then statement, which means at the machine-code level using cmp and jmp to analyze the Sign argument. Such comparison and jump instructions I think may be slower than directly using the value of sign in the equation (though as it's only used in initially calculating the twiddle factors, not during the main loop, it may not be that significant of a speedup).




> The "speed up" in the second implementation approach is dubious provided that you constantly (like in the innermost loop) call GetTwiddle which calls both Sin and Cos which are orders of magnitude slower that f-point multiplication (or division).


Look more carefully at the code. GetTwiddle is only used in the loop to initialize the Twiddles array. All use of twiddle factors in the main processing loop is done by accessing the values already stored in the Twiddles array. And that was even before I converted the math done in the main loop to being entirely fixed-point math. So yes, this speed-up by switching to fixed-point math is actually huge. And it becomes an even bigger speed-up after compiling to an EXE, especially if you disable integer overflow checks and array bounds checks. I would not be surprised at all if this sped up version, at least in a compiled EXE, could handle the real-time processing of high sample rate data (like the 2.4 MSPS data stream taken from an RTL-SDR, software defined radio device).




> You have to concentrate on optimizing part of the code where speed gains are possible. This most likely includes precomputing LUTs for the GetTwiddle and stop calling this function in the innermost loop.


That was already done in my code. Please check it again more carefully. Twiddles is not the same as GetTwiddle. GetTwiddle is a function. Twiddles is the array that is initialized in a separate loop using the GetTwiddle function. GetTwiddle is never called in the main processing loop.

----------


## Ben321

The only thing I see I could optimize more is to make it so that Twiddles only ever needs to get calculated once period, instead of once with each call to the FFT function. I had it that way because Twiddles will be different between performing an FFT and an IFFT, and so if you are doing an FFT first and then an IFFT, it will need to calculate the FFT twiddles during the first call, but calculate the IFFT twiddles to the next call of the function. I'll see if I can figure out a way to avoid doing that, and have it store the twiddles in a way that doesn't need to be re-calculated for every call to the FFT function.

----------


## Ben321

Here I've made that improvement to the fixed point FFT.


```
Public Const Pi As Double = 3.14159265358979


Dim FFTSize As Long
Dim FFTSize_Log2 As Long
Dim BOR_LUT() As Long
Dim Pow2_LUT() As Long


Public Sub SetupFFT(ByVal Size As Long, ByRef TwiddlesFFT() As FPComplex, ByRef TwiddlesIFFT() As FPComplex)
    Dim n As Long
    Dim n1 As Long
    Dim n2 As Long
    Dim i As Long
    
    FFTSize = Size
    FFTSize_Log2 = Log(FFTSize) / Log(2)
    
    ReDim BOR_LUT(FFTSize - 1)
    
    For n = 0 To FFTSize - 1
        n1 = n
        n2 = 0
        For i = 0 To FFTSize_Log2 - 1
            n2 = (n2 * 2) Or (n1 And 1)
            n1 = n1 \ 2
        Next i
        BOR_LUT(n) = n2
    Next n
    
    ReDim Pow2_LUT(FFTSize_Log2 - 1)
    For n = 0 To FFTSize_Log2 - 1
        Pow2_LUT(n) = 2 ^ n
    Next n
    
    ReDim TwiddlesFFT(FFTSize - 1)
    For n = 0 To FFTSize - 1
        TwiddlesFFT(n) = GetTwiddle(n, -1)
    Next n
    
    ReDim TwiddlesIFFT(FFTSize - 1)
    For n = 0 To FFTSize - 1
        TwiddlesIFFT(n) = GetTwiddle(n, 1)
    Next n
End Sub


Public Sub FFT(ByRef DataIn() As FPComplex, ByRef DataOut() As FPComplex, ByRef Twiddles() As FPComplex)
    Dim Layer As Long
    Dim BlockCount As Long
    Dim BlockSize As Long
    Dim HalfBlockSize As Long
    Dim BlockNumber As Long
    Dim ButterflyNumber As Long
    Dim n0 As Long
    Dim n1 As Long
    Dim n2 As Long
    Dim DataTemp() As FPComplex
    Dim n As Long
    
    
    DataTemp() = DataIn()
    
    For Layer = 0 To FFTSize_Log2 - 1
        BlockCount = Pow2_LUT(Layer)
        BlockSize = FFTSize \ BlockCount
        HalfBlockSize = BlockSize \ 2
        For BlockNumber = 0 To BlockCount - 1
            n0 = BlockNumber * BlockSize
            For ButterflyNumber = 0 To HalfBlockSize - 1
                n1 = n0 + ButterflyNumber
                n2 = n1 + HalfBlockSize
                
                'do butterfly calculation
                DataOut(n1) = FPCAdd(DataTemp(n1), DataTemp(n2))
                DataOut(n2) = FPCMult(FPCSub(DataTemp(n1), DataTemp(n2)), Twiddles(ButterflyNumber * BlockCount))
                
            Next ButterflyNumber
        Next BlockNumber
        DataTemp() = DataOut()
    Next Layer
    
    For n = 0 To FFTSize - 1
        DataOut(n) = DataTemp(BOR_LUT(n))
    Next n
End Sub


Private Function GetTwiddle(ByVal k As Long, ByVal Sign As Long) As FPComplex
    With GetTwiddle
        .Re = DoubleToFP(Cos(2 * Pi * k / FFTSize))
        .Im = DoubleToFP(Sin(2 * Pi * k / FFTSize) * Sign)
    End With
End Function
```

I've changed SetupFFT to take 2 user supplied arrays to be filled with the twiddles for both the FFT and IFFT processes. I also made the corresponding changes to the FFT function. Whether it performs and FFT or IFFT is now determined by which array of twiddle factors is supplied to its 3rd argument. This way it also no longer takes an argument for Sign like it had previously.

----------


## Ben321

For those of you who don't want to have to use functions to convert between floating point and fixed point, before and after performing an FFT on your data, I've gone back and made the same changes to my floating point version of the function (the one where the Re and Im parts of the Complex user defined type are of type Double instead of type FP1000). Here's the code for the improved floating point FFT module.


```
Public Const Pi As Double = 3.14159265358979

Public Type Complex
    Re As Double
    Im As Double
End Type


Dim FFTSize As Long
Dim FFTSize_Log2 As Long
Dim BOR_LUT() As Long
Dim Pow2_LUT() As Long


Public Sub SetupFFT(ByVal Size As Long, ByRef TwiddlesFFT() As Complex, ByRef TwiddlesIFFT() As Complex)
    Dim n As Long
    Dim n1 As Long
    Dim n2 As Long
    Dim i As Long
    
    FFTSize = Size
    FFTSize_Log2 = Log(FFTSize) / Log(2)
    
    ReDim BOR_LUT(FFTSize - 1)
    
    For n = 0 To FFTSize - 1
        n1 = n
        n2 = 0
        For i = 0 To FFTSize_Log2 - 1
            n2 = (n2 * 2) Or (n1 And 1)
            n1 = n1 \ 2
        Next i
        BOR_LUT(n) = n2
    Next n
    
    ReDim Pow2_LUT(FFTSize_Log2 - 1)
    For n = 0 To FFTSize_Log2 - 1
        Pow2_LUT(n) = 2 ^ n
    Next n
    
    ReDim TwiddlesFFT(FFTSize - 1)
    For n = 0 To FFTSize - 1
        TwiddlesFFT(n) = GetTwiddle(n, -1)
    Next n
    
    ReDim TwiddlesIFFT(FFTSize - 1)
    For n = 0 To FFTSize - 1
        TwiddlesIFFT(n) = GetTwiddle(n, 1)
    Next n
End Sub


Public Sub FFT(ByRef DataIn() As Complex, ByRef DataOut() As Complex, ByRef Twiddles() As Complex)
    Dim Layer As Long
    Dim BlockCount As Long
    Dim BlockSize As Long
    Dim HalfBlockSize As Long
    Dim BlockNumber As Long
    Dim ButterflyNumber As Long
    Dim n0 As Long
    Dim n1 As Long
    Dim n2 As Long
    Dim DataTemp() As Complex
    Dim n As Long
    
    
    
    DataTemp() = DataIn()
    
    For Layer = 0 To FFTSize_Log2 - 1
        BlockCount = Pow2_LUT(Layer)
        BlockSize = FFTSize \ BlockCount
        HalfBlockSize = BlockSize \ 2
        For BlockNumber = 0 To BlockCount - 1
            n0 = BlockNumber * BlockSize
            For ButterflyNumber = 0 To HalfBlockSize - 1
                n1 = n0 + ButterflyNumber
                n2 = n1 + HalfBlockSize
                
                'do butterfly calculation
                DataOut(n1) = ComplexAdd(DataTemp(n1), DataTemp(n2))
                DataOut(n2) = ComplexMult(ComplexSub(DataTemp(n1), DataTemp(n2)), Twiddles(ButterflyNumber * BlockCount))
                
            Next ButterflyNumber
        Next BlockNumber
        DataTemp() = DataOut()
    Next Layer
    
    For n = 0 To FFTSize - 1
        DataOut(n) = DataTemp(BOR_LUT(n))
    Next n
End Sub


Private Function ComplexAdd(ByRef Value1 As Complex, ByRef Value2 As Complex) As Complex
    ComplexAdd.Re = Value1.Re + Value2.Re
    ComplexAdd.Im = Value1.Im + Value2.Im
End Function

Private Function ComplexSub(ByRef Value1 As Complex, ByRef Value2 As Complex) As Complex
    ComplexSub.Re = Value1.Re - Value2.Re
    ComplexSub.Im = Value1.Im - Value2.Im
End Function

Private Function ComplexMult(ByRef Value1 As Complex, ByRef Value2 As Complex) As Complex
    ComplexMult.Re = Value1.Re * Value2.Re - Value1.Im * Value2.Im
    ComplexMult.Im = Value1.Re * Value2.Im + Value1.Im * Value2.Re
End Function


Private Function GetTwiddle(ByVal k As Long, ByVal Sign As Long) As Complex
    With GetTwiddle
        .Re = Cos(2 * Pi * k / FFTSize)
        .Im = Sin(2 * Pi * k / FFTSize) * Sign
    End With
End Function
```

Note that if you want, you can change go to the type definition for Complex and change both Re and Im to type Single if you want. This may give a slight speed up, as 32 bit programs (like those written in VB6) may handle 64bit floating point values (doubles) slower than 32bit floating point values (singles).

----------


## Elroy

> This may give a slight speed up, as 32 bit programs (like those written in VB6) may handle 64bit floating point values (doubles) slower than 32bit floating point values (singles).


Doesn't VB6 serve up 80-bits (and then round upon return) to the FPU for floating point operations, regardless of Single or Double?  If this is the case, then I doubt Single vs Double makes much difference (except for memory requirements).

----------

