Coding FFT using Excel VBA: no limitation of data size and O(nlogn) performance

The FFT Excel program is updated to v2.0.

You can download it on the bottom of this page or the new page. 
The new version performs FFT for 64 kBytes data in 1 second. (235ms for 16 kBytes data)

 

FFT(Fast Fourier Transformation) can be done in the Excel's internal function. However it has limitation of data size 4096 and it works a little slow.

Therefore, I decided to make FFT code with Excel VBA.

   - created excel file is attached in the end of this article.  

   - FFT running time is under 1 seccond for 16K data size.

 

 

Algorithm

Cooley–Tukey algorithm is used for FFT implementation.

Below is modified algorithm equation for Excel VBA coding:

Xm+1[i]            =Xm[i]+Xm[i+half]W[r]

Xm+1[i+half]=Xm[i]-Xm[i+half]W[r]

 

Usage

1. Call macro "LoadFFTForm" : View > Macros > View Macros > LoadFFTForm

2. Select input data range for "Input Range" and output cell for "Output Range" 

3. Click "Run" button

 

 

 

VBA Code

'Author: HJ Park
'Date  : 2019.5.18
Option Explicit
Public Const myPI As Double = 3.14159265358979
Public Function Log2(X As Long) As Double
  Log2 = Log(X) / Log(2)
End Function
Public Function Ceiling(ByVal X As Double, Optional ByVal Factor As Double = 1) As Double
    ' X is the value you want to round
    ' is the multiple to which you want to round
    Ceiling = (Int(X / Factor) - (X / Factor - Int(X / Factor) > 0)) * Factor
End Function
Public Function Floor(ByVal X As Double, Optional ByVal Factor As Double = 1) As Double
    ' X is the value you want to round
    ' is the multiple to which you want to round
    Floor = Int(X / Factor) * Factor
End Function
Public Function IsPowerOfTwo(N As Long) As Boolean
  If N = 0 Then GoTo EXIT_FUNCTION
  'If Application.WorksheetFunction.Ceiling_Math(Log2(N)) = Application.WorksheetFunction.Floor_Math(Log2(N)) Then
  If Ceiling(Log2(N)) = Floor(Log2(N)) Then
    IsPowerOfTwo = True
  Else
    IsPowerOfTwo = False
  End If
  Exit Function
EXIT_FUNCTION:
  IsPowerOfTwo = False
End Function
' Dec2Bin(13,4) -> {1,1,0,1} -> {1,0,1,1}
Public Function Dec2Bin_Inverse(ByVal num As Long, numOfBits As Integer) As Byte()
  Dim arr() As Byte
  Dim i As Integer, j As Integer
  Dim K As Long
  If (2 ^ numOfBits) <= num Then GoTo EXIT_FUNCTION
  ReDim arr(0 To numOfBits - 1) As Byte
  For i = 0 To numOfBits - 1
    j = (numOfBits - 1) - i
    K = WorksheetFunction.RoundDown(num / (2 ^ j), 0)
    'arr(i) = (k And 1)
    arr(j) = (K And 1)
  Next i
  Dec2Bin_Inverse = arr
EXIT_FUNCTION:
End Function
'{1,0,1,1} -> 1*2^3 + 0*2^2 + 1*2&1 + 1 = 11
Public Function Bin2Dec(arr() As Byte, numOfBits As Integer) As Long
  Dim i As Integer
  Bin2Dec = 0
  For i = 0 To numOfBits - 1
    Bin2Dec = Bin2Dec + arr(i) * 2 ^ (numOfBits - 1 - i)
  Next i
End Function
''''''''''''''''''''''''''''
'''''''''''''''''''''''''''''
' read the first row value in the range and return it as complex value array
Function Range2Array(r As Range) As Complex()
  Dim i As Long, size As Long
  Dim arr() As Complex
  size = r.Rows.Count
  ReDim arr(0 To size - 1) As Complex
  Dim re As Double, im As Double
  For i = 1 To size
    re = WorksheetFunction.ImReal(r.Rows(i).Value)
    im = WorksheetFunction.imaginary(r.Rows(i).Value)
    arr(i - 1) = Cplx(re, im)
  Next i
  Range2Array = arr
End Function
Function SortedNum(num As Long, numOfBits As Integer) As Long
  Dim arr() As Byte
  ' Decimal number -> reversed binary number -> decimal number
  arr = Dec2Bin_Inverse(num, numOfBits)
  SortedNum = Bin2Dec(arr, numOfBits)
