如何实现实数傅立叶变换的算法?

Question

现有一些采样数据,要进行傅立叶变换该如何实现

Answer 1

http://www.21ic.com/new_info/news/files/news/2003126140744.asp

Answer 2

Option ExplicitPublic Const AngleNumerator = 6.283185   ' 2 * Pi = 2 * 3.14159265358979
Public Const NumSamples = 1024 '1024
Public Const NumBits = 10Private ReversedBits(0 To NumSamples - 1) As LongSub DoReverse()
    Dim i As Long
    For i = LBound(ReversedBits) To UBound(ReversedBits)
        ReversedBits(i) = ReverseBits(i, NumBits)
    Next
End SubFunction 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
'        我们可以把FFT简单地看作一个变换器，输入N+1个数，输出N+1个数
'        ，但他们对应的意义不同，如果把输入当作时域，则输出为频域，表
'        怔了其对应域的变化快慢?
'
'        假设输入信号本身的频率为fc(或者说频带宽为fc)，被频率为fs的冲击
'        串采样（由采样定理，fs >= 2*fc），则变换前的N＋1个数字对应的x
'        轴为{t0,t1,…tN}={0,Ts,2*Ts,....,N*Ts} (其中Ts为1/fs,为采样周期)
'
'        则变换后的N+1个数对应的x轴变为频率，范围为0~fs，以fs/N为间隔，
'        既为频率点{0,fs/N,2*fs/N,……,fs}，在matlab中如果用fftshift(fft(data))
'        ，则变换后对应x轴为-fs/2~fs/2，如果满足采样定理的化，信号频带-fc~fc
'        就包含在转换后的频谱里面了，就不会有失真。
'
'        注意：变换后的数字为复数，因为其中包括了幅度的信
'        息，abs(fftshift(fft(data)))为幅度，angle(fftshift(fft(data)))为相位Sub FFTAudio(RealIn() As Integer, RealOut() As Single)    Static ImagOut(0 To NumSamples - 1) As Single
    
    Static i As Long, j As Long, k As Long, n As Long, BlockSize As Long, BlockEnd As Long
    Static DeltaAngle As Single, DeltaAr As Single
    Static Alpha As Single, Beta As Single
    Static TR As Single, TI As Single, AR As Single, AI As Single
    
    For i = 0 To (NumSamples - 1)
        j = ReversedBits(i)
        RealOut(j) = RealIn(i)
        ImagOut(j) = 0
    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
End Sub