Дискретное преобразование Фурье на VB.NET

''' <summary>
''' Комплексное число.
''' </summary>
Public Class ComplexNumber

    ''' <summary>
    ''' Действительная часть комплексного числа.
    ''' </summary>
    Public Real As Double = 0

    ''' <summary>
    ''' Мнимая часть комплексного числа.
    ''' </summary>
    Public Imaginary As Double = 0

    Public Sub New()
        Real = 0
        Imaginary = 0
    End Sub

    ''' <summary>
    ''' Создаёт комплексное число.
    ''' </summary>
    ''' <param name="r">Действительная часть комплексного числа.</param>
    ''' <param name="im">Мнимая часть комплексного числа.</param>
    Public Sub New(ByVal r As Double, Optional ByVal im As Double = 0)
        Real = r
        Imaginary = im
    End Sub

    Private usCult As New Globalization.CultureInfo("en-US") 'используем культуру "en-US" чтобы целая и дробная части разделялись точкой, а не запятой

    ''' <summary>
    ''' Возвращает строку, состоящую из действительной и мнимой части, разделённых символом табуляции.
    ''' </summary>
    Public Overrides Function ToString() As String
        Return (Real.ToString(usCult) & ControlChars.Tab & Imaginary.ToString(usCult))
    End Function

End Class

''' <summary>
''' Рассчитывает спектр сигнала методом быстрого преобразования Фурье. Использовать только (N/2+1) возвращаемых значений (до половины частоты дискретизации).
''' </summary>
''' <param name="signal">Сигнал, содержащий количество отсчётов, кратное степени двойки, и состоящий из действительной и мнимой частей. Все мнимые части сигнала заполнены нулями.</param>
''' <returns>Возвращает массив комплексных чисел спектра. 
''' Значимы только первые N/2+1, остальные - симметричная часть, соответствующая отрицательным частотам.
''' Первое значение спектра - это постоянная составляющая, последнее - соответствует половине частоты дискретизации (частота Найквиста).
''' Значения выше половины частоты дискретизации - не использовать.
''' </returns>
Public Shared Function FFT(ByVal signal As ComplexNumber()) As ComplexNumber()

    Dim order As Integer = signal.Length 'порядок ДПФ
    CheckFftOrder(order) 'Проверяем, что порядок равен степени двойки

    Dim spectrumLen As Integer = order \ 2
    Dim j As Integer = spectrumLen

    'Бит-реверсная сортировка:
    For i As Integer = 1 To order - 2
        If (i < j) Then
            Dim tmpRe As Double = signal(j).Real
            Dim tmpIm As Double = signal(j).Imaginary
            signal(j).Real = signal(i).Real
            signal(j).Imaginary = signal(i).Imaginary
            signal(i).Real = tmpRe
            signal(i).Imaginary = tmpIm
        End If
        Dim k As Integer = spectrumLen
        Do Until (k > j)
            j -= k
            k \= 2
        Loop
        j += k
    Next

    'Цикл по уровням разложения:
    For level As Integer = 1 To CInt(Math.Log(order) / Math.Log(2))
        Dim lvl As Integer = CInt(2 ^ level)
        Dim lvl2 As Integer = lvl \ 2
        Dim tmp As Double = Math.PI / lvl2
        Dim sr As Double = Math.Cos(tmp)
        Dim si As Double = -Math.Sin(tmp)

        Dim tr As Double = 0
        Dim ur As Double = 1
        Dim ui As Double = 0
        For jj As Integer = 1 To lvl2 'Цикл по спектрам внутри уровня
            For i As Integer = (jj - 1) To (order - 1) Step lvl 'Цикл по отдельным "бабочкам"
                Dim ip As Integer = i + lvl2
                tr = signal(ip).Real * ur - signal(ip).Imaginary * ui 'Операция "бабочка"
                Dim ti As Double = signal(ip).Real * ui + signal(ip).Imaginary * ur
                signal(ip).Real = signal(i).Real - tr
                signal(ip).Imaginary = signal(i).Imaginary - ti
                signal(i).Real = signal(i).Real + tr
                signal(i).Imaginary = signal(i).Imaginary + ti
            Next
            tr = ur
            ur = tr * sr - ui * si
            ui = tr * si + ui * sr
        Next
    Next

    'Заполняем массив комплексных чисел, обработанных БПФ:
    Dim spectrum(order - 1) As ComplexNumber
    For i As Integer = 0 To order - 1 
        With signal(i)
            spectrum(i) = New ComplexNumber(.Real, .Imaginary)
        End With
    Next

    Return spectrum

End Function

''' <summary>
''' Восстанавливает сигнал по его спектру методом обратного быстрого преобразования Фурье.
''' </summary>
''' <param name="spectrum">Спектр сигнала, содержащий количество отсчётов, кратное степени двойки, и состоящий из действительной и мнимой частей.</param>
Public Shared Function InverseFFT(ByVal spectrum As ComplexNumber()) As ComplexNumber()

    Dim order As Integer = spectrum.Length 'Порядок обратного ДПФ.
    CheckFftOrder(order)

    'Изменение арифметического знака элементов мнимой части:
    For i As Integer = 0 To spectrum.Length - 1
        spectrum(i).Imaginary = -spectrum(i).Imaginary
    Next

    'Вычисление прямого БПФ:
    Dim directFFT As ComplexNumber() = FFT(spectrum)

    'Деление на order во временной области со сменой арифметического знака мнимой части:
    Dim signal(directFFT.Length - 1) As ComplexNumber
    For i As Integer = 0 To directFFT.Length - 1
        Dim ReX As Double = directFFT(i).Real / order
        Dim ImX As Double = -directFFT(i).Imaginary / order
        signal(i) = New ComplexNumber(ReX, ImX)
    Next

    Return signal

End Function