Attribute VB_Name = "VBFFT"
'--------------------------------------------------------------------
' VB FFT Release 2-B
' by Shearzheng (shearzheng@126.com)
' 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.

' 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.
'--------------------------------------------------------------------

Option Explicit
Const Pi = 3.14159265358979

Function NumberOfBitsNeeded(PowerOfTwo As Long) 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
End Function


Function IsPowerOfTwo(X As Long) As Boolean
    If (X < 2) Then IsPowerOfTwo = False: Exit Function
    If (X And (X - 1)) = False Then IsPowerOfTwo = True
End Function


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


Sub FourierTransform(NumSamples As Long, RealIn() As Double, ImageIn() As Double, RealOut() As Double, ImagOut() As Double, Optional InverseTransform As Boolean = False)
    Dim AngleNumerator As Double
    Dim NumBits As Byte, I As Long, j As Long, K As Long, n As Long, BlockSize As Long, BlockEnd As Long
    Dim DeltaAngle As Double, DeltaAr As Double
    Dim Alpha As Double, Beta As Double
    Dim TR As Double, TI As Double, AR As Double, AI As Double
    
    If InverseTransform Then
        AngleNumerator = -2# * Pi
    Else
        AngleNumerator = 2# * Pi
    End If

    If (IsPowerOfTwo(NumSamples) = False) Or (NumSamples < 2) Then
        Call MsgBox("Error in procedure Fourier:" + vbCrLf + " NumSamples is " + CStr(NumSamples) + ", which is not a positive integer power of two.", , "Error!")
        Exit Sub
    End If
   
    NumBits = NumberOfBitsNeeded(NumSamples)
    For I = 0 To (NumSamples - 1)
        j = ReverseBits(I, NumBits)
        RealOut(j) = RealIn(I)
        ImagOut(j) = ImageIn(I)
    Next
    
    BlockEnd = 1
    BlockSize = 2
    
    Do While BlockSize <= NumSamples
        DeltaAngle = AngleNumerator / BlockSize
        Alpha = Sin(0.5 * DeltaAngle)
        Alpha = 2# * Alpha * Alpha
        Beta = Sin(DeltaAngle)
        
        I = 0
        Do While I < NumSamples
            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
            
            I = I + BlockSize
        Loop
        
        BlockEnd = BlockSize
        BlockSize = BlockSize * 2
    Loop

    If InverseTransform Then
        'Normalize the resulting time samples...
        For I = 0 To NumSamples - 1
            RealOut(I) = RealOut(I) / NumSamples
            ImagOut(I) = ImagOut(I) / NumSamples
        Next
    End If
End Sub


Function FrequencyOfIndex(NumberOfSamples As Long, ByVal Index As Long) As Double
    'Based on IndexToFrequency().  This name makes more sense to me.
    
    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


Sub CalcFrequency(NumberOfSamples As Long, FrequencyIndex As Long, 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