Attribute VB_Name = "VBFFT"
'--------------------------------------------------------------------
' VB FFT Release 2-B
' by Murphy McCauley (MurphyMc@Concentric.NET)
' 10/01/99
'--------------------------------------------------------------------
' About:
' This code is very, very heavily based on Don Cross's fourier.pas
' Turbo Pascal Unit for calculating the Fast Fourier Transform.
' I've not implemented all of his functions, though I may well do
' so in the future.
' For more info, you can contact me by email, check my website at:
' http://www.fullspectrum.com/deeth/
' or check Don Cross's FFT web page at:
' http://www.intersrv.com/~dcross/fft.html
' You also may be intrested in the FFT.DLL that I put together based
' on Don Cross's FFT C code.  It's callable with Visual Basic and
' includes VB declares.  You can get it from either website.
'--------------------------------------------------------------------
' History of Release 2-B:
' Fixed a couple of errors that resulted from me mucking about with
'   variable names after implementation and not re-checking.  BAD ME.
'  --------
' History of Release 2:
' Added FrequencyOfIndex() which is Don Cross's Index_to_frequency().
' FourierTransform() can now do inverse transforms.
' Added CalcFrequency() which can do a transform for a single
' frequency.
'--------------------------------------------------------------------
' Usage:
' The useful functions are:
' FourierTransform() performs a Fast Fourier Transform on an pair of
'  Double arrays -- one real, one imaginary.  Don't want/need
'  imaginary numbers?  Just use an array of 0s.  This function can
'  also do inverse FFTs.
' FrequencyOfIndex() can tell you what actual frequency a given index
'  corresponds to.
' CalcFrequency() transforms a single frequency.
'--------------------------------------------------------------------
' Notes:
' All arrays must be 0 based (i.e. Dim TheArray(0 To 1023) or
'  Dim TheArray(1023)).
' The number of samples must be a power of two (i.e. 2^x).
' FrequencyOfIndex() and CalcFrequency() haven't been tested much.
' Use this ENTIRELY AT YOUR OWN RISK.
'--------------------------------------------------------------------

'xFX JumpStart (http://www.xfx.net)
'The code has been modified and optimized to be used in conjunction
'with the DXVU Meter control

Option Explicit
Public Const Pi = 3.14159265358979
Public Const PI2 = 2# * Pi

Private Function NumberOfBitsNeeded(PowerOfTwo As Integer) As Byte

    Dim i As Byte

    For i = 0 To 16
        If (PowerOfTwo And (2 ^ i)) <> 0 Then
            NumberOfBitsNeeded = i
            Exit Function
        End If
    Next i
    
End Function

Private Function IsPowerOfTwo(x As Long) As Boolean

    IsPowerOfTwo = Not (x And (x - 1)) And (x >= 2)

End Function

Private Function ReverseBits(ByVal Index As Long, ByVal NumBits As Byte) As Long

    Dim i As Byte
    Dim Rev As Long
    
    For i = 0 To NumBits - 1
        Rev = (Rev * 2) Or (Index And 1)
        Index = Index \ 2
    Next i
    
    ReverseBits = Rev
    
End Function

