본문 바로가기

푸리에 변환, 신호/푸리에 변환의 모든 것

[07-4 Code] Fastest FFT code for Excel VBA

반응형

This page explain FFT program made for Excel VBA. It is no limit of input data size, therefore you can calculate over 4096 size data and it will perform very fast. 

 

FFT Execution time in the Excel by self-made FFT program 

 

  • Computer Specification: 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz
Input data count FFT execution time [ms]
4 kBytes(4096) 62
8 kBytes 115
16 kBytes 235
64 kBytes 984
  • I think, it's almost 100 times faster than the internal FFT function in the Excel.
  • The program made by pure VBA code not using a DLL library. 
    (Only one DLL is used to measure time for program execution, but it's not related on the FFT code)

Guide for calculating FFT and IFFT in the Excel by the program

1. Select y(n) from B2 to B8193 and press "Launch Dialog for FFT/IFFT" button. Press "Run" button.
    You will obtain Y(k) data.
    Press "Close" button and check whether the Y(k) data is generated well.

  8 kBytes data is calculated in 125ms.

 

2. Now, we will calculate y(n) form the Y(k) using IFFT functionality.
    Select the range of C2 to C8193 and press "Launch Dialog".
    In the dialog, check "Inverse FFT" and press "Run" button

  • y(n) data is generated from the Y(k) by IFFT and the data is same with original y(n) data. 
  • There is some gap between the original y(n) and calculated y(n). It is caused from the round function. 
  • If you raise up the "Round Num Digit" value, the gap will be getting smaller.

 


Source code

'Module: fftProgram
'Author: HJ Park
'Date  : 2019.5.18(v1.0), 2022.8.1(v2.0)

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
    ' Factor 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
    ' Factor is the multiple to which you want to round
    Floor = Int(X / Factor) * Factor
End Function

' return 0 if N is 2^n value,
' return (2^n - N) if N is not 2^n value. 2^n is Ceiling value.
' return -1, if error
Public Function IsPowerOfTwo(N As Long) As Long
  If N = 0 Then GoTo EXIT_FUNCTION
  
  Dim c As Long, F As Double
  c = Ceiling(Log2(N)) 'Factor=0, therefore C is an integer number
  F = Floor(Log2(N))
  
  If c = F Then
    IsPowerOfTwo = 0
  Else
    IsPowerOfTwo = (2 ^ c - N)
  End If
  Exit Function
  
EXIT_FUNCTION:
  IsPowerOfTwo = -1
End Function




''''''''''''''''''''''''''''
'''''''''''''''''''''''''''''

Function MakePowerOfTwoSize(ByRef r As Range, ByVal fillCount As Long) As Boolean
  Dim arr() As Integer
  On Error GoTo ERROR_HANDLE
  
  '1)make a array with zero
  ReDim arr(0 To fillCount - 1) As Integer
    
  '2)set a range to be filled with zero
  Dim fillRowStart As Long
  Dim fillRange As Range
  
  fillRowStart = r.Row + r.Rows.Count
  Set fillRange = Range(Cells(fillRowStart, r.Column), Cells(fillRowStart + fillCount - 1, r.Column))
  
  '3)fill as zero
  fillRange = arr
  
  '4)update range area to be extended
  Set r = Union(r, fillRange)
  
  MakePowerOfTwoSize = True
  Exit Function
  
ERROR_HANDLE:
  MakePowerOfTwoSize = False
End Function

' read 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
  
  On Error GoTo ERROR_HANDLE
  For i = 1 To size
    arr(i - 1) = String2Complex(r.Rows(i).Value)
  Next i
  
  Range2Array = arr
  
  Exit Function
ERROR_HANDLE:
  MsgBox "Error: " & i
End Function

Function ArrangedNum(num As Long, numOfBits As Integer) As Long
  Dim arr() As Byte
  Dim i As Integer, j As Integer
  Dim k As Long
  
  If (2 ^ numOfBits) <= num Then GoTo EXIT_FUNCTION
  
  '1) Decimal number -> Reversed Binary array : (13,4) -> {1,1,0,1} -> {1,0,1,1}
  ReDim arr(0 To numOfBits - 1) As Byte
  For i = 0 To numOfBits - 1
    j = (numOfBits - 1) - i
    k = Int((num / (2 ^ j)))
    arr(j) = (k And 1)
  Next i
  
  '2) Reversed Binary -> Decimal: {1,0,1,1} -> 1*2^3 + 0*2^2 + 1*2&1 + 1 = 11
  Dim d As Long
  For i = 0 To numOfBits - 1
    d = d + (arr(i) * 2 ^ (numOfBits - 1 - i))
  Next i
  
  ArrangedNum = d
  Exit Function
  
EXIT_FUNCTION:
  ArrangedNum = 0
End Function

' rangeArr[1 to n, 1]
Function arrangeToFFTArray(arr() As Complex, size As Long, numOfBits As Integer) As Complex()
  Dim i As Long, j As Long
  Dim arrangedArr() As Complex
    
  ReDim arrangedArr(0 To size - 1) As Complex
  For i = 0 To size - 1
    j = ArrangedNum(i, numOfBits)  '{000,001,010, 011, 100, 101, 110, 111} -> {0, 4, 2, 6, 1, 5, 3, 7}
    arrangedArr(j) = arr(i)
  Next i
  
  arrangeToFFTArray = arrangedArr
End Function

' calculate convolution ring W
' W[k] = cos(theta) - isin(theta)
'   theta = (2pi*k/N)
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
  Dim N As Long, N2 As Long
  
  N = cnt
  N2 = N / 2
  ReDim arr(0 To N2 - 1) As Complex  'enough to calculate 0 to (N/2 -1)
  T = 2 * myPI / CDbl(N)
  
  If isInverse Then
    For i = 0 To N2 - 1
      theta = -(T * i)
      arr(i) = Cplx(Cos(theta), -Sin(theta))
    Next i
  Else
    For i = 0 To N2 - 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
' kJump : k'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, kJump As Long, W() As Complex)
  Dim i As Long, e As Long
  Dim half As Long
  Dim k As Long
  Dim T As Complex
  
  ' Xm+1[i] = Xm[i] + Xm[i+half]W[k]
  ' Xm+1[i+half] = Xm[i] - Xm[i+half]W[k]
  k = 0
  e = s + (size / 2) - 1
  half = size / 2
  For i = s To e
    T = CMult(X(src, i + half), W(k))
    X(tgt, i) = CAdd(X(src, i), T)
    X(tgt, i + half) = CSub(X(src, i), T)
    k = k + kJump
  Next i
