Размещено на http://www.allbest.ru/

Быстрая свертка на основе быстрого преобразования Фурье

Рассмотрены практические вопросы использования быстрого преобразования Фурье комплексных сигналов длины N для повышения эффективности вычислений дискретной свертки за счет одновременной обработки двух действительных последовательностей каждая длины N. Представлены условия выбора минимального периода для кольцевой свертки, и приведена оценка погрешности вычисления свертки из-за сокращения длины ядра. Оценена сложность реализации свертки в элементарных операциях. Показаны преимущества вычисления длинных сверток в частотной области.

Процесс томографической реконструкции объекта по его проекциям связан с большим объемом вычислений, среди которых свертка занимает заметное место. Оптимизация вычислений свертки важна как для томографической реконструкции, так и практически для любой области цифровой обработки сигналов [1].

Рассмотрим случай одномерной свертки сигнала (исходной проекции) с соответствующим ядром. В терминах линейных систем ядро представляет собой импульсный отклик системы на -функцию и характеризует ее динамические свойства. Ядро свертки - функция с четной симметрией.

Пусть сигнал представляет собой входную последовательность длиной N1 отсчетов. Ядро содержит нечетное число отсчетов N2. При обработке входной последовательности длиной N1 отсчетов число отсчетов полного ядра равно величине N2 = 2N1 - 1. При апериодической свертке [1] двух конечных последовательностей длиной N1 и N2 отличных от нуля отсчетов свертки будет не более

M = N1 + N2 - 1 = 3N1 - 2.

Апериодическая свертка может быть заменена периодической (круговой или циклической), если период повторения сигнала будет, по крайней мере, столь же большим как длина отсчетов результата свертки (3N1 - 2). При таких условиях не происходит никакого наложения сигналов, что облегчает реализацию качественной интерполяции. На практике требуемое число отсчетов результатов свертки обычно равно длине входной последовательности N1. Это позволяет несколько сократить минимальный период круговой свертки. Проведенный анализ показал, что при четной симметрии ядра и нечетном числе его членов величина M может быть выбрана из менее жесткого условия:

M N1 + (N2 - 1)/2 = < 2N1 - 1.

сигнал дискретный свертка

При этом возникает некоторое наложение на отдельных участках результатов свертки. Тем не менее, интересующий участок круговой свертки длиной N1 остается без влияния наложений и содержит правильные результаты. Это важно, потому что мы можем примерно в 1,5 раза сократить минимальный период круговой свертки, получая при этом правильный результат для N1 интересующих значений апериодической свертки.

Для сверток с коротким ядром N2 < 2N1 - 1 могут возникнуть проблемы с использованием интерполяции с ограниченным спектром (например, sinc-интерполяции). Дело в том, что ход полученной после такой интерполяции плавной кривой, проходящей через все интересующие точки N1, будет также зависеть и от остальных M - N1 точек, претерпевших искажения от наложения сигналов из-за сокращения периода круговой свертки.

Известно, что трудоемкая операция дискретной свертки во «временной» области (области сигнала) эквивалентна в частотной области простому умножению одноименных гармоник спектров сигнала и ядра. Однако при этом необходимо дополнительно выполнить переход из области сигналов в область спектров и, после умножения одноименных гармоник спектров, дополнительно сделать обратный переход в область сигналов. Переход из временной области в частотную и обратно требует прямого дискретного преобразования Фурье (ДПФ) входных последовательностей по M отсчетов, умножения одноименных гармоник спектра с последующим обратным преобразованием результатов их произведения во временную область с помощью обратного ДПФ (ОДПФ). Теоретические основы вычисления свертки в частотной области известны весьма давно. Прошло почти два столетия с тех пор, как Жан Батист Жозеф Фурье (1807 г.) представил доклад о возможности представления непрерывного периодического сигнала в виде суммы синусоид.

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

Практическое применение вычисления свертки с переходом в область спектров имеет смысл только при условии применения эффективных методов ДПФ и ОДПФ. В этом плане особого внимания заслуживает быстрое преобразование Фурье (БПФ), представляющее собой частный случай ДПФ при числе отсчетов N = pq, где p и q - целые числа. В настоящее время, с полным на то основанием, БПФ связывают с именем Кули и Тьюки (1960 г.) [2]. (Справедливости ради следует отметить что, пожалуй, первым, кто использовал свойства ДПФ для усовершенствования алгоритма его расчета, столетием ранее, был немецкий математик Гаусс. Свой вклад внесли Рунге (1903 г.), Ланцош (1942 г.) и др.).

Возможность применения БПФ длинны N для вычисления круговой свертки вместо ДПФ длины M, ограничена единственным требованием:

N ? M

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

M 3N1 - 2

При вычислении ДПФ методом БПФ по основанию 2 достаточно увеличить период до N = 2q ? M. В общем, для выполнения свертки необходимо выполнить следующие операции.

