fft.bas

自编FFT谱分析中的栅栏效应演示

Attribute VB_Name = "Module1"
Public Declare Sub CopyMemory Lib "kernel32" Alias "RtlMoveMemory" (Destination As Any, Source As Any, ByVal Length As Long)

Public Type Complex
  Re As Double
  Im As Double
End Type

Public Const PI = 3.1415926535

Public TD(1023) As Complex, FD(1023) As Complex
Public y(1023) As Double

Public Function Complex_Compue(a As Complex, b As Complex, Mark As Byte) As Complex
  '复数运算
  Select Case Mark
    Case 1 '复数加
      Complex_Compue.Re = a.Re + b.Re
      Complex_Compue.Im = a.Im + b.Im
    Case 2 '复数减
      Complex_Compue.Re = a.Re - b.Re
      Complex_Compue.Im = a.Im - b.Im
    Case 3 '复数乘
      Complex_Compue.Re = a.Re * b.Re - a.Im * b.Im
      Complex_Compue.Im = a.Re * b.Im + a.Im * b.Re
  End Select
  
End Function

Public Sub FFT(TD() As Complex, FD() As Complex, m As Long)

  Dim i As Integer, j As Integer, k As Integer, bfSize As Integer, p As Integer
  
  Dim Angle As Double
    
  Dim N As Long
  
  N = 2 ^ m
  
  ReDim X1(N - 1) As Complex, X2(N - 1) As Complex, X(N - 1) As Complex, W(N \ 2 - 1) As Complex ', TD(n - 1) As Complex, FD(n - 1) As Complex

  For i = 0 To N \ 2 - 1
  
    Angle = -i * PI * 2 / N
    
    With W(i)
      .Re = Cos(Angle)
      .Im = Sin(Angle)
    End With
    
  Next i
  
  CopyMemory X1(0), TD(0), 16 * N
  
  For k = 0 To m - 1 '进行蝶形图运算
  
    For j = 0 To 2 ^ k - 1
    
      bfSize = 2 ^ (m - k)
      
      For i = 0 To bfSize \ 2 - 1
      
        p = j * bfSize
        
        X2(i + p) = Complex_Compue(X1(i + p), X1(i + p + bfSize \ 2), 1)
           
        X2(i + p + bfSize \ 2) = Complex_Compue(Complex_Compue(X1(i + p), X1(i + p + bfSize \ 2), 2), W(i * (2 ^ k)), 3)
      
      Next i
      
    Next j
    
    '交换X1和X2的数据
    CopyMemory X(0), X1(0), 16 * N
    CopyMemory X1(0), X2(0), 16 * N
    CopyMemory X2(0), X(0), 16 * N
    
  Next k
  
  For j = 0 To N - 1 '倒序
  
    p = 0
    
    For i = 0 To m - 1
    
      If (j And (2 ^ i)) = (2 ^ i) Then p = p + 2 ^ (m - i - 1)
      
    Next i
    
    FD(j) = X1(p)
    
  Next j
  
End Sub