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

uttx8gqw  于 2023-01-10  发布在  其他
关注(0)|答案(4)|浏览(692)

我在尝试实现快速傅立叶变换(Radix-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
u5rb5r59

u5rb5r591#

我回顾了整个过程,确定问题是我给旋转因子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
9bfwbjaz

9bfwbjaz2#

函数调用不好调用调用_核心. DecimationInTime(WS,n,2 ^x,x,TFactor_N1,TFactor_N2)
它应该是:调用及时抽取(WS,n,2 ^x,x,T系数_N1,T系数_N2)

iyfjxgzm

iyfjxgzm3#

在ExcelVBA中实现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
gab6jxml

gab6jxml4#

作为已经发布的VBA解决方案的替代方案,最新版本的Excel允许将FFT实现为具有LAMBDA函数的纯公式(即,没有任何VBA代码)。
https://github.com/altomani/XL-FFT就是这样一种实现。
对于2的幂长度,它使用递归基数-2 Cooley-Tukey algorithm,对于其他长度,使用Bluestein's algorithm的版本,该版本将计算减少到2的幂的情况。

相关问题