Fast Fourier Transform Using Excel's VBA

Question

I'm trying to implement a Fast Fourier Transform (Radix-2) in MS's Excel VBA. The code I'm using pulls data from a range in the worksheet, does the calculations, then dumps the results in the adjacent columns. What I'm having trouble with is 1) know what to do with the resulting X[k] arrays, and 2) matching these results with the results from Excel's built in FFT (they do not currently match). The code is shown below. Thanks in advance for your help.

Sub Enforce_DecimationInTime()

On Error GoTo ERROR_HANDLING
Dim SubName As String
SubName = "Enforce_DecimationInTime()"

Dim WS As Worksheet
Dim n As Long, v As Long, LR As Long, x As Long

Set WS = Worksheets("FFT")
LR = WS.Range("A" & Rows.Count).End(xlUp).Row
n = LR - 1
Do Until 2 ^ x <= n And 2 ^ (x + 1) > n     'locates largest power of 2 from size of input array
    x = x + 1
Loop
n = n - (n - 2 ^ x) 'calculates n using the largest power of 2
If n + 1 <> WS.Range("A" & Rows.Count).End(xlUp).Row Then
    WS.Range("A" & 2 ^ x + 2 & ":A" & LR).Delete xlUp   'deletes extra input data
End If
v = WorksheetFunction.Log(n, 2)     'calculates number of decimations necessary

Application.ScreenUpdating = False
For x = 1 To v
    Call Called_Core.DecimationInTime(WS, n, 2 ^ x, x)  'calls decimation in time subroutine
Next x
Application.ScreenUpdating = True

Exit Sub
ERROR_HANDLING:
    MsgBox "Error encountered in " & SubName & ": exiting subroutine." _
    & vbNewLine _
    & vbNewLine & "Error description: " & Err.Description _
    & vbNewLine & "Error number: " & Err.Number, vbCritical, Title:="Error!"
    End

End Sub

The above subroutine calls the below subroutine through a For/Next loop to the count of "v".

Sub DecimationInTime(WS As Worksheet, n As Long, Factor As Integer, x As Long)

On Error GoTo ERROR_HANDLING
Dim SubName As String
SubName = "DecimationInTime()"

Dim f_1() As Single, f_2() As Single
Dim i As Long, m As Long, k As Long
Dim TFactor_N1 As String, TFactor_N2 As String, X_k() As String
Dim G_1() As Variant, G_2() As Variant

ReDim f_1(0 To n / Factor - 1) As Single
ReDim f_2(0 To n / Factor - 1) As Single
ReDim G_1(0 To n / 1 - 1) As Variant
ReDim G_2(0 To n / 1 - 1) As Variant
ReDim X_k(0 To n - 1) As String

TFactor_N1 = WorksheetFunction.Complex(0, -2 * WorksheetFunction.Pi / (n / 1))  'twiddle factor for N
TFactor_N2 = WorksheetFunction.Complex(0, -2 * WorksheetFunction.Pi / (n / 2))  'twiddle factor for N/2

For i = 0 To n / Factor - 1
    f_1(i) = WS.Range("A" & 2 * i + 2).Value    'assign input data
    f_2(i) = WS.Range("A" & 2 * i + 3).Value    'assign input data
Next i

WS.Cells(1, 1 + x).Value = "X[" & x & "]"   'labels X[k] column with k number
For k = 0 To n / 2 - 1
    For m = 0 To n / Factor - 1
        G_1(m) = WorksheetFunction.ImProduct(WorksheetFunction.ImPower(TFactor_N2, k * m), WorksheetFunction.Complex(f_1(m), 0))    'defines G_1[m]
        G_2(m) = WorksheetFunction.ImProduct(WorksheetFunction.ImPower(TFactor_N2, k * m), WorksheetFunction.Complex(f_2(m), 0))    'defines G_2[m]
    Next m
    X_k(k) = WorksheetFunction.ImSum(WorksheetFunction.ImSum(G_1), WorksheetFunction.ImProduct(WorksheetFunction.ImSum(G_2), WorksheetFunction.ImPower(TFactor_N1, k)))  'defines X[k] for k
    If k <= n / 2 Then X_k(k + n / 2) = WorksheetFunction.ImSum(WorksheetFunction.ImSum(G_1), WorksheetFunction.ImProduct(WorksheetFunction.ImSum(G_2), WorksheetFunction.ImPower(TFactor_N1, k), WorksheetFunction.Complex(-1, 0)))  'defines X[k] for k + n/2
    WS.Cells(k + 2, 1 + x).Value = X_k(k)
    WS.Cells(k + 2 + n / 2, 1 + x).Value = X_k(k + n / 2)
Next k

Exit Sub
ERROR_HANDLING:
    MsgBox "Error encountered in " & SubName & ": exiting subroutine." _
    & vbNewLine _
    & vbNewLine & "Error description: " & Err.Description _
    & vbNewLine & "Error number: " & Err.Number, vbCritical, Title:="Error!"
    End