End Function
' rangeArr[1 to n, 1]
Function sortToFFTArray(arr() As Complex, size As Long, numOfBits As Integer) As Complex()
  Dim i As Long, j As Long
  Dim sortedArr() As Complex
  ReDim sortedArr(0 To size - 1) As Complex
  For i = 0 To size - 1
    j = SortedNum(i, numOfBits)  '{000,001,010, 011, 100, 101, 110, 111} -> {0, 4, 2, 6, 1, 5, 3, 7}
    sortedArr(j) = arr(i)
  Next i
  sortToFFTArray = sortedArr
End Function
' calculate convolution ring W
' W[r] = cos(theta) + isin(theta)
'   theta = (2pi/N) * r
Function CalculateW(cnt As Long, isInverse As Boolean) As Complex()
  Dim arr() As Complex
  Dim i As Long
  Dim T As Double, theta As Double
  ReDim arr(0 To cnt) As Complex
  T = 2 * myPI / CDbl(cnt)
  If isInverse Then
    For i = 0 To cnt - 1
      theta = -(T * i)
      arr(i) = Cplx(Cos(theta), -Sin(theta))
    Next i
  Else
    For i = 0 To cnt - 1
      theta = T * i
      arr(i) = Cplx(Cos(theta), -Sin(theta))
    Next i
  End If
  CalculateW = arr
End Function
' X({0,1}, [0,n-1]): 2d array.  (0, n) <--> (1,n)
' src: src index of the array. 0 or 1
' tgt: tgt index of the array. 1 or 0
' s : starting index of the data in the array
' size: region size to be calculated
' rGap : r's jumping value
' W(0 ~ n-1) : Convolution ring
Sub RegionFFT(X() As Complex, src As Integer, tgt As Integer, _
            s As Long, size As Long, rGap As Long, W() As Complex)
  Dim i As Long, e As Long
  Dim half As Long
  Dim r As Long
  Dim K As Complex
  ' Xm+1[i] = Xm[i] + Xm[i+half]W[r]
  ' Xm+1[i+half] = Xm[i] - Xm[i]W[r]
  r = 0
  e = s + (size / 2) - 1
  half = size / 2
  For i = s To e
    K = CMult(X(src, i + half), W(r))
    X(tgt, i) = CAdd(X(src, i), K)
    X(tgt, i + half) = CSub(X(src, i), K)
    r = r + rGap
  Next i
End Sub
Sub WriteToTarget(tgtRange As Range, X() As Complex, tgtIdx As Integer, N As Long)
  Dim i As Long
  Dim arr() As Variant
  ReDim arr(0 To N - 1) As Variant
  For i = 0 To N - 1
    'If X(tgtIdx, i).re = 0 And X(tgtIdx, i).im = 0 Then
    '  arr(i) = 0
    'Else
      arr(i) = WorksheetFunction.Complex(X(tgtIdx, i).re, X(tgtIdx, i).im, "j")
    'End If
  Next i
  tgtRange.Rows = Application.Transpose(arr)
End Sub
Public Sub FFT(xRange As Range, tgtRange As Range, isInverse As Boolean)
  Dim i As Long, N As Long
  Dim totalLoop As Integer, curLoop As Integer
  Dim xArr() As Complex, xSortedArr() As Complex
  Dim W() As Complex
  Dim X() As Complex
  '1) check whether 2^r count data
  N = xRange.Rows.Count
  If IsPowerOfTwo(N) = False Then GoTo ERROR_HANDLING
  '2) calculate loop count for FFT
  '   2->1  4->2  8->3
  totalLoop = Log2(N)
  '3) sort x for 2's FFT : converted to reversed binary and then convert to decimal
  'xArr = xRange.Value  'xArr[1,n]
  xArr = Range2Array(xRange)  'xArr[0,n-1]
  xSortedArr = sortToFFTArray(xArr, N, totalLoop) 'xSortedArr[0,n-1]
  '4) calculate W
  W = CalculateW(N, isInverse)
  '5) convert real number to complex number
  ReDim X(0 To 1, 0 To N - 1) As Complex
  For i = 0 To N - 1
    X(0, i) = xSortedArr(i)
  Next i
  '6) Do 2's FFT with sorted x
  Dim srcIdx As Integer, tgtIdx As Integer
  Dim rGap As Long, regionSize As Long
  tgtIdx = 0
  For curLoop = 0 To totalLoop - 1
    tgtIdx = (tgtIdx + 1) Mod 2
    srcIdx = (tgtIdx + 1) Mod 2
    regionSize = 2 ^ (curLoop + 1)
    rGap = 2 ^ (totalLoop - curLoop - 1)
    i = 0
    Do While i < N
      Call RegionFFT(X, srcIdx, tgtIdx, i, regionSize, rGap, W)
      i = i + regionSize
    Loop
  Next curLoop
  '7) Write the result into the tgt region
  Dim resultIdx As Integer
  If (totalLoop Mod 2) = 0 Then resultIdx = 0 Else resultIdx = 1
  If isInverse Then
    For i = 0 To N - 1
      X(resultIdx, i) = CDivR(X(resultIdx, i), N)
    Next i
  End If
  'make zero if the imaginary value is too low < 10^(-15)
  For i = 0 To N - 1
    If X(resultIdx, i).re > -(10 ^ (-10)) And X(resultIdx, i).re < (10 ^ (-10)) Then X(resultIdx, i).re = 0
    If X(resultIdx, i).im > -(10 ^ (-10)) And X(resultIdx, i).im < (10 ^ (-10)) Then X(resultIdx, i).im = 0
  Next i
  Call WriteToTarget(tgtRange, X, resultIdx, N)
  Exit Sub
