Vadim Kudryavtsev home page

 

Алгоритмы цифровой обработки сигналов

В этом разделе я буду выкладывать описания алгоритмов и их реализации. Постараюсь обязательно выкладывать их реализации на Матлабе. Если они реализованы на каком-либо языке программирования (то есть ближе к железу) - постараюсь выкладывать тоже. Пока этот раздел несистематизирован. С ростом количества алгоритмов придет и необходимость в систематизации.


Решение тёплицевой системы вида y=Px
R. Kumar. A fast algorithm for solving a toeplitz system of equations. IEEE Trans. on Acous., Speech, and Sign. Proc., 33:254--267, February, 1985

Page 1-2 Page 3-4 Page 5-6 Page 7-8 Page 9-10 Page 11-12 Page 13-14

Реализация на МатЛабе

Данный алгоритм является быстрой реализацией решения уравнения y=Px, где P является теплицевой матрицей. Такие уравнения часто приходится решать, например, при нахождении инверсных фильтров в задачах эквализации трактов. Сложность данного алгоритма - O(n*log2(n)^2). Сравните эту оценку со сложностью прямого инвертирования матрицы O(n^3) и со сложностью алгоритма Левинсона-Дурбина O(n^2).

Алгоритм состоит из трех шагов.
Шаг 1:
- дополнение теплицевой матрицы до матрицы-циркулянта и нахождение матрицы, инверсной к данной матрице-циркулянту. Используется преобразование Фурье.
Шаг 2:
- находится первая строка инверсной матрицы P^(-1). Используется алгоритм деления полиномов.
Шаг 3:
- находится решение уравнения - вектор x. Используется знание первой строки инверсной матрицы и особая структура матрицы P.

Коментарии:
Для меня не очевиден второй шаг. В статье в качестве материала, содержащего теоретическое обоснование, приводятся статьи
[4] W.F. Trench. "An algorithm for the inversion of finite Toeplitz matrices", J.SIAM, vol.12, pp.512-522, Sept. 1964.
[5] S. Zohar. "The solution of a Toeplitz set of linear equation", J. Ass. Comput. Mach., vol.21, pp.272-276, Apr. 1974.
[15] H.H. Rosenbrock. State Space and Multivariable Theory. New-York: Wiley, 1970.
 


КИХ-фильтрация на основе бэнк-фильтров
M.Vetterli. Running FIR and IIR Filtering Using Multirate techniques. IEEE Trans. On Acoustics, Speech and Signal Processing, 36:730-738, May 1988.

Page 1-2 Page 3-4 Page 5-6 Page 7-8 Page 9-10

Неплохую реализацию на Матлабе можно посмотреть в отчете Jelen Tesic

Краткий комментарий:
Как известно, при длинах фильтров меньше 100 спектральные методы вычисления сверток не имеют преимущества или проигрывают по вычислительной сложности прямым методам. В этом случае имеет смысл пользоваться алгоритмами, основанными на китайской теореме об остатках [1] Однако, при длине фильтра, равной несколько десятков, хотя и возможно построить быстрый алгоритм методом вложений (nested algorithms), получившийся алгоритм является замысловатым и вдобавок при изменении длины фильтра его приходится переделывать.
Идея предлагаемого алгоритма состоит в разбиении импульсной характеристики фильтра и входной последовательности на полифазные компоненты. Свертку можно выразить через комбинацию (псевдоциклическую свертку) этих полифазных компонент. Быстрый алгоритм получается при помощи использования китайской теоремы об остатках.
Преимущества алгоритма:
1. Свертка делится на несколько сверток меньшей длины, причем каждая свертка обрабатывается на пониженной тактовой частоте.
2. За счет этого общее количество вычислений может падать в 25%.
Недостатки:
1. Больший размер кода.
2. Больший размер данных.
Послесловие:
1. В оригинальной статье Мартина Веттерли не указано использование китайской теоремы об остатках, поэтому пришлось указать на это в своей статье "Теория анализа и синтеза бэнк-фильтров и их применение"
2. Если Вы интересуетесь развитием этой идеи на случай конечных полей, можете также ознакомится со статьями Вычисление линейной свертки с использованием МТЧП в структуре бэнк-фильтра" и "Вычисление линейной свертки с использованием полиномиального преобразования в структуре бэнк-фильтра"
3. Спустя много лет после публикации своей статьи (см. пункт 1), в 2005 году я случайно обнаружил, что ту же самую идею двумя годами раньше высказал профессор Митра в статье
Ing-Song Lin, and Sanjit K. Mitra, "Overlapped Block Digital Filtering", IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS-II, VOL. 43, NO. 8, AUGUST 1996
Признаю его первенство в этом вопросе.


Быстрый алгоритм 50-точечного преобразования Фурье

Описание алгоритма

Реализация алгоритма на Матлабе

“На самом же деле хорошие БПФ-алгоритмы существуют практически для произвольной длины блока” (Ричард Блейхут)

