modfft.bas

我自己编写的VB的FFT程序

Attribute VB_Name = "modFFT"
Option Explicit
Public Const PI = 3.1415926
Type COMPLEX
    re As Double
    im As Double
End Type
Public TD() As COMPLEX
Public FD() As COMPLEX
Public W(511) As COMPLEX
Public X1(1023) As COMPLEX
Public Ts As Double
Public Fs As Double
'复数加运算
Public Function COMPLEX_ADD(c1 As COMPLEX, c2 As COMPLEX) As COMPLEX
    Dim c As COMPLEX
    c.re = c1.re + c2.re
    c.im = c1.im + c2.im
    COMPLEX_ADD.re = c.re
    COMPLEX_ADD.im = c.im
End Function
'复数减运算
Public Function COMPLEX_DEC(c1 As COMPLEX, c2 As COMPLEX) As COMPLEX
    Dim c As COMPLEX
    c.re = c1.re - c2.re
    c.im = c1.im - c2.im
    COMPLEX_DEC.re = c.re
    COMPLEX_DEC.im = c.im
End Function
'复数乘运算
Public Function COMPLEX_MUL(c1 As COMPLEX, c2 As COMPLEX) As COMPLEX
    Dim c As COMPLEX
    c.re = c1.re * c2.re - c1.im * c2.im
    c.im = c1.re * c2.im + c1.im * c2.re
    COMPLEX_MUL.re = c.re
    COMPLEX_MUL.im = c.im
End Function

Public Sub FFT(ByVal power As Integer, ByVal mode As Boolean)
    Dim count As Integer
    Dim angle As Double
    Dim i As Integer
    Dim j As Integer
    Dim k As Integer
    Dim p As Integer
    Dim bfsize As Integer
    Dim wCount As Integer
    count = 2 ^ power
    wCount = count \ 2
 '   ReDim FD(count - 1) As COMPLEX
  '  Dim W(wCount) As COMPLEX
   ' Dim X1(count - 1) As COMPLEX
    'Dim X2(count - 1) As COMPLEX
    Dim X As COMPLEX
    Dim Xa As COMPLEX
    Dim Xb As COMPLEX
    
    '计算加权系数
    For i = 0 To count \ 2 - 1
        angle = i * PI * 2 / count
        W(i).re = Cos(angle)
        W(i).im = Sin(angle)
    Next
    If mode = True Then
    For i = 0 To count - 1
        X1(i).im = TD(i).im
        X1(i).re = TD(i).re
    Next
    End If
    '蝶形运算
    For k = 0 To power - 1      '迭代次数
        For j = 0 To 2 ^ k - 1  '蝶形单元数
            bfsize = 2 ^ (power - k)
            For i = 0 To bfsize \ 2 - 1     '每个蝶形单元循环次数
                p = j * bfsize
                'X = COMPLEX_MUL(X1(i + p + bfsize \ 2), W(i * (2 ^ k)))
                'Xa = COMPLEX_DEC(X1(i + p), X)
                'Xb = COMPLEX_ADD(X1(i + p), X)
                'X1(i + p + bfsize \ 2) = Xa
                'X1(i + p) = Xb
                Xa = COMPLEX_ADD(X1(i + p), X1(i + p + bfsize \ 2))
                Xb = COMPLEX_MUL(COMPLEX_DEC(X1(i + p), X1(i + p + bfsize \ 2)), W(i * (2 ^ k)))
                X1(i + p) = Xa
                X1(i + p + bfsize \ 2) = Xb
                'X2(i + p) = COMPLEX_ADD(X1(i + p), COMPLEX_MUL(X1(i + p + bfsize / 2), W(i * (2 ^ k))))
                'X2(i + p + bfsize / 2) = COMPLEX_DEC(X1(i + p), COMPLEX_MUL(X1(i + p + bfsize / 2), W(i * (2 ^ k))))
            Next
        Next
    Next
    '重新排序
    If mode = True Then
    For j = 0 To count - 1
        p = 0
        k = j
        For i = power - 1 To 0 Step -1
            '位测试运算
            If k >= 2 ^ i Then
                p = p + 2 ^ (power - i - 1)
                k = k - 2 ^ i
            End If
        Next
        FD(j).im = X1(p).im
        FD(j).re = X1(p).re
    Next
    Else
    For j = 0 To count - 1
        p = 0
        k = j
        For i = power - 1 To 0 Step -1
            '位测试运算
            If k >= 2 ^ i Then
                p = p + 2 ^ (power - i - 1)
                k = k - 2 ^ i
            End If
        Next
        If j < p Then
        X.im = X1(p).im
        X.re = X1(p).re
        X1(p).re = X1(j).re
        X1(p).im = X1(j).im
        X1(j).re = X.re
        X1(j).im = X.im
        End If
    Next
    End If
End Sub
Public Sub IFFT(ByVal power As Integer)
    Dim count As Integer
    Dim i As Integer
    '计算傅里叶反变换点数
    count = 2 ^ power
    '求频域点的共扼
    For i = 0 To count - 1
        X1(i).re = FD(i).re
        X1(i).im = -FD(i).im
    Next
    '调用快速傅里叶变换
    Call FFT(power, False)
    '求时域点的共扼
    For i = 0 To count - 1
        X1(i).re = X1(i).re / count
        X1(i).im = -X1(i).im / count
    Next
End Sub
Public Sub Test()
    ReDim TD(1023) As COMPLEX
    ReDim FD(1023) As COMPLEX
    Dim i As Integer
    For i = 0 To 1023
        TD(i).re = 5 * Sin(2 * PI * i / 128) + 4 * Sin(2 * PI * i / 64) + 2 * Sin(2 * PI * i / 32) + 1 * Sin(2 * PI * i / 16)
        TD(i).im = 0
    Next
    Ts = 1 / 8000
    Fs = 8000 / 1024
    Call FFT(10, True)
    Call IFFT(10)
End Sub