End Sub

Sub WriteToTarget(tgtRange As Range, X() As Complex, tgtIdx As Integer, N As Long, roundDigit As Integer)
  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).im < 0 Then
     arr(i) = Round(X(tgtIdx, i).re, roundDigit) & Round(X(tgtIdx, i).im, roundDigit) & "i"
   Else
     arr(i) = Round(X(tgtIdx, i).re, roundDigit) & "+" & Round(X(tgtIdx, i).im, roundDigit) & "i"
   End If
  Next i
  
  tgtRange.Rows = Application.Transpose(arr)
End Sub

' xRange: input data
' tgtRange: output range
' isInverse: FFT or IFFT
Public Function FFT_Forward(xRange As Range, tgtRangeStart As Range, roundDigit As Integer, isInverse As Boolean) As Complex()
  Dim i As Long, N As Long
  Dim totalLoop As Integer, curLoop As Integer 'enough as Integer b/c it is used for loop varoable
  Dim xArr() As Complex, xSortedArr() As Complex
  Dim W() As Complex 'convolution ring
  Dim X() As Complex 'output result
  Dim errMsg As String
  
  errMsg = "Uncatched error"
    
  '1) check whether 2^r count data, if not pad to zero
  Dim fillCount As Long
  N = xRange.Rows.Count
  fillCount = IsPowerOfTwo(N)
  If fillCount = -1 Then
    errMsg = "No input data. Choose input data"
    GoTo ERROR_HANDLE
  End If
  If fillCount <> 0 Then
    If MakePowerOfTwoSize(xRange, fillCount) = False Then  'xRange's size will be chnaged
      errMsg = "Error while zero padding"
      GoTo ERROR_HANDLE
    End If
  End If
  
  '2) calculate loop count for FFT: 2->1  4->2  8->3 ...
  N = xRange.Rows.Count 'xRange's size can be changed so read one more...
  totalLoop = Log2(N)
  
  '3) sort x for 2's FFT : convert to reversed binary and then convert to decimal
  xArr = Range2Array(xRange)  'xArr[0,n-1]
  xSortedArr = arrangeToFFTArray(xArr, N, totalLoop) 'xSortedArr[0,n-1]
  
  '4) calculate W
  W = CalculateW(N, isInverse)
  
  '5) use 2-dimensional array to save memory space. X[0, ] <-> X[1, ]
  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 kJump 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)          ' if N=8: 2 -> 4 -> 8
    kJump = 2 ^ (totalLoop - curLoop - 1)    ' if N=8: 4 -> 2 -> 1
    i = 0
    Do While i < N
      Call RegionFFT(X, srcIdx, tgtIdx, i, regionSize, kJump, W)
      i = i + regionSize
    Loop
  Next curLoop
  
  '7)return the value
  Dim resultIdx As Integer
  If (totalLoop Mod 2) = 0 Then resultIdx = 0 Else resultIdx = 1
  
  Dim result() As Complex
  ReDim result(0 To N - 1) As Complex
  If isInverse = True Then
    For i = 0 To N - 1
      result(i) = CDivR(X(resultIdx, i), N)
    Next i
  Else
    For i = 0 To N - 1
      result(i) = X(resultIdx, i)
    Next i
  End If
  
  FFT_Forward = result
  
  Exit Function
  
