vb-fft.txt

FFT实现序列离散傅立叶变换

'************************************************************************
'FFT0 数组下标以0开始  FFT1 数组下标以1开始
'AR() 数据实部         AI() 数据虚部
'N 数据点数，为2的整数次幂
'NI 变换方向 1为正变换，-1为反变换
'************************************************************************
Public Sub FFT(AR() As Double, AI() As Double, N As Integer)
Dim I As Integer, J As Integer, K As Integer, L As Integer, M As Integer
Dim P As Integer,IP as Integer, LE As Integer
Dim B as integer, N1 As Integer, N2 As Integer
dim TR As Double,TI As Double,TAR As Double, TAI As Double
Dim WR As Double, WI As Double
Dim Pi As Double
Pi = 3.141592654
M = NTOM(N)
N2 = N / 2
N1 = N - 2
'***************************************************************************
'序列倒序如：01234567－04261537;0123-0213
'***************************************************************************
J = N2

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
    If J >= K Then
       Do While (J >= K)
          J = J - K
          K = K / 2
       Loop
    End If
  J = J + K
Next I
'*******************************************************************************
'N点FFT
'XL(J)=XL-1(J)+XL-XL-1(J+B)WpN
'XL(J+B)=XL-1(J)-XL-1(J+B)WpN
'式中:p=J*2^(M-L)；J=0,1,…，2 L-1-1；L=1,2,…，M
'令IP＝J+B
'XL-1(J+B)WpN=XL-1(IP)WpN=TR+TI
'TR=AR(IP) * WR - AI(IP) * WI
'TI=AR(IP) * WI + AI(IP) * WR
'********************************************************************************
For L = 1 To M
    LE = 2 ^ L
    B = LE / 2

    For J = 0 To B-1
        P=J*2^ (M-L)
        WR=cos(Pi*2*p/N)
        WI=-sin(Pi*2*p/N)

        For K = J To N-1 Step LE
            IP=k+B
            TR = AR(IP) * WR - AI(IP) * WI
            TI = AR(IP) * WI + AI(IP) * WR
            TAR=AR(K)
            TAI=AI(K)
            AR(K)=TAR+TR
            AI(k)=TAI+TI
            AR(IP)=TAR-TR
            AI(IP)=TAI-TI
        Next K
   Next J
Next L
End Sub
'*******************************************************************************
'NTOM(N)返回log2(N)如N＝1024则NTOM(1024)=10
'*******************************************************************************
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