Быстрая свертка дискретных сигналов, как? Опции

[phantom]

Надо написать алгоритм КИХ фильтрации. Особенность процесса такова, что сигнал выдается фрагментами (начиная от 64 отсчетов).Светрку реализовили через БПФ в нескольких вариантах (overlap&save, overlap&add), но необходимой производительности достичь не удалось (получается в 2-8раз медленнее чем у существующих программ). Причем у программ-аналогов загрузка процессора практически не зависит от длины импульсной характеристики, а у нас - чем больше длина ИХ, тем больше времени требуется на обработку. Думали может БПФ медленное, заменили на Intel FFT которое в 10 раз быстрее общепринятых алгоритмов, но это существенно не помогло. Осталось попаробовать метод блоков переменной длины, но не совсем понятно как это работает. Что можно было бы сделать?

[thermit]

phantom:
130000 у меня есть до начала процесса обработки, FFT над ним я делаю один раз и далее все время использую при обработке каждого входного блока. Основная вычислительная работа - это умножение и суммирование в частотной области, время вычисления FFT по сравнению с ним - всего ничего...

Я решал похожую задачу (корреляция 2-х последовательностей длиной 16к отсчетов) следующим образом:
- Общую корреляцию разбивал на короткие (1к) блоки
- Вычислял корреляции этих блоков через бпф-перемножение спектров-обпф
- Результаты суммировал с перекрытием
С точки зрения числа операций выигрыша никакого. Только экономия памяти.

В вашем случае можно поступить аналогично исходя из размера входного блока:
Скажем, для 64-х входных отсчетов и длиной фильтра 130000 нужно вычислить
ceil(130000/64) = 2032 сверток. бпф-обпф размера 128 по основанию 2 займет ~ 1792 комплексных умножений-сложений. (За одно комплексное бпф вычисляется спектры 2-х вещественных последовательностей). Вычисление кроссспектра 6 вещ * 4 вещ + (если спектры не разделять) на точку. Или 390 вещ * 260 вещ +. Это явно меньше, чем затраты на бпф. Ну и на перекрытие с суммированием потребуется 63 вещ сложения.
Итого на всю свертку
2032*1792 = 3641344 комплексных */+ прямое-обратное дпф
2032*390 = 792480 вещественных * для вычисления кросспектра
2032*(260+63) = 656336 вещественных + для вычисления кросспектра и 64-х отсчетов результата.
Очевидно, что основные вычисления приходятся на преобразования. На перемножение спектров нужно не больше 5-10% всех вычислений.


Конечно, можно представить эту свертку набором совсем коротких сверток и выполнять все вычисления без дпф через алгоритмы быстрых сверток. Это даст дополнительную экономию времени только если операция умножения выполняется дольше сложения (нет аппаратного умножителя). Но скорее всего выигрыш будет незначительным...

Сообщение отредактировал thermit - Feb 1 2010, 11:34

[phantom]

Скажем, для 64-х входных отсчетов и длиной фильтра 130000 нужно вычислить
ceil(130000/64) = 2032 сверток. бпф-обпф размера 128 по основанию 2 займет ~ 1792 комплексных умножений-сложений. (За одно комплексное бпф вычисляется спектры 2-х вещественных последовательностей). Вычисление кроссспектра 6 вещ * 4 вещ + (если спектры не разделять) на точку. Или 390 вещ * 260 вещ +. Это явно меньше, чем затраты на бпф. Ну и на перекрытие с суммированием потребуется 63 вещ сложения.
Итого на всю свертку
2032*1792 = 3641344 комплексных */+ прямое-обратное дпф
2032*390 = 792480 вещественных * для вычисления кросспектра
2032*(260+63) = 656336 вещественных + для вычисления кросспектра и 64-х отсчетов результата.
Очевидно, что основные вычисления приходятся на преобразования. На перемножение спектров нужно не больше 5-10% всех вычислений.

Спасибо за ответ по сути!
Однако же, Вы не совсем правы
Во-первых, прямое БПФ размера 2*М, где М-длина блока, у меня вычисляется всего один раз - для каждого входного блока, т.к. БПФ импульсной характеристики я считаю перед началом процесса (импульсная характеристика разбивается на перекрывающиеся блоки длины 2*М).
Во-вторых, обратное БПФ того же размера у меня вычисляется, как ни удивительно это звучит, тоже всего один раз - за счет суммирования перекрывающихся результатов сверток в частотной области (см. статью: http://www.ramsete.com/Public/Presentations/mohonk_2003.pdf)
Таким образом, для каждого входного блока имеем при М=64:
1) БПФ+ОБПФ размера 128 по основанию 2: 1792 комплексных умножений-сложений
2) Умножение спектров: 2032*128*(6 вещ *, 4 вещ+) = 1560576*, 1040384+
3) Суммирование перекрытий в частотной области: 2031*128=259968+
Ну так и сколько там того БПФ?

Пояснение. Схема overlap-discard для частичных сверток входного блока с перекрывающимися кусочками ИХ используется "внутри" общей схемы overlap-add для свертки входных блоков с целой ИХ (блоки входного сигнала берутся без перекрытия).

ЗЫ. Схема работает (это если вдруг у кого появятся подозрения что это мухлеж и нельзя обойтись одной парой коротеньких БПФ-ОБПФ на блок )