FFT，快速傅里也变换

fft程序， Dim REX(512) 'REX[ ] holds the real part of the frequency domain Dim IMX(512) 'IMX[ ] holds the imaginary part of the frequency domain Const N = 512 Public Sub fft() pi = 3.14159265 'Set constants 1000 'THE FAST FOURIER TRANSFORM 'copyright © 1997-1999 by California Technical Publishing 'published with permission from Steven W Smith, www.dspguide.com 'GUI by logix4u , www.logix4u.net 'modified by logix4u, www.logix4.net 1010 'Upon entry, N% contains the number of points in the DFT, REX[ ] and 1020 'IMX[ ] contain the real and imaginary parts of the input. Upon return, 1030 'REX[ ] and IMX[ ] contain the DFT output. All signals run from 0 to N%-1. 1060 NM1% = N% - 1 1070 ND2% = N% / 2 1080 M% = CInt(Log(N%) / Log(2)) 1090 J% = ND2% 1100 ' 1110 For i% = 1 To N% - 2 'Bit reversal sorting 1120 If i% >= J% Then GoTo 1190 1130 TR = REX(J%) 1140 TI = IMX(J%) 1150 REX(J%) = REX(i%) 1160 IMX(J%) = IMX(i%) 1170 REX(i%) = TR 1180 IMX(i%) = TI 1190 K% = ND2% 1200 If K% > J% Then GoTo 1240 1210 J% = J% - K% 1220 K% = K% / 2 1230 GoTo 1200 1240 J% = J% + K% 1250 Next i% 1260 ' 1270 For L% = 1 To M% 'Loop for each stage 1280 LE% = CInt(2 ^ L%) 1290 LE2% = LE% / 2 1300 UR = 1 1310 UI = 0 1320 SR = Cos(pi / LE2%) 'Calculate sine & cosine values 1330 SI = -Sin(pi / LE2%) 1340 For J% = 1 To LE2% 'Loop for each sub DFT 1350 JM1% = J% - 1 1360 For i% = JM1% To NM1% Step LE% 'Loop for each butterfly 1370 IP% = i% + LE2% 1380 TR = REX(IP%) * UR - IMX(IP%) * UI 'Butterfly calculation 1390 TI = REX(IP%) * UI + IMX(IP%) * UR 1400 REX(IP%) = REX(i%) - TR 1410 IMX(IP%) = IMX(i%) - TI 1420 REX(i%) = REX(i%) + TR 1430 IMX(i%) = IMX(i%) + TI 1440 Next i% 1450 TR = UR 1460 UR = TR * SR - UI * SI 1470 UI = TR * SI + UI * SR 1480 Next J% 1490 Next L% 1500 ' End Sub Private Sub btn_clr_Click() ipplot.Cls FuncGen.Clear End Sub Private Sub btn_generate_Click() FuncGen.SamplingRate = CInt(txtsamp.Text) FuncGen.amplitude = Slideramp.Value If opt_sin.Value = True Then FuncGen.GenSine CInt(txtfreq.Text) End If If opt_cos.Value = True Then FuncGen.GenCos CInt(txtfreq.Text) End If If opt_impulse.Value = True Then FuncGen.GenImpulse Val(txt_impulsepos.Text) End If If Opt_sqr.Value = True Then FuncGen.GenSquare CInt(txtfreq.Text) End If ipplot.Cls plotip End Sub Private Sub Command1_Click() output.Cls End Sub Private Sub Command2_Click() On Error Resume Next If Optreal.Value = True Then output.ForeColor = vbRed For cnt = 0 To 156 outputarray(cnt) = REX(cnt) Next cnt End If If Optimg.Value = True Then output.ForeColor = vbWhite For cnt = 0 To 156 outputarray(cnt) = IMX(cnt) Next cnt End If If Optmag.Value = True Then output.ForeColor = vbYellow For cnt = 0 To 156 outputarray(cnt) = Sqr((IMX(cnt) * IMX(cnt)) + (REX(cnt) * REX(cnt))) Next cnt End If If Optphase.Value = True Then output.ForeColor = vbCyan For cnt = 0 To 156 outputarray(cnt) = Atn(IMX(cnt) / REX(cnt)) Next cnt End If If Optpower.Value = True Then output.ForeColor = vbGreen For cnt = 0 To 156 outputarray(cnt) = (IMX(cnt) * IMX(cnt)) + (REX(cnt) * REX(cnt)) Next cnt End If For cnt = 0 To 256 If (maxvalue < outputarray(cnt)) Then maxvalue = outputarray(cnt) End If Next cnt oldval = (2600 / (2 * Screen.TwipsPerPixelY)) - (outputarray(1) * ((950 / Screen.TwipsPerPixelY) / maxvalue)) For cnt = 1 To 256 output.Line (cnt - 1, oldval)-(cnt, (2600 / (2 * Screen.TwipsPerPixelY)) - (outputarray(cnt) * ((950 / Screen.TwipsPerPixelY) / maxvalue))) oldval = (2600 / (2 * Screen.TwipsPerPixelY)) - (outputarray(cnt) * ((950 / Screen.TwipsPerPixelY) / maxvalue)) Next cnt End Sub Private Sub Command3_Click() MsgBox "Developed and published by : " + vbNewLine + "LOGIX4U" + vbNewLine + "www.logix4u.net" + vbNewLine + vbNewLine + "Algorithms by : " + vbNewLine + "Steven W Smith " + vbNewLine + "www.dspguide.com", vbInformation, "About" End Sub Private Sub Command4_Click() On Error Resume Next Label5.Visible = True DoEvents If Opthamm.Value = True Then FuncGen.ApplyHamming ipplot.Cls plotip End If If Opthann.Value = True Then FuncGen.ApplyHanning ipplot.Cls plotip End If If Optblack.Value = True Then FuncGen.ApplyBlackman ipplot.Cls plotip End If For cnt = 1 To 512 REX(cnt) = FuncGen.Samples(cnt) IMX(cnt) = 0 Next cnt fft If Optreal.Value = True Then output.ForeColor = vbRed For cnt = 0 To 156 outputarray(cnt) = REX(cnt) Next cnt End If If Optimg.Value = True Then output.ForeColor = vbWhite For cnt = 0 To 156 outputarray(cnt) = IMX(cnt) Next cnt End If If Optmag.Value = True Then output.ForeColor = vbYellow For cnt = 0 To 156 outputarray(cnt) = Sqr((IMX(cnt) * IMX(cnt)) + (REX(cnt) * REX(cnt))) Next cnt End If If Optphase.Value = True Then output.ForeColor = vbCyan For cnt = 0 To 156 outputarray(cnt) = Atn(IMX(cnt) / REX(cnt)) Next cnt End If If Optpower.Value = True Then output.ForeColor = vbGreen For cnt = 0 To 156 outputarray(cnt) = Sqr((IMX(cnt) * IMX(cnt)) + (REX(cnt) * REX(cnt))) Next cnt End If For cnt = 0 To 256 If (maxvalue < outputarray(cnt)) Then maxvalue = outputarray(cnt) End If Next cnt oldval = (2600 / (2 * Screen.TwipsPerPixelY)) - (outputarray(1) * ((950 / Screen.TwipsPerPixelY) / maxvalue)) For cnt = 1 To 256 output.Line (cnt - 1, oldval)-(cnt, (2600 / (2 * Screen.TwipsPerPixelY)) - (outputarray(cnt) * ((950 / Screen.TwipsPerPixelY) / maxvalue))) oldval = (2600 / (2 * Screen.TwipsPerPixelY)) - (outputarray(cnt) * ((950 / Screen.TwipsPerPixelY) / maxvalue)) Next cnt Label5.Visible = False End Sub Private Sub Form_Activate() ipplot.Line (0, 1450 / Screen.TwipsPerPixelY)-(10000, 1450 / Screen.TwipsPerPixelY), vbRed End Sub Private Sub Form_Paint() On Error Resume Next Dim maxvalue, oldval As Double ipplot.Cls plotip End Sub Private Sub opt_impulse_Click() If opt_impulse.Value = True Then txt_impulsepos.Enabled = True End If End Sub 快速傅里叶变换（FFT）是一种高效的算法，用于计算离散傅里叶变换（DFT）和其逆变换。在给定的代码中，我们看到了一个VBA实现的FFT算法，用于处理复数序列的DFT。以下是这个FFT算法的工作原理和关键步骤的详细解释： 1. **初始化和变量声明**： - `REX` 和 `IMX` 数组分别存储输入信号的实部和虚部。 - `N` 定义了DFT的点数，即处理的样本数量。 - `fft` 子程序执行FFT运算。 2. **位反序排序**（Bit Reversal）： - 这个过程将输入序列重新排列，以便于后续的蝶形运算。 - `i%` 和 `J%` 用于遍历数组，`K%` 控制当前的位反序级别。 - 如果 `i% >= J%`，则交换 `REX[i%]` 和 `REX[J%]` 的值以及它们对应的虚部。 3. **蝶形运算**（Butterfly Operations）： - 这是FFT的核心部分，通过一系列乘以复数因子（sine 和 cosine 值）的加减运算，逐步构建出频域表示。 - `L%` 表示阶段数，与二进制位数相关。 - `LE%` 和 `LE2%` 分别是每个阶段处理的子序列长度及其一半。 - `UR` 和 `UI` 是初始的cosine和-sine值。 - `TR` 和 `TI` 用于临时存储中间结果。 - 在每个子DFT中，使用`J%` 和 `i%` 遍历并执行蝶形运算，更新 `REX` 和 `IMX`。 4. **其他辅助功能**： - `btn_generate_Click` 事件处理函数生成不同类型的信号（如正弦、余弦、脉冲或方波），并设置采样率和幅度。 - `Command1_Click` 清除输出。 - `Command2_Click` 显示DFT的结果，包括实部、虚部、模值、相位和功率谱。 - `Command3_Click` 显示关于FFT的信息。 - `opt_impulse_Click` 控制是否启用脉冲信号生成。 5. **显示和输出**： - `ipplot` 和 `output` 两个组件可能用于图形化显示输入信号和DFT结果。 - 通过不同的颜色和线型，用户可以选择查看实部、虚部、模值、相位或功率谱。 这段代码提供了一个完整的FFT实现，包括输入信号的生成、位反序排序、蝶形运算和结果展示。这个实现适用于教育和简单分析，但可能不包含更复杂的情况，如窗函数应用或其他优化。在实际的工程应用中，通常会使用更为优化的库函数，如MATLAB的`fft`或Python的`numpy.fft`。