Excel的VBA实现快速傅里叶变换

Question

我正在尝试在微软的Excel VBA中实现一个快速傅里叶变换（基数-2）。我使用的代码从工作表中的范围提取数据，进行计算，然后将结果转储到相邻的列中。我遇到的问题是：1）不知道如何处理得到的X [k]数组；2）将这些结果与Excel内置的FFT结果匹配（它们目前不匹配）。下面是代码。感谢您的帮助。

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

上面的子程序通过For / Next循环调用下面的子程序来计数“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

- GollumInATuxedo

---

由于内置的FFT可能是用C实现的，因此它比VBA中的任何东西快至少一个数量级。您不想使用它的原因是什么？可以从VBA调用它。请参见：http://www.cpearson.com/excel/atp.htm - John Coleman

是的，我知道相比于其他很多编程语言，VBA运行速度会非常慢。内置FFT的输入数组限制为4096，但我想要通过FFT运行更大的数组（尽管时间明显增加）。此外，我更愿意学习实际算法，而不是让程序替我完成工作。 - GollumInATuxedo

我理解你的意思 - 我只是认为你可能不知道可以直接在VBA中使用它。自己实现这样的算法肯定是一次学习经历。从更实际的角度来看 - 也许你可以看看这个：https://newtonexcelbach.wordpress.com/2015/05/05/xlscipy-1-01/（即使你不选择这条路，包含此帖子的博客对于任何想在Excel中进行硬核科学计算的人来说都是一个好的阅读材料） - John Coleman

谢谢，我一定会仔细阅读的。 - GollumInATuxedo

Answer 1

我重新回顾了整个流程，发现我的问题在于我将错的值赋给了旋转因子TFactor_N1和TFactor_N2。在修复这个问题并调整显示哪些值后，我能够获得与Excel内置FFT相同的结果。修正后的代码如下所示。

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

- GollumInATuxedo

Answer 2

在Excel VBA中实现FFT有点复杂，但并不太难。对于典型应用程序，输入信号通常是实值的，而不是复值的。如果您正在从速度或加速度传感器或麦克风测量动态信号，则会出现这种情况。 这里展示了一个DFT，它将转换任意数量的输入对（例如时间和速度）。它们不必是2 ^ N个数据（FFT所需）。通常，时间均匀分割，因此您只需要DeltaTime（数据之间的时间间隔）。让我在这里放入代码：

    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

- user19361336

Answer 3

函数调用不正确 应该是： Call DecimationInTime(WS, n, 2 ^ x, x, TFactor_N1, TFactor_N2)

- Bob

---

这并没有回答问题。一旦您拥有足够的声望，您将能够评论任何帖子；相反，提供不需要询问者澄清的答案。- 来自审核 - Ike

Answer 4

除了已发布的VBA解决方案之外，最近版本的Excel允许使用LAMBDA函数作为纯公式来实现FFT（即不需要任何VBA代码）。

其中一个这样的实现是https://github.com/altomani/XL-FFT。

对于2的幂长度，它使用递归基数-2 Cooley-Tukey算法，对于其他长度，使用Bluestein算法的一个版本将计算减少到2的幂情况。

- AndreA