ERROR_HANDLE:
  Err.Raise Number:=vbObjectError, Description:=("FFT calculation error: " & errMsg)
End Function


Public Sub FFT(xRange As Range, tgtRangeStart As Range, roundDigit As Integer)
  Dim X() As Complex
  Dim tgtRange As Range
  
  '1. calculate FFT_forward value
  On Error GoTo ERROR_HANDLE
  X = FFT_Forward(xRange, tgtRangeStart, roundDigit, False)
  
  '2. write to the worksheet
  Dim N As Long
  N = UBound(X) - LBound(X) + 1
  
  Dim i As Long
  Dim arr() As Variant
  ReDim arr(0 To N - 1) As Variant
  For i = 0 To N - 1
   If X(i).im < 0 Then
     arr(i) = Round(X(i).re, roundDigit) & Round(X(i).im, roundDigit) & "i"
   Else
     arr(i) = Round(X(i).re, roundDigit) & "+" & Round(X(i).im, roundDigit) & "i"
   End If
  Next i
  
  Set tgtRange = Range(Cells(tgtRangeStart.Row, tgtRangeStart.Column), Cells(tgtRangeStart.Row + N - 1, tgtRangeStart.Column))
  tgtRange.Rows = Application.Transpose(arr)
  Exit Sub
  
ERROR_HANDLE:
  
End Sub

Public Sub IFFT(xRange As Range, tgtRangeStart As Range, roundDigit As Integer)
  Dim X() As Complex
  Dim tgtRange As Range
  
  '1. calculate FFT_forward value
  On Error GoTo ERROR_HANDLE
  X = FFT_Forward(xRange, tgtRangeStart, roundDigit, True)
  
  '2.write to the worksheet
  Dim N As Long
  N = UBound(X) - LBound(X) + 1
  
  Dim arr() As Variant
  ReDim arr(0 To N - 1) As Variant
  Dim i As Long
  For i = 0 To N - 1
    arr(i) = Round(X(i).re, roundDigit)
  Next i
  
  Set tgtRange = Range(Cells(tgtRangeStart.Row, tgtRangeStart.Column), Cells(tgtRangeStart.Row + N - 1, tgtRangeStart.Column))
  tgtRange.Rows = Application.Transpose(arr)
  Exit Sub
  
ERROR_HANDLE:

End Sub


Sub LoadFFTForm()
  FFT_Form.Show
End Sub

 

'Module: 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

Public Function CMagnitude(ByRef z As Complex) As Double
    CMagnitude = Sqr((z.re ^ 2) + (z.im ^ 2))
End Function