End Sub

Answer 1

In alternative to the VBA solutions already posted, recent versions of Excel allow to implement the FFT as a pure formula with LAMBDA functions (i.e. without any VBA code).

One such implementation is https://github.com/altomani/XL-FFT.

For power of two length it uses a recursive radix-2 Cooley-Tukey algorithm and for other length a version of Bluestein's algorithm that reduces the calculation to a power of two case.

Answer 2

Implementing FFT in ExcelVBA is kind of involved, but not too bad. For typical applications, the input signal is usually real-valued, not coplex-valued. This would be the case if you were measuring a dynamic signal from a velocity or acceleration transducer, or microphone. Shown here is a DFT that will convert any number of input pairs, (eg. time and velocity). They do not have to be 2^N number of data (required for FFT). Usually the time is evenly divided so that all you need is DeltaTime (the time interval between your data). Let me drop in the code here:

    Sub dft()
    Dim ytime(0 To 18000) As Double 'Time history values such as velocity or acceleration
    Dim omega(0 To 8096) As Double 'Discreet frequency values used in transform
    Dim yfreqr(0 To 8096) As Double 'Real valued component of transform
    Dim yfreqi(0 To 8096) As Double 'Imaginary component of transform
    Dim t As Double, sumr As Double, sumi As Double, sum As Double 'Cumulative sums
    Dim omegadt As Double, omegat As Double, deltime As Double 'More constants self explanitory
    Dim wksInData As Worksheet 'This is the Excel worksheet where the data is read from and written to
    Dim s As Integer, i As Integer 'Counters for the transform loops
    Dim transdim As Integer 'Dimension of the transform
    'Read number of values to read, delta time
    'Read in dimension of transform
    
    Set wksInData = Worksheets("DFT Input") 'This is what I named the worksheet
    numval = wksInData.Cells(5, 2)
    deltime = wksInData.Cells(6, 2)
    transdim = wksInData.Cells(5, 4)
    For i = 0 To numval - 1    'Read in all the input data, its just a long column
     ytime(i) = wksInData.Cells(i + 8, 2) 'So the input starts on row 8 column 2 (time values on column 1 for plotting)
    Next i 'Loop until you have all the numbers you need

   'Start the transform outer loop...for each discreet frequency
   'Value s is the counter from 0 to 1/2 transform dimension
   'So if you have 2000 numbers to convert, transdim is 2000

    For s = 0 To transdim / 2  'Since transform is complex valued, use only 1/2 the number of transdim
     sumr = 0#  'Set the sum of real values to zero
     sumi = 0#  'Set the sum of imaginary values to zero
     omega(s) = 2# * 3.14159265 * s / (transdim * deltime) 'These are the discreet frequencies
     omegadt = omega(s) * deltime  'Just a number used in computations
   ' Start the inner loop for DFT
     For i = 0 To numval - 1
       sumr = sumr + ytime(i) * Cos(omegadt * i) 'This is the real valued sum
       sumi = sumi + ytime(i) * Sin(omegadt * i) 'This is the complex valued sum
     Next i  ' and back for more
     yfreqr(s) = sumr * 2# / transdim  'This is what is called the twiddle factor, just a constant
     yfreqi(s) = -sumi * 2# / transdim  'Imaginary component is negative 
   Next s

   'One last adjustment for the first and last transform values
   'They are only 1/2 of the rest, but it is easiest to do this now after the inner loop is done

   yfreqr(0) = yfreqr(0) / 2#  'Beginning factor
   yfreqi(0) = yfreqi(0) / 2#
   yfreqr(transdim / 2) = yfreqr(transdim / 2) / 2# 'End factor
   yfreqi(transdim / 2) = yfreqi(transdim / 2) / 2#
   wksInData.Cells(2, 8) = "Output" 'Just a column text header
  For s = 0 To transdim / 2   'And write the output to columns 3, 4, 5 to the worksheet
     wksInData.Cells(s + 8, 3) = omega(s)  'remember that magnitude is sqrt(real ^2 + imaginary ^2 ) 
     wksInData.Cells(s + 8, 4) = yfreqr(s) 'but you can do this with an Excel formula on the worksheet
     wksInData.Cells(s + 8, 5) = yfreqi(s) 'same with phase angle = arctan(Imaginary/Real)
   Next s  'End of writeout loop.
   'This is the inverse DFT
   'I like to check my calculation, 
   'Should get the original time series back
   For i = 0 To numval - 1
     sum = 0
     t = deltime * i
     For s = 0 To transdim / 2
       omegat = omega(s) * t
       sum = sum + yfreqr(s) * Cos(omegat) - yfreqi(s) * Sin(omegat)
     Next s
     ytime(i) = sum
  Next i

Answer 3

The function call is not good Call Called_Core.DecimationInTime(WS, n, 2 ^ x, x, TFactor_N1, TFactor_N2)

