|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
本帖最后由 VibInfo 于 2016-11-9 15:09 编辑
FFT的VB实现- -Tag: FFT VB
'***************************************************************
'FFT0 数组下标以0开始 FFT1 数组下标以1开始
'AR() 数据实部 AI() 数据虚部
'N 数据点数,为2的整数次幂
'NI 变换方向 1为正变换,-1为反变换
'***************************************************************
Public Const Pi = 3.1415926
Public Sub FFT0(AR() As Double, AI() As Double, N As Integer, NI As Integer)
Dim I As Integer, J As Integer, K As Integer, L As Integer, M As Integer
Dim IP As Integer, LE As Integer
Dim L1 As Integer, N1 As Integer, N2 As Integer
Dim SN As Double, TR As Double, TI As Double, WR As Double, WI As Double
Dim UR As Double, UI As Double, US As Double
M = NTOM(N)
N2 = N / 2
N1 = N - 1
SN = NI
J = 1
For I = 1 To N1
If I < J Then
TR = AR(J - 1)
AR(J - 1) = AR(I - 1)
AR(I - 1) = TR
TI = AI(J - 1)
AI(J - 1) = AI(I - 1)
AI(I - 1) = TI
End If
K = N2
While (K < J)
J = J - K
K = K / 2
Wend
J = J + K
Next I
For L = 1 To M
LE = 2 ^ L
L1 = LE / 2
UR = 1#
UI = 0#
WR = Cos(Pi / L1)
WI = SN * Sin(Pi / L1)
For J = 1 To L1
For I = J To N Step LE
IP = I + L1
TR = AR(IP - 1) * UR - AI(IP - 1) * UI
TI = AI(IP - 1) * UR + AR(IP - 1) * UI
AR(IP - 1) = AR(I - 1) - TR
AI(IP - 1) = AI(I - 1) - TI
AR(I - 1) = AR(I - 1) + TR
AI(I - 1) = AI(I - 1) + TI
Next I
US = UR
UR = US * WR - UI * WI
UI = UI * WR + US * WI
Next J
Next L
If SN <> -1 Then
For I = 1 To N
AR(I - 1) = AR(I - 1) / N
AI(I - 1) = AI(I - 1) / N
Next I
End If
End Sub
Public Sub FFT1(AR() As Double, AI() As Double, N As Integer, NI As Integer)
Dim I As Integer, J As Integer, K As Integer, L As Integer, M As Integer
Dim IP As Integer, LE As Integer
Dim L1 As Integer, N1 As Integer, N2 As Integer
Dim SN As Double, TR As Double, TI As Double, WR As Double, WI As Double
Dim UR As Double, UI As Double, US As Double
M = NTOM(N)
N2 = N / 2
N1 = N - 1
SN = NI
J = 1
For I = 1 To N1
If I < J Then
TR = AR(J)
AR(J) = AR(I)
AR(I) = TR
TI = AI(J)
AI(J) = AI(I)
AI(I) = TI
End If
K = N2
While (K < J)
J = J - K
K = K / 2
Wend
J = J + K
Next I
For L = 1 To M
LE = 2 ^ L
L1 = LE / 2
UR = 1#
UI = 0#
WR = Cos(Pi / L1)
WI = SN * Sin(Pi / L1)
For J = 1 To L1
For I = J To N Step LE
IP = I + L1
TR = AR(IP) * UR - AI(IP) * UI
TI = AI(IP) * UR + AR(IP) * UI
AR(IP) = AR(I) - TR
AI(IP) = AI(I) - TI
AR(I) = AR(I) + TR
AI(I) = AI(I) + TI
Next I
US = UR
UR = US * WR - UI * WI
UI = UI * WR + US * WI
Next J
Next L
If SN <> -1 Then
For I = 1 To N
AR(I) = AR(I) / N
AI(I) = AI(I) / N
Next I
End If
End Sub
Private Function NTOM(N As Integer) As Integer
Dim ND As Double
ND = N
NTOM = 0
While (ND > 1)
ND = ND / 2
NTOM = NTOM + 1
Wend
End Function
Public Sub FFT(INr#(), INi#(), n%, Mm%, TT#, NPP%)
Dim FFTn1%, FFTn2%, FFTsn%, FFTj%, FFTi%, FFTip%, tmpR#, tmpI#, FFTk%, FFTL%, FFTLE#, FFTL1#, FFTus#, FFTur#, FFTui#, FFTwr#, FFTwi#
FFTn2 = n / 2
FFTn1 = n - 1
FFTsn = NPP
FFTj = 1
For FFTi = 1 To FFTn1
If (FFTi >= FFTj) Then GoTo X25
tmpR = INr(FFTj)
INr(FFTj) = INr(FFTi)
INr(FFTi) = tmpR
tmpI = INi(FFTj)
INi(FFTj) = INi(FFTi)
INi(FFTi) = tmpI
X25: FFTk = FFTn2
X45: If (FFTk >= FFTj) Then GoTo X35
FFTj = FFTj - FFTk
FFTk = FFTk / 2
GoTo X45
X35: FFTj = FFTj + FFTk
Next FFTi
For FFTL = 1 To Mm
FFTLE = 2 ^ FFTL
FFTL1 = FFTLE / 2
FFTur = 1
FFTui = 0
FFTwr = Cos(3.1415926 / FFTL1)
FFTwi = FFTsn * Sin(3.1415926 / FFTL1)
For FFTj = 1 To FFTL1
For FFTi = FFTj To n Step FFTLE
FFTip = FFTi + FFTL1
tmpR = INr(FFTip) * FFTur - INi(FFTip) * FFTui
tmpI = INi(FFTip) * FFTur + INr(FFTip) * FFTui
INr(FFTip) = INr(FFTi) - tmpR
INi(FFTip) = INi(FFTi) - tmpI
INr(FFTi) = INr(FFTi) + tmpR
INi(FFTi) = INi(FFTi) + tmpI
Next FFTi
FFTus = FFTur
FFTur = FFTus * FFTwr - FFTui * FFTwi
FFTui = FFTui * FFTwr + FFTus * FFTwi
Next FFTj
Next FFTL
If (FFTsn = -1) Then
For FFTi = 1 To n
INr(FFTi) = INr(FFTi) * TT
INi(FFTi) = INi(FFTi) * TT
Next FFTi
Else
For FFTi = 1 To n
INr(FFTi) = INr(FFTi) / n / TT
INi(FFTi) = INi(FFTi) / n / TT
Next FFTi
End If
End Sub
[ 本帖最后由 zhlong 于 2007-6-4 21:50 编辑 ]
|
|