반응형
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:
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 구현해보기 |
반응형
'푸리에 변환, 신호 > 푸리에 변환의 모든 것' 카테고리의 다른 글
[07-5 Code] Fastest FFT code for Java (0) | 2022.08.05 |
---|---|
07-5. 가장 빠른 Java용 FFT 구현해보기 (2) | 2022.08.05 |
07-4. 엑셀에서 직접 FFT 프로그램 만들기 (7) | 2022.07.28 |
07-3. 엑셀에서 FFT 계산하기 (엑셀 자체 FFT 기능) (0) | 2022.07.28 |
07-2. FFT 예제를 손으로 풀어보며 이해하기 (5) | 2022.07.28 |