It should be: Call DecimationInTime(WS, n, 2 ^ x, x, TFactor_N1, TFactor_N2)

Answer 4

I went back through the process and determined my problem was that I had assigned the wrong values to the twiddle factors, TFactor_N1 and TFactor_N2. After fixing this problem and adjusting which values are displayed, I was able to get the same results as Excel's built in FFT. The fixed code is show below.

Sub Enforce_DecimationInTime()

On Error GoTo ERROR_HANDLING
Dim SubName As String
SubName = "Enforce_DecimationInTime()"

Dim WS As Worksheet
Dim n As Long, v As Long, LR As Long, x As Long
Dim TFactor_N1 As String, TFactor_N2 As String

Set WS = Worksheets("FFT")
LR = WS.Range("A" & Rows.Count).End(xlUp).Row
n = LR - 1
Do Until 2 ^ x <= n And 2 ^ (x + 1) > n                                                                     'locates largest power of 2 from size of input array
    x = x + 1
Loop
n = n - (n - 2 ^ x)                                                                                         'calculates n using the largest power of 2
If n + 1 <> WS.Range("A" & Rows.Count).End(xlUp).Row Then
    WS.Range("A" & 2 ^ x + 2 & ":A" & LR).Delete xlUp                                                       'deletes extra input data
End If
v = WorksheetFunction.Log(n, 2)                                                                             'calculates number of decimations necessary

TFactor_N1 = WorksheetFunction.ImExp(WorksheetFunction.Complex(0, -2 * WorksheetFunction.Pi / (n / 1)))     'twiddle factor for N
TFactor_N2 = WorksheetFunction.ImExp(WorksheetFunction.Complex(0, -2 * WorksheetFunction.Pi / (n / 2)))     'twiddle factor for N/2

Application.ScreenUpdating = False
For x = 1 To v
    Call Called_Core.DecimationInTime(WS, n, 2 ^ x, x, TFactor_N1, TFactor_N2)                              'calls decimation in time subroutine
Next x
Application.ScreenUpdating = True

Exit Sub
ERROR_HANDLING:
    MsgBox "Error encountered in " & SubName & ": exiting subroutine." _
    & vbNewLine _
    & vbNewLine & "Error description: " & Err.Description _
    & vbNewLine & "Error number: " & Err.Number, vbCritical, Title:="Error!"
    End

End Sub


Sub DecimationInTime(WS As Worksheet, n As Long, Factor As Integer, x As Long, TFactor_N1 As String, TFactor_N2 As String)

On Error GoTo ERROR_HANDLING
Dim SubName As String
SubName = "DecimationInTime()"

Dim f_1() As String, f_2() As String
Dim i As Long, m As Long, k As Long
Dim X_k() As String
Dim G_1() As Variant, G_2() As Variant

ReDim f_1(0 To n / Factor - 1) As String
ReDim f_2(0 To n / Factor - 1) As String
ReDim G_1(0 To n / 1 - 1) As Variant
ReDim G_2(0 To n / 1 - 1) As Variant
ReDim X_k(0 To n - 1) As String

For i = 0 To n / Factor - 1
    f_1(i) = WS.Cells(2 * i + 2, 1).Value                                                                   'assign input data
    f_2(i) = WS.Cells(2 * i + 3, 1).Value                                                                   'assign input data
Next i
For k = 0 To n / 2 - 1
    For m = 0 To n / Factor - 1                                                                             'defines G_1[m] and G_2[m]
        G_1(m) = WorksheetFunction.ImProduct(WorksheetFunction.ImPower(TFactor_N2, k * m), f_1(m))
        G_2(m) = WorksheetFunction.ImProduct(WorksheetFunction.ImPower(TFactor_N2, k * m), f_2(m))
    Next m                                                                                                  'defines X[k] for k and k + n/2
    X_k(k) = WorksheetFunction.ImSum(WorksheetFunction.ImSum(G_1), WorksheetFunction.ImProduct(WorksheetFunction.ImSum(G_2), WorksheetFunction.ImPower(TFactor_N1, k)))
    If k <= n / 2 Then X_k(k + n / 2) = WorksheetFunction.ImSub(WorksheetFunction.ImSum(G_1), WorksheetFunction.ImProduct(WorksheetFunction.ImSum(G_2), WorksheetFunction.ImPower(TFactor_N1, k)))
    If x = 1 Then
        WS.Cells(k + 2, 1 + x).Value = X_k(k)
        WS.Cells(k + 2 + n / 2, 1 + x).Value = X_k(k + n / 2)
    End If
Next k

Exit Sub
ERROR_HANDLING:
    MsgBox "Error encountered in " & SubName & ": exiting subroutine." _
    & vbNewLine _
    & vbNewLine & "Error description: " & Err.Description _
    & vbNewLine & "Error number: " & Err.Number, vbCritical, Title:="Error!"
    End

End Sub