Свертка через БПФ основано на том, что перемножение в частотной области соответствует свертке во временной области. Входной сигнал преобразуется с помощью ДПФ в частотную область, перемножается на частотный отклик фильтра, и затем преобразуется назад во временную область с помощью инверсного ДПФ. Этот подход был известен еще со времен Фуре, но на практике не применялся. Это объясняется тем, что время вычисления ДПФ было больше, чем время прямого вычисления свертки. Ситуация изменилась с 1965 года с открытием быстрого преобразования Фурье (БПФ). Использование алгоритмов БПФ сделало вычисление ДПФ, выполненное через частотную область, более быстрым, чем прямое вычисление свертки во временной области. Окончательный результат тот же самый; изменилось только число вычислений благодаря более эффективным алгоритмам. По этой причине свертка через БПФ так же называется высокоскоростной сверткой.
РИСУНОК 18-1
Метод обработки сигнала по сегментам (overlap-add method). Цель – осуществление свертки входного сигнала (а) с ядром фильтра (b). Входной сигнал разбивается на сегменты (с), (d) и (е), добавляется необходимое для осуществления свертки число нулей. Свертка каждого входного сегмента с ядром фильтра производит выходные сегменты (f), (g) и (h). Выходной сигнал (i) находится сложением перекрывающихся выходных сегментов.
Свертка через БПФ (дальше просто БПФ свертка) использует метод обработки сигнала по сегментам, показанный на рисунке 18-1, но только способ преобразования входного сегмента в выходной изменен. Рисунок 18-2 показывает, как входной сегмент преобразуется в выходной сегмент с помощью БПФ свертки. Вначале находится частотная характеристика фильтра с помощью ДПФ ядра фильтра (используя БПФ). Например, (а) показывает ядро фильтра, это низкочастотный оконный sinc фильтр. БПФ преобразует его в реальную и мнимую часть частотной характеристики (частотного отклика), показанных в (b) и (с). Эти сигналы частотной области могут не выглядеть похожими на характеристику низкочастотного фильтра потому, что они показаны в прямоугольной форме. Вспомните, полярная форма обычно лучше подходит для человеческого восприятия, а прямоугольная - для математических вычислений. Эти реальная и мнимая части хранятся в памяти компьютора для использования при вычислении каждого сегмента.
Рисунок (d) показывает обработку входного сегмента. БПФ используется для нахождения частотного спектра, он показан в (е) и (f). Частотный спектр выходного сегмента (h) и (i) находится перемножением частотного отклика фильтра (b) и (с) на спектр входного сегмента (е) и (f). Поскольку эти спектры состоят из реальной и мнимой части, они перемножаются согласно формуле 9-1 из главы 9. Затем используется инверсное БПФ для нахождения выходного сегмента (g) из его частотного спектра (h) и (i). Выходной сегмент точно такой же, как если бы он был получен прямой сверткой входного сегмента (d) и ядра фильтра (а).
БПФ должен иметь достаточную длину, чтобы не было циклической свертки (тоже обсуждалось в главе 9). Это означает, что БПФ должно иметь ту же длину, что и выходной сегмент (g). Например, на рисунке 18-2 ядро фильтра содержит 129 точек, и каждый сегмент содержит 128 точек, это делает длину выходного сегмента равной 256 точкам. В этом случае необходимо использовать 256-точечное БПФ. Следовательно, к ядру фильтра (а) надо добавить 127 нулей для получения общей длины в 256 точек. Аналогично, к каждому входному сегменту надо добавить 128 нулей. Как другой пример, представьте, что вам необходимо свернуть очень длинный сигнал с ядром фильтра, имеющего 600 отсчетов. Одна возможность заключается в использовании сегментов по 425 точек и 1024-точечного БПФ. Другая – использование сегмента из 1449 точек и 2048-точечного БПФ.
Таблица 18-1 показывает пример программы выполнения БПФ свертки. Эта программа фильтрует 10 миллионов точек сигнала, путем свертки его с 400 точками ядра фильтра. Она выполняется разбивкой сигнала на 16 000 сегментов по 625 точек. Когда каждый из этих сегментов свертывается с ядром фильтра, получается выходной сегмент из 625+400-1=1024 точек. Таким образом, требуется 1024-точечное БПФ. После определения и инициализации всех массивов (линии 130-230), на первом шаге вычисляется и запоминается частотная характеристика фильтра (линии 250-310). Линия 260 вызывает воображаемую подпрограмму для загрузки ядра фильтра в XX[0]…XX[399] и устанавливает XX[400]…XX[1023] в ноль. Подпрограмма в линии 270 – БПФ, преобразующее 1024 точек, содержащихся в XX[ ], в 513 отсчетов, содержащихся в REX[ ] и IMX[ ] - реальной и мнимой частях частотной характеристики. Эти величины передаются в массивы REFR[ ] и IMER[ ] (для REal и IMaginary Frequency Response) для последующего использования в программе.
РИСУНОК 18-2
БПФ свертка. Ядро фильтра (а) и сегмент сигнала (d) преобразуются в свои спектры (b), (c) и (e), (f) через БПФ. Эти спектры перемножаются и дают спектр выходного сегмента (h), (i). Затем через инверсное БПФ находится выходной сегмент (g).
Цикл FOR-NEXT между линиями 340 и 580 контролирует обработку 16 000 сегментов. В линии 360 подпрограмма загружает следующий сегмент для обработки в XX[0]…XX[624] и устанавливает XX[625]…XX[1023] в ноль. В линии 370 подпрограмма БПФ используется для нахождения частотного спектра сегмента с реальной частью из 513 точек, расположенной REX[ ], и с мнимой частью из 513 точек, расположенных в IMX[ ]. Линии 390…430 показывают перемножение частотного спектра сегмента, содержащегося в REX[ ] и IMX[ ] на частотный отклик фильтра, содержащийся в REFR[ ] и IMER[ ]. Результат перемножения записывается в REX[ ] и IMX[ ] на место предыдущих данных. Поскольку это теперь частотный спектр выходного сегмента, может быть использовано инверсное БПФ для нахождения выходного сегмента. Это выполняет воображаемая подпрограмма инверсного БПФ в линии 450, которая преобразует 513 точек, содержащихся в REX[ ] и IMX[ ], в 1024 точки, содержащихся в XX[ ], выходном сегменте.
Линии 470…550 управляют перекрытием сегментов. Каждый выходной сегмент делится на две секции. Первые 625 точек (от 0 до 624) необходимо скомбинировать с перекрывающимися точками из предыдущего выходного сегмента и записать в выходной сигнал. Последние 399 точек (от 625 до 1023) необходимо запомнить, так как они могут перекрываться со следующим выходным сегментом.
Что бы это понять, посмотрите назад на рисунок 18-1. Отсчеты от 100 до 199 в (g) необходимо скомбинировать с перекрывающимися отсчетами из предыдущего выходного отсчета (f), и затем их можно поместить в выходной сигнал (i). Для сравнения, отсчеты от 200 до 299 в (g) необходимо запомнить, что бы их можно было скомбинировать с отсчетами из следующего сегмента (h).
Вернемся к программе. Массив OLAP[ ] используется для хранения 399 отсчетов, которые перекрываются с отсчетами следующего сегмента. В линиях 470…490 399 величин в этом массиве (из предыдущего выходного сегмента) складываются с отсчетами текущего сегмента, содержащиеся в XX[ ]. Мифическая подпрограмма в линии 550 направляет выходные 625 отсчетов от XX[0] до XX[624] в файл, хранящий выходной сигнал. 399 отсчетов текущего выходного сегмента, которые необходимы для следующего выходного сегмента, записываются в OLAP[ ] в линиях 510…530.
После обработки всех сегментов от 0 до 15 999, массив OLAP[ ] будет содержать 399 отсчетов из сегмента 15999, которые должны перекрываться с сегментом 16 000. Поскольку сегмент 16 000 не существует (или может рассматриваться как состоящий из одних нулей), эти 399 отсчетов записываются в выходной сигнал в линии 600. Это делает длину выходного сигнала равной 16000´625+399=10 000 399 точкам. Это совпадает с длиной входного сигнала плюс длина ядра фильтра минус единица.
