vb 傅里叶变换

VB FFT 快速傅里叶变换/反变换2008-10-19 22:11模块****************************************** **************
'FFT0 数组下标以0开始FFT1 数组下标以1开始
'AR() 数据实部AI() 数据虚部
'N 数据点数，为2的整数次幂
'NI 变换方向1为正变换，-1为反变换
'****************************************** *********************
Public Const Pi = 3.1415926
Public Function 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 Function
Public Function 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 Function
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
'*使用**********
Const fftIn = 128
Dim i As Integer
Dim xr(128) As Double
Dim xi(128) As Double
'赋值，IaIn(i)是采得的数据。

For i = 0 To 128
xr(i) = 100 * IaIn(i)
xi(i) = 0
Next
'FFT变换
Call FFT0(xr(), xi(), 128, 1)
'绘图
picI_FFT.Scale (0, 100)-(fftIn - 1, -10)
picI_FFT.DrawWidth = 2
For i = 0 To fftIn - 1
picI_FFT.Line (i, Abs(xr(i)))-(i + 1, Abs(xr(i + 1))), vbBlue
Next i
谁会FFT(快速傅里叶变换)??VB编的模块!谁看得懂源代码??假如以SIN（X)为例，画频谱图，在这个模块中的realin和imagin分别代表什么？？realout和imagout又是代表什么??怎么用这个模块了!!以画Y=sin(X)为例的话,x=1--1024.是不是realin=sin(X), 而imagin=0呢??输出的时候横坐标是什么呢??纵坐标是什么呢??我只画实频图,也就是只要变化后的实部吗??调用的时候是这样的吗??FourierTransform(16,realin(),imageout(),realout(),imagout()) 假设事
先定义了dim realin(1000),imagin(1000),realout(1000),imagout(1000).希望哪位高手能指点小弟我了！！小弟我只会一点点了！请留言的时候，详细点了。

谢谢了！！我的Email 地址是:quanquan3031@. 下面是FFT（快速傅里叶变换）的VB源代码：'--------------------------------------------------------------------
' VB FFT Release 1.0
' by Murphy McCauley (MurphyMc@)
' 08/01/99
'--------------------------------------------------------------------
' This code is very, very heavily based on Don Cros s 's fourier.pas
' Turbo Pascal Unit for calculating the Fast Fourier Tr ansform.
' I 've not implemented all of his functions, though I may well do
' so in the future.
' For more info, you can contact me by email, chec k my website at:
' /deeth/
' or check Don Cross 's FFT web page at:
' /~dcross/fft.html
' You also may be intrested in the FFT.DLL that I put together based
' on Don Cross 's FFT C code. It 's callable with Visual Basic and
' includes VB declares. You can get it from his we bsite.
'--------------------------------------------------------------------
Option Explicit
Const Pi = 3.14159265358979
Function NumberOfBitsNeeded(PowerOfTwo As Long) As Byt e
Dim I As Byte
For I = 0 To 16
If (PowerOfTwo And (2 ^ I)) <> 0 Then
NumberOfBitsNeeded = I
Exit Function
End If
Next
End Function
Function IsPowerOfTwo(X As Long) As Boolean
If (X < 2) Then IsPowerOfTwo = False: Exit F unction
If (X And (X - 1)) = False Then IsPowerOfTw o = True
End Function
Function ReverseBits(ByV al Index As Long, NumBits As Byte) As Long
Dim I As Byte, Rev As Long
For I = 0 To NumBits - 1
Rev = (Rev * 2) Or (Index And 1)
Index = Index \ 2
Next
ReverseBits = Rev
End Function
Sub FourierTransform(NumSamples As Long, realin() As Do uble, ImageIn() As Double, RealOut() As Double, ImagOu t() As Double)
Dim AngleNumerator As Double: AngleNumerator = 2# * Pi
Dim NumBits As Byte, I As Long, j As Long, k As Long, n As Long, BlockSize As Long, Block End As Long
Dim DeltaAngle As Double, DeltaAr As Double
Dim Alpha As Double, Beta As Double
Dim TR As Double, TI As Double, AR As Do uble, AI As Double
If (IsPowerOfTwo(NumSamples) = False) Or (NumSam ples < 2) Then
Call MsgBox( "Error in procedure Fourier: " + vbCrLf + " NumSamples is " + CStr(NumSamples ) + ", which is not a positive integer power of tw o. ", , "Error! ")
Exit Sub
NumBits = NumberOfBitsNeeded(NumSamples)
For I = 0 To (NumSamples - 1)
j = ReverseBits(I, NumBits)
RealOut(j) = realin(I)
ImagOut(j) = ImageIn(I)
Next
BlockEnd = 1
BlockSize = 2
Do While BlockSize <= NumSamples
DeltaAngle = AngleNumerator / BlockSize
Alpha = Sin(0.5 * DeltaAngle)
Alpha = 2# * Alpha * Alpha
Beta = Sin(DeltaAngle)
I = 0
Do While I < NumSamples
AR = 1#
AI = 0#
j = I
For n = 0 To BlockEnd - 1
k = j + BlockEnd
TR = AR * RealOut(k) - AI * ImagOut(k)
TI = AI * RealOut(k) + AR * ImagOut(k)
RealOut(k) = RealOut(j) - TR
ImagOut(k) = ImagOut(j) - TI
RealOut(j) = RealOut(j) + TR
ImagOut(j) = ImagOut(j) + TI
DeltaAr = Alpha * AR + Beta * AI
AI = AI - (Alpha * AI - B eta * AR)
AR = AR - DeltaAr
j = j + 1
Next
I = I + BlockSize
Loop
BlockEnd = BlockSize
BlockSize = BlockSize * 2 Loop
End Sub。