Всем известна быстрая реализация дискретного преобразования Фурье по основанию, кратному степени двойки, с помощью алгоритма Кули-Тьюки. Разумеется, это один из наиболее эффективных и однотипно структурированных алгоритмов (сам часто выбираю преобразования равными степени двойки). Однако, приводимый здесь 50-точечный алгоритм быстрого преобразования Фурье очень наглядно иллюстрирует вышеприведенное высказывание Блейхута. Путем надлежащего использования методов разбиения многоточечных преобразований на короткие вычислительная сложность алгоритма по любому основанию может быть существенно снижена.
При решении этой задачи (кстати, возникшей на форуме телесистем) использовались большинство способов снижения вычислительных затрат, описанных в книге Р.Блейхута "Быстрые алгоритмы цифровой обработки сигналов".


Комбинационный алгоритм вычисления квадратного корня
Описание алгоритма Блок-схема алгоритма

Реализация алгоритма на языке Verilog

Я даже не подозревал, что квадратный корень можно вычислять так быстро и с такими низкими затратами. Результат вычисляется за n/2 итераций, где n - количество двоичных бит, которым может быть представлено исходное число. Причем на каждой итерации из арифметических операций используется только одна операция вычитания и одна операция сравнения. Остальное - это сдвиги и операции конкатенации. Поэтому этот алгоритм очень подходит для реализации на ПЛИС.

Аналогичные алгоритмы разработаны для большинства распространенных математических функций.
Рекомендую в этом направлении книги
:

Также даю ссылку на очень интересную статью
Auger, F. Zhen Lou Feuvrie, B. Feng Li. Multiplier-Free Divide, Square Root, and Log Algorithms [DSP Tips and Tricks]
Матлаб файлы к ней можно скачать по этому адресу.


Дискретные тригонометрические преобразования
Быстрый способ вычисления дискретных косинусных и синусных преобразований друг через друга и через БПФ
Markus Puschel and Jose Moura The Algebraic Approach to the Discrete Cosine and Sine Transforms and their Fast Algorithms SIAM Journal of Computing 2003, Vol. 32, No. 5, pp. 1280-1316
R. Gluth, “ Regular FFT-related transform kernels for DCT/DST- based polyphase filter banks,” Proc. IEEE ICASSP 1991, pp.2205-2208, Toronto, Canada, May 1991.
Математические соотношения между дискретными тригонометрическими преобразованиями (файл в формате Word)
Матлаб код, иллюстрирующий связь некоторых преобразований

Как известно, существуют 8 видов дискретного косинусного и столько же видов дискретного синусного преобразований. В указанной по ссылке статье Маркуса Пюшеля и Хосе Моуры разработана красивая теория, следствием которой является очень простая связь косинусных и синусных преобразования вида 1...4 между собой, а также аналогичная связь косинусных и синусных преобразования вида 5...8. Дополнительно в статье Рольфа Глюта (см. вторую ссылку) указан эффективный способ вычисления преобразований вида 4 (косинусного и синусного) через преобразование Фурье половинной длины. Таким образом, косинусные и синусные преобразования вида 1...4 и преобразование Фурье оказываются связанными строгими математическими соотношениями, позволяющими найти эффективный способ вычисления одного преобразования через другое. Иными словами, если в DSP-библиотеке есть код только одного из этих преобразований (как правило в наличии имеется только БПФ), мы можем вычислить любое другое ценой небольшого увеличения вычислений на пред- и постпреобразование результата.
Если вас более глубоко интересует этот вопрос, вы можете заглянуть на страничку Маркуса Пюшеля.

Соотношения между дискретными тригонометрическими преобразованиями
Отчеты сгенерированы на матлабе

DCTI
DCTI<->DSTII
DCTI<->DCTII
DCTI<->DSTI

DCTIII
DCTIII<->DCTIV
DCTIII<->DCTIV(alt)

DCTV
DCTV<->DSTV
DCTV<->DCTVI
DCTV<->DSTVI

DCTVII
DCTVII<->DSTVII
DCTVII<->DCTVIII
DCTVII<->DSTVIII

DSTIII
DSTIII<->DSTIV
DSTIII<->DSTIV(alt)

DSTI
DSTI<->DSTII
DSTI<->DCTII
DSTI<->DCTI

DSTVII
DCTVII<->DSTVII
DSTVII<->DCTVIII
DSTVII<->DSTVIII

DSTV
DCTV<->DSTV
DSTV<->DCTVI
DSTV<->DSTVI

DCTVI
DCTV<->DCTVI
DSTV<->DCTVI

DCTVIII
DCTVII<->DCTVIII
DSTVII<->DCTVIII

DCTII
DCTII<->DCTIV
DCTII<->DCTIV(alt)
DCTII<->DSTII
DCTII<->DCTI
DCTII<->DSTI

DCTIV
DCTIV<->DCTIII
DCTIV<->DCTIII(alt)
DCTIV<->DCTII
DCTIV<->DCTII(alt)

DSTVIII
DCTVII<->DSTVIII
DSTVII<->DSTVIII