Дополнить последовательности сигнала N1 и ядра N2 нулевыми отсчетами до общей длины периода N.

Выполнить прямое комплексное БПФ сигнала длиной N отсчетов (мнимая составляющая равна нулю).

Выполнить прямое комплексное БПФ ядра длиной N отсчетов (мнимая составляющая равна нулю).

Произвести почленное умножение одноименных гармоник спектров сигнала и ядра. Для этого в общем случае необходимо выполнить N комплексных умножений.

Выполнить обратное БПФ (ОБПФ) результатов N комплексных умножений.

Выбрать из N результатов действительной части ОБПФ N1 значений свертки (мнимая составляющая результатов ОБПФ очевидно равна нулю).

В результате выполнения пунктов 1-6, формально используя теорему об эквивалентности свертки во временной области умножению в частотной области, можно получить апериодическую свертку исходного сигнала N1 с ядром N2 = = 2N1 - 1.

Оптимизация вычислений свертки. Рассмотрим пути повышения эффективности вычисления свертки. Быстрое преобразование Фурье дает возможность преобразования комплексного сигнала длины N, представляющего собой пару последовательностей (действительной и мнимой), длиной N отсчетов каждая [3]. Это открывает возможность одновременной обработки двух действительных последовательностей [4]. Но так как программа БПФ предполагает наличие действительной и мнимой составляющей, то в данном случае на выходе вместо ожидаемых спектров пары входных действительных последовательностей мы получим «нечто» другое. В принципе, если нас интересуют спектры, то путем некоторых дополнительных операций их можно получить. При вычислении же свертки в этом нет нужды. Оказалось, что достаточно выполнить почленное умножение комплексных гармоник спектров «нечто» и ядра. Из области спектров при помощи ОБПФ N комплексных результатов умножения преобразуется в две круговых свертки двух входных сигналов с ядром. Такой подход позволяет примерно вдвое сократить число операций на вычисление свертки.

Есть также возможность несколько сократить объем вычислений при реализации БПФ за счет использования прямого порядка обработки отсчетов со стороны сигнала (т.н. «прореживание по частоте»), так как при выполнении покомпонентного умножения спектров порядок следования гармоник не имеет значения. Это позволяет исключить операции, связанные с созданием требуемого порядка следования членов обрабатываемой последовательности.

Еще одна возможность сокращения количества операций - переход от операций комплексных умножений к действительным умножениям. Переход на умножении спектров позволяет в 6 раз уменьшить требуемое количество операций с плавающей точкой. Условие правомочности перехода к действительным умножениям спектров - отсутствие мнимой составляющей в спектре ядра. Последнее, учитывая четность функции ядра (h[i] = h[-i]), может быть достигнуто за счет соответствующего его кольцевого сдвига (на 0 или N/2 отсчетов). Рассмотрим пример. Пусть длина периода N = 8, а исходное ядро задано в 5-ти дискретных точках:

сигнал дискретный свертка

При сдвиге на N/2 исходные 5 отсчетов следует дополнить нулями слева и справа, причем центральный отсчет исходной последовательности h[0] должен быть расположен в точке N/2 (N/2 = 4) участка 0…N - 1:

Показанные выше оба расположения исходных отсчетов ядра на кольце исключают появление в спектре ядра мнимой составляющей. Заметим, что оба расположения эквивалентны с точки зрения результатов свертки. Однако если не предполагается последующая интерполяция свертки, удобнее воспользоваться расположением (5a). В этом случае все N1 результатов свертки идут первыми. Полученный после одного прямого БПФ спектр ядра запоминается и далее при выполнении свертки используется для всего массива данных. Рассмотрим на простом примере вычисление свертки двух входных последовательностей длиной N1 = 4 ([0, 0, 0, 1] и [1, 0, 0, 0]) с ядром длиной N2 = 7, ([-3, -7, -35, 105, -35, -7, -3]). После дополнения нулями входных последовательностей до одинаковой длины (N1 + (N2 - 1)/2 ? 2P = 8) получим:

Ядро может быть представлено двумя различными способами, при которых его спектр не содержит мнимой составляющей. В нашем случае (N = 8) возможны два эквивалентных представления ядра:

Единственно, чем отличаются указанные последовательности (h), так это кольцевым сдвигом по отношению к входным последовательностям (sx, sy).

Две входные действительные последовательности (sx, sy), каждая длиной N отсчетов, рассматриваем как одну комплексную последовательность длиной N. После выполнения комплексного БПФ длиной N получим действительную (Sx) и мнимую (Sy) составляющие спектра. Выполняем их покомпонентное умножение на полученный ранее действительный спектр ядра (Hr):

После выполнения обратного БПФ результатов умножения спектров во временной области получим значения двух кольцевых сверток. Результаты кольцевых сверток с ядром h сигналов sx и sy будут расположены соответственно на месте действительной и «мнимой» составляющих результата обратного БПФ.