Public Sub FourierTransform(NumSamples As Integer, RealIn() As Double, RealOut() As Double, DoInverse As Boolean)
    
    Static NumBits As Byte
    Static i As Long, j As Long, k As Long, n As Long
    Static BlockSize As Long, BlockEnd As Long
    
    Static DeltaAngle As Double, DeltaAr As Double
    Static Alpha As Double, Beta As Double
    Static TR As Double, TI As Double, AR As Double, AI As Double
    Static ImagOut() As Double
    ReDim ImagOut(0 To NumSamples - 1)
    
    Static rBits() As Long
    Static LastNumSamples
    Static nIter As Long
    
    ' This is the original code to calculate the reversed bits
    'NumBits = NumberOfBitsNeeded(NumSamples)
    'For i = 0 To NumSamples - 1
    '    j = ReverseBits(i, NumBits)
    '    RealOut(j) = CSng(RealIn(i))
    '    ImagOut(j) = 0
    'Next i

    ' This is the enhanced version of the code above that recalculates
    ' the reversed bits only when the NumSamples changes
    ' This saves a couple of milliseconds on the FFT calculation
    If LastNumSamples <> NumSamples Then
        LastNumSamples = NumSamples
        ReDim rBits(0 To NumSamples - 1)
        NumBits = NumberOfBitsNeeded(NumSamples)
        For i = 0 To NumSamples - 1
            rBits(i) = ReverseBits(i, NumBits)
            RealOut(rBits(i)) = CSng(RealIn(i))
            ImagOut(rBits(i)) = 0
        Next i
    Else
        For i = 0 To NumSamples - 1
            RealOut(rBits(i)) = CSng(RealIn(i))
            ImagOut(rBits(i)) = 0
        Next i
    End If
    
    BlockEnd = 1
    BlockSize = 2
    
    Do
        DeltaAngle = PI2 / BlockSize * IIf(DoInverse, -1#, 1#)
        Alpha = 2# * Sin(0.5 * DeltaAngle) ^ 2
        Beta = Sin(DeltaAngle)
        
        For i = 0 To NumSamples - 1 Step BlockSize
            AR = 1#
            AI = 0#
            
            j = i
            For n = 0 To BlockEnd - 1
                k = j + BlockEnd
                TR = AR * RealOut(k) - AI * ImagOut(k)
                TI = AI * RealOut(k) + AR * ImagOut(k)
                RealOut(k) = RealOut(j) - TR
                ImagOut(k) = ImagOut(j) - TI
                RealOut(j) = RealOut(j) + TR
                ImagOut(j) = ImagOut(j) + TI
                DeltaAr = Alpha * AR + Beta * AI
                AI = AI - (Alpha * AI - Beta * AR)
                AR = AR - DeltaAr
                j = j + 1
            Next n
        Next i
        
        BlockEnd = BlockSize
        BlockSize = BlockSize * 2
    Loop While BlockSize <= NumSamples
    
    If DoInverse Then
        For i = 0 To NumSamples - 1
            RealOut(i) = RealOut(i) / NumSamples
        Next i
    End If
    
End Sub

Public Function FrequencyOfIndex(NumberOfSamples As Integer, Index As Integer) As Double
    
    If Index >= NumberOfSamples Then
        FrequencyOfIndex = 0#
        Exit Function
    ElseIf Index <= NumberOfSamples / 2 Then
        FrequencyOfIndex = CDbl(Index) / CDbl(NumberOfSamples)
        Exit Function
    Else
        FrequencyOfIndex = -CDbl(NumberOfSamples - Index) / CDbl(NumberOfSamples)
        Exit Function
    End If

End Function

Public Sub CalcFrequency(NumberOfSamples As Integer, FrequencyIndex As Integer, RealIn() As Double, ImagIn() As Double, RealOut As Double, ImagOut As Double)
    
    Dim k As Long
    Dim Cos1 As Double, Cos2 As Double, Cos3 As Double, Theta As Double, Beta As Double
    Dim Sin1 As Double, Sin2 As Double, Sin3 As Double
    
    Theta = 2 * Pi * FrequencyIndex / CDbl(NumberOfSamples)
    Sin1 = Sin(-2 * Theta)
    Sin2 = Sin(-Theta)
    Cos1 = Cos(-2 * Theta)
    Cos2 = Cos(-Theta)
    Beta = 2 * Cos2
    
    For k = 0 To NumberOfSamples - 2
        'Update trig values
        Sin3 = Beta * Sin2 - Sin1
        Sin1 = Sin2
        Sin2 = Sin3
        
        Cos3 = Beta * Cos2 - Cos1
        Cos1 = Cos2
        Cos2 = Cos3
        
        RealOut = RealOut + RealIn(k) * Cos3 - ImagIn(k) * Sin3
        ImagOut = ImagOut + ImagIn(k) * Cos3 + RealIn(k) * Sin3
    Next
End Sub
