|
| 一段VB编的FFT，有疑问。急！！！ |
| 新闻出处：21ic
发布时间： 2007-08-14 |
rockzone 发布于 2007-8-13 20:39:00
程序如下：
Option Explicit
Const Pi = 3.14159265358979
Function NumberOfBitsNeeded(PowerOfTwo As Long) As Byte
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 Function
If (X And (X - 1)) = False Then IsPowerOfTwo = True
End Function
Function ReverseBits(ByVal 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 Double, ImageIn() As Double, RealOut() As Double, ImagOut() As Double, Optional InverseTransform As Boolean = False)
Dim AngleNumerator As Double
Dim NumBits As Byte, I As Long, j As Long, K As Long, n As Long, BlockSize As Long, BlockEnd 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 Double, AI As Double
If InverseTransform Then
AngleNumerator = -2# * Pi
Else
AngleNumerator = 2# * Pi
End If
If (IsPowerOfTwo(NumSamples) = False) Or (NumSamples < 2) Then
Call MsgBox("Error in procedure Fourier:" + vbCrLf + " NumSamples is " + CStr(NumSamples) + ", which is not a positive integer power of two.", , "Error!")
Exit Sub
End If
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 - Beta * AR)
AR = AR - DeltaAr
j = j + 1
Next
I = I + BlockSize
Loop
BlockEnd = BlockSize
BlockSize = BlockSize * 2
Loop
If InverseTransform Then
'Normalize the resulting time samples...
For I = 0 To NumSamples - 1
RealOut(I) = RealOut(I) / NumSamples
ImagOut(I) = ImagOut(I) / NumSamples
Next
End If
End Sub
Function FrequencyOfIndex(NumberOfSamples As Long, ByVal Index As Long) As Double
'Based on IndexToFrequency(). This name makes more sense to me.
If Index >= NumberOfSamples Then
FrequencyOfIndex = 0#
Exit Function
ElseIf Index <= NumberOfSamples / 2 Then
FrequencyOfIndex = CDbl(Index) / CDbl(NumberOfSamples)
Exit Function
Else
FrequencyOfIndex = -CDbl(NumberOfSamples - Index) / CDbl(NumberOfSamples)
Exit Function
End If
End Function
Sub CalcFrequency(NumberOfSamples As Long, FrequencyIndex As Long, RealIn() As Double, ImagIn() As Double, RealOut As Double, ImagOut As Double)
Dim K As Long
Dim Cos1 As Double, Cos2 As Double, Cos3 As Double, Theta As Double, Beta As Double
Dim Sin1 As Double, Sin2 As Double, Sin3 As Double
Theta = 2 * Pi * FrequencyIndex / CDbl(NumberOfSamples)
Sin1 = Sin(-2 * Theta)
Sin2 = Sin(-Theta)
Cos1 = Cos(-2 * Theta)
Cos2 = Cos(-Theta)
Beta = 2 * Cos2
For K = 0 To NumberOfSamples - 2
'Update trig values
Sin3 = Beta * Sin2 - Sin1
Sin1 = Sin2
Sin2 = Sin3
Cos3 = Beta * Cos2 - Cos1
Cos1 = Cos2
Cos2 = Cos3
RealOut = RealOut + RealIn(K) * Cos3 - ImagIn(K) * Sin3
ImagOut = ImagOut + ImagIn(K) * Cos3 + RealIn(K) * Sin3
Next
End Sub
以上程序各函数功能是什么，请高手讲解一下。
computer00 发布于 2007-8-13 21:16:00
NumberOfBitsNeeded 计算需要多少个bit，因为要做倒序排列，所以要知道有多少个bit，才能高低位掉转。
IsPowerOfTwo 判断它是不是2的整数次方。
ReverseBits 位反转用的，前面那个NumberOfBitsNeeded计算出需要多少位，然后再用这个函数将位反转。
FourierTransform 这个不用说了吧？就是FFT了。
FrequencyOfIndex 从名字上来看，是计算FFT结果中某点的频率值。
最后一个子过程， CalcFrequency，计算频率，不知道干啥用的...
rockzone 发布于 2007-8-13 21:32:00
RealIn()，Imagin（）是输入信号的实部和虚部数组
我做了一个数组正弦信号1024个数，但都是浮点的十进制数，怎么变成
输入信号的实部和须部呢？
rockzone 发布于 2007-8-13 21:38:00
请高手举一个例子，多谢了！！！！！
将军令 发布于 2007-8-14 9:25:00
而应该是慢速傅立叶变换】】
为什么不用C
在21ic，用basic作型号处理，是稀有物种哦！！！！！！！
rockzone 发布于 2007-8-14 14:26:00
C语言如何嵌套到VB中呢？
rockzone 发布于 2007-8-14 16:32:00
楼上的BOSS，能不能帮我说说这段程序我应该怎么用啊，
比如说我有1024个十进制的数，是个正弦波，
我想画出频域谱线，应该怎么引用这段程序呢？
computer00 发布于 2007-8-14 20:21:00
关于变换后的结果是什么意思，你可以去DSP版块搜索我以前发过的帖子。
rockzone 发布于 2007-8-14 20:45:00
结果计算的结果是实部＋虚部，
这样的数怎么画成谱线呢，
BOSS,你在DSP的帖子好多，我该看哪个呢，请明示！
|
| 【关闭】 【打印】 |
|
|
|
|