DSTVI
DCTV<->DSTVI
DSTV<->DSTVI

DSTIV
DSTIV<->DSTIII
DSTIV<->DSTIII(alt)
DSTIV<->DSTII
DSTIV<->DSTII(alt)

DSTII
DSTII<->DSTIV
DSTII<->DSTIV(alt)
DSTII<->DCTII
DSTII<->DCTI
DSTII<->DSTI

Дополнительные плюшки:


EMQF фильтры
Иллюстрации применений EMQF фильтров
Matlab файлы, используемые для расчета EMQF фильтров

EMQF фильтры (elliptic filters with minimal Q-factor) являются частным случаем эллиптических фильтров, обладающих одинаковой неравномерностью в полосе пропускания и в полосе заграждения. Несмотря на это ограничение, данный класс фильтров обладает интересными свойствами, на основе которых можно получить фильтры, быстрые в реализации (на основе параллельного соединения всепропускающих фильтров), получить эффективные структуры БИХ дециматоров и интерполяторов, а также пар БИХ фильтров, сопряженных по Гильберту. Помимо демонстрационной странички, иллюстрирующей эти и другие применения EMQF фильтров, Вы можете посетить сайт Miroslav D. Lutovac: список публикаций , где сможете более глубоко познакомиться с данной темой.


Расчет БИХ-фильтров при помощи оптимизации по методу средних квадратов

Иллюстрации использования программы

Matlab файлы, используемые для расчета фильтров

Я также восстановил все матлаб-скрипта из диссертации мистера Ланга. Вы можете найти их по ссылке diss_lang.zip.

Данный пакет предназначен для расчета БИХ-фильтров по заданным амплитудной и фазовой характеристикам. При расчете используется оптимизация по методу средних квадратов.  Одним из параметров является максимальный радиус полюсов фильтра. Это позволяет гарантировать устойчивость фильтра при дискретизации коэффициентов, а также контролировать динамический диапазон переменных внутреннего состояния фильтра. Пакет заимствован из диссертации Mathias Lang, "Algorithms for the Constrained Design of Digital Filters with Arbitrary Magnitude and Phase Responses", Vienna, June 1999.


Симметрично-сопряженное секвентивно-упорядоченное преобразование Адамара

Иллюстративная страничка CS-SCHT

Matlab пакет для применения CS-SCHT

Этот пакет был написан по статье "Conjugate Symmetric Sequency-Ordered Complex Hadamard Transform", Aye Aung   Boon Poh Ng   Rahardja, S.. Авторы утверждают, что быстрый алгоритм для вещественного варианта этого преобразования требует на 25% меньше операций, чем традиционное преобразование Адамара. К тому же базисные функции этого преобразования являются аппроксимацией базисных функции ДПФ, что делает его привлекательным для быстрой (хотя и неточной) оценки спектра сигнала. По ссылке на иллюстративную страничку представлены свойства этого преобразования и применение его к сжатию изображений, а также сравнение качества сжатия при использовании различных преобразований.



Filter visualization tool (fvtool) для Octave

Иллюстративная страничка для fvtool for Octave

Матлаб скрипт (fvtool.m и demo-script)


Иногда я использую Octave (open source аналог Матлаба) на своюх мобильных устройствах. Отсутствие такого удобного инструмента как fvtool доставляет определенные неудобства. Поэтому однажды я написал скрипт fvtool for Octave. Он не такой мощный как в Матлабе (в Матлабе он интерактивный), но с его помощью вы сможете построить различные характеристики: частотную, фазовую, фазовой задержки, групповой задержки, а также импульсную характеристику. Частотная и амплитудная оси отображаться в линейном и логарифмическом масштабе. Для удобства анализаможет быть указана частота дискретизации. В этом случае оси будут отображать значения, соответствующие данной частоте дискретизации.


Использование БПФ длины N для блочной свертки по методу перекрытия с суммированием и перекрытия с накоплением

Иллюстративная страничка для метода перекрытия с суммированием

Иллюстративная страничка для метода перекрытия с накоплением

Матлаб реализации предложенных методов


Классический методы блочной свертки "перекрытие с суммированием" и "перекрытие с накоплением" требуют преобразования Фурье длины 2N. Предлагаемые методы позволяют добиться уменьшения количества вычислений путем замены преобразования Фурье длины 2N на собственное преобразование относительно косой свертки mod(x^N-j). Упомянутое собственное преобразование вычисляется последовательным применением умножения на диагональную матрицу поворачивающих множителей exp(j*pi/(2N)*k), k=0..N-1, с последующим преобразованием Фурье длины N, Таким образом, удаётся значительно снизить количество вычислений, необходимых для блочной свертки. Оригинальная идея была предложена в статье Jung Gap Kuk, Seyun Kim, and Nam Ik Cho, "A New Overlap Save Algorithm for Fast Block Convolution and Its Implementation Using FFT," Journal of Signal Processing Systems(Online), Feb 2010.