Public Function String2Complex(ByVal s As String) As Complex
  Dim sLen As Integer
  Dim re As Double, im As Double
  
  s = Trim(s)
  sLen = Len(s)
  
  If sLen = 0 Then: GoTo ERROR_HANDLE
  
  'support a+aj style
  If InStr(s, "j") Then
    s = Replace(s, "j", "i")
  End If
  
  'total cases: (1){a, -a}  (2){-a+bi  +a+bi} {-a-bi  +a-bi} {-bi  +bi} (3){a+bi, a-bi, bi}
  Dim pos As Integer
  If InStr(s, "i") = 0 Then ' (1){a, -a}
    re = CDbl(s)
    im = 0
  ElseIf InStr(s, "+") = 1 Or InStr(s, "-") = 1 Then '(2){-a+bi  +a+bi} {-a-bi  +a-bi} {-bi  +bi}
    pos = InStr(2, s, "+")
    If pos > 0 Then '-a+bi  +a+bi
      re = CDbl(Mid(s, 1, pos - 1))
      im = CDbl(Mid(s, pos + 1, sLen - pos - 1))
    Else
      pos = InStr(2, s, "-")
      If pos > 0 Then '-a-bi  +a-bi
        pos = InStr(2, s, "-")
        re = CDbl(Mid(s, 1, pos - 1))
        im = CDbl(Mid(s, pos, sLen - pos))
      Else ' -bi  +bi
        re = 0
        im = CDbl(Left(s, sLen - 1))
      End If
    End If
  Else '(3){a+bi a-bi, bi}
    pos = InStr(s, "+")
    If pos > 0 Then 'a+bi
      re = CDbl(Mid(s, 1, pos - 1))
      im = CDbl(Mid(s, pos + 1, sLen - pos - 1))
    Else 'a-bi  bi
      pos = InStr(s, "-")
      If pos > 0 Then 'a-bi
        re = CDbl(Mid(s, 1, pos - 1))
        im = CDbl(Mid(s, pos, sLen - pos))
      Else 'bi
        re = 0
        im = CDbl(Left(s, sLen - 1))
      End If
    End If
  End If
      
  String2Complex = Cplx(re, im)
  Exit Function
  
ERROR_HANDLE:
  String2Complex = Cplx(0, 0)
End Function

 

' "FFT_Form" Code

Option Explicit

Private mStartTime As Long
Private Declare PtrSafe Function GetTickCount Lib "kernel32" () As Long


Public Sub StartTimer()
    mStartTime = GetTickCount
End Sub

Public Function EndTimer() As Long
    EndTimer = (GetTickCount - mStartTime)
End Function

' When close button is pressed.
Private Sub CloseButton_Click()
  Unload Me
End Sub

Private Sub RoundSpin_SpinUp()
  If Me.RoundTextBox = "" Then Me.RoundTextBox = 0
  
  If RoundTextBox >= RoundSpin.Max Then
    RoundTextBox = RoundSpin.Max
  Else
    Me.RoundTextBox = Me.RoundTextBox + 1
  End If
End Sub


Private Sub RoundSpin_SpinDown()
  If Me.RoundTextBox = "" Then Me.RoundTextBox = 0
  
  If RoundTextBox <= RoundSpin.Min Then
    RoundTextBox = RoundSpin.Min
  Else
    Me.RoundTextBox = Me.RoundTextBox - 1
  End If
End Sub



' Run button event
Private Sub RunButton_Click()
  Dim inputRange As Range, outputRangeStart 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 a start cell for output data range"
    Exit Sub
  End If
  
  If RoundTextBox > RoundSpin.Max Or RoundTextBox < RoundSpin.Min Then
    ErrorLabel.Caption = "Round num_digit should be between 0 to 15"
    RoundTextBox = 2
    Exit Sub
  End If
  On Error GoTo ERROR_HANDLE
  
  Set inputRange = Range(InputEdit.Value)
  Set outputRangeStart = Range(OutputEdit.Value)
    
  StartTimer
  If InverseFFTCheck.Value = True Then
    Call IFFT(inputRange, outputRangeStart, RoundTextBox)
    SuccessLabel.Caption = "IFFT work is done (" & EndTimer & "ms)"
  Else
    Call FFT(inputRange, outputRangeStart, RoundTextBox)
    SuccessLabel.Caption = "FFT work is done (" & EndTimer & "ms)"
  End If
  
  Exit Sub
  
ERROR_HANDLE:
  ErrorLabel.Caption = "Choose right data"
End Sub

Private Sub TextBox1_Change()

End Sub

Private Sub UserForm_Initialize()
  Dim inputRange As Range, outputRange As Range
  
  FFT_Form.Height = 350
  FFT_Form.Width = 450

  '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

Design of FFT_Form

 


Excel file containing whole program code: 

fft_v2.0 - distribution - encrypted.xlsm
2.91MB

This file is protected by password.
If someone want to use this file, please request password by writing at this page's comment and I'll reply you. Otherwise send e-mail to p14jeffwest@gmail.com. Then I'll let you know the password freely.

The password is: 2.718282

 

 

 Before: 07-4. 엑셀에서 직접 FFT 프로그래밍 작성하기
 Next: 07-5. 가장 빠른 Java용 FFT 구현해보기

 

반응형