Рейтинг@Mail.ru

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

Print Friendly, PDF & Email

Даётся программный код для прямого и обратного преобразования Фурье. Рассматривается быстрое преобразование Фурье.

Дискретное преобразование Фурье (ДПФ) – это мощный инструмент анализа, который широко используется в области цифровой обработки сигналов (ЦОС). Существуют прямое и обратное преобразования Фурье. Прямое дискретное преобразование Фурье переводит сигнал из временной области в частотную и служит для анализа частотного спектра сигнала. Обратное преобразование делает ровно противоположное: по частотному спектру сигнала восстанавливает сигнал во временной области.

Для расчёта преобразования Фурье обычно используется ускоренная процедура расчёта – т.н. быстрое преобразование Фурье (БПФ). Это позволяет в значительной мере сократить процессорное время на достаточно сложные и ресурсоёмкие математические расчёты.

1Комплексныечисла

Для начала нам потребуется вспомогательный класс, который будет описывать комплексные числа. Комплексные числа – это особый вид чисел в математике. Каждое комплексное число состоит из двух частей – действительной и мнимой. Сейчас нам достаточно знать о комплексных числах применительно к ДПФ то, что действительная часть комплексного числа хранит информацию об амплитуде сигнала, а мнимая – о фазе.

Код класса для описания комплексных чисел (разворачивается)
''' <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

2Прямое дискретное быстрое преобразование Фурье

Теперь приведу код функции, которая реализует расчёт прямого дискретного преобразования Фурье. Здесь встречается термин «бабочка». Не пугайтесь, это общепринятое название одной из элементарных операций, которые входят в алгоритм БПФ. Подробнее с термином «бабочка» можно ознакомиться здесь.

На вход функции передаётся массив комплексных чисел. Действительная часть которого представляет произвольный дискретный сигнал, с отсчётами через равные промежутки времени. Мнимая часть содержит нули. Число отсчётов в сигнале должно равняться степени двойки. Если ваш сигнал короче, то дополните его нулями до числа, кратного степени 2: 256, 512, 1024 и т.д. Чем длиннее сигнал, тем у рассчитанного спектра будет выше разрешение по частоте.

Код для расчёта прямого быстрого преобразования Фурье на VB.NET (разворачивается)
''' <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

3Обратное дискретное быстрое преобразование Фурье

Обратное дискретное преобразование Фурье (ОДПФ) одним из этапов расчёта включает в себя прямое ДПФ на массиве комплексных чисел, где мнимая часть – это инверсия относительно оси X мнимой части спектра.

Код для расчёта обратного быстрого преобразования Фурье на VB.NET (разворачивается)
''' <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

Ну и конечно же, опишем использовавшийся метод, который проверяет число элементов переданного массива:

''' <summary>
''' Проверяет, является ли порядок БПФ степенью двойки, и если нет - вызывает исключение.
''' </summary>
''' <param name="order">Порядок БПФ.</param>
Private Shared Sub CheckFftOrder(ByVal order As Integer) 
    Dim chk As Double = Math.Abs(Math.Floor(Math.Log(order, 2)) - Math.Log(order, 2))
    If (chk > 0.0001) Then
        Throw New ArgumentException(String.Format("Длина массива ({0}) не кратна степени двойки.", order))
    End If
End Sub

4Проверка прямого и обратного преобразования Фурье

Теперь давайте проверим, что наши функции работают. Для этого пропустим произвольный сигнал через механизм прямого преобразования Фурье, а затем «соберём» его обратно с помощью обратного преобразования Фурье. Восстановленный сигнал должен практически совпадать с исходным. Ошибки округления, возникающие при работе с числами в компьютере, имеют место быть, поэтому сигналы не будут идентичны полностью, но их отклонение друг от друга должно быть пренебрежимо малым.

Для примера в качестве исходного сигнала возьмём функцию синуса и сформируем данные длиной 128 отсчётов вот таким образом:

Dim cn(127) As ComplexNumber
For i As Integer = 0 To cn.Length - 1
    cn(i) = New ComplexNumber(Math.Sin(i * 3 * Math.PI / 180))
Next

Получим вот такой сигнал:

Исходный сигнал во временной области до преобразования Фурье
Исходный сигнал во временной области до преобразования Фурье

Здесь по оси X – номера отсчётов во временной области, по оси Y – амплитуда. Обратим внимание, что сигнал состоит только из действительных частей, а мнимая часть на всём отрезке равна "0".

Теперь передадим этот сигнал на вход функции FFT(). По полученным в ходе прямого преобразования Фурье массивам комплексных чисел построим два графика – действительной (Re) и мнимой (Im) частей спектра:

Действительная и мнимая части сигнала в частотной области
Действительная и мнимая части сигнала в частотной области

Здесь по оси X – отсчёты в частотной области, по оси Y – амплитуда. Чтобы получить реальные значения частоты, необходимо рассчитать их, учитывая, что "0" оси Y соответствует нулевой частоте, максимум оси Y соответствует частоте дискретизации.

Полученный спектр сигнала передадим функции обратного преобразования Фурье IFFT(). Получим массив комплексных чисел, где действительная часть будет содержать восстановленный сигнал:

Действительная и мнимая части восстановленного с помощью обратного преобразования Фурье сигнала
Действительная и мнимая части восстановленного с помощью обратного преобразования Фурье сигнала

Как видно, восстановленный сигнал полностью повторяет исходный.

Последнее изменениеПятница, 11 Июнь 2021 12:25 Прочитано 8398 раз
Теги :

Поблагодарить автора:

Поделиться

Print Friendly, PDF & Email

Оставить комментарий