ERROR_HANDLING:
  MsgBox "should be power of 2 data"
End Sub
Sub LoadFFTForm()
  FFT_Form.Show
End Sub

 

 

' module Name: complexNum

Option Explicit
Public Type Complex
    re As Double
    im As Double
End Type
Public Function Cplx(ByVal real As Double, ByVal imaginary As Double) As Complex
    Cplx.re = real
    Cplx.im = imaginary
End Function
Public Function CAdd(ByRef z1 As Complex, ByRef z2 As Complex) As Complex
    CAdd.re = z1.re + z2.re
    CAdd.im = z1.im + z2.im
End Function
Public Function CSub(ByRef z1 As Complex, ByRef z2 As Complex) As Complex
    CSub.re = z1.re - z2.re
    CSub.im = z1.im - z2.im
End Function
Public Function CMult(ByRef z1 As Complex, ByRef z2 As Complex) As Complex
    CMult.re = z1.re * z2.re - z1.im * z2.im
    CMult.im = z1.re * z2.im + z1.im * z2.re
End Function
Public Function CDivR(ByRef z As Complex, ByVal r As Double) As Complex
' Divide complex number by real number
    CDivR.re = z.re / r
    CDivR.im = z.im / r
End Function

 

' FFT_Form code

Option Explicit
Private Sub CloseButton_Click()
  Unload Me
End Sub
Private Sub RunButton_Click()
  Dim inputRange As Range, outputRange As Range
  ErrorLabel.Caption = ""
  SuccessLabel.Caption = ""
  If InputEdit.Value = "" Then
    ErrorLabel.Caption = "Please choose input data range"
    Exit Sub
  End If
  If OutputEdit.Value = "" Then
    ErrorLabel.Caption = "Please choose output data range"
    Exit Sub
  End If
  On Error GoTo ERROR_HANDLE
  Set inputRange = Range(InputEdit.Value)
  Set outputRange = Range(OutputEdit.Value)
  'check output range
  If outputRange.Rows.Count <> inputRange.Rows.Count Then
    Set outputRange = Range(Cells(outputRange.Row, outputRange.Column), Cells(outputRange.Row + inputRange.Rows.Count - 1, outputRange.Column))
  End If
  Call FFT(inputRange, outputRange, InverseFFTCheck.Value)
  SuccessLabel.Caption = "FFT work is done"
  Exit Sub
ERROR_HANDLE:
  ErrorLabel.Caption = "Choose right data"
End Sub
Private Sub UserForm_Initialize()
  Dim inputRange As Range, outputRange As Range
  'input: selected range
  Set inputRange = Selection
  InputEdit.Value = inputRange.Address
  'output: shifted to one column right from the input range
  Set outputRange = inputRange.Rows(1).Offset(0, 1)
  OutputEdit.Value = outputRange.Address
  'set focus to input edit box
  InputEdit.SetFocus
End Sub

 

 

Excel File

 

fft_v1.1.xlsm

5.12MB

---

 

Please download the new version of below.

The new version of fft_v2.0 file is protected with password.

The password is: 2.718282

 

fft_v2.0 - distribution - encrypted.xlsm

2.91MB

-end-

 

fft_v1.1.xlsm

5.12MB