Алгоритм быстрого преобразования Фурье.
БПФ с прореживанием по времени.
Рассмотрим идею БПФ с прореживанием (децимацией, происхождение слова “децимация”) (decimation in time, DIT).
В зависимости от способов деления отсчётов на части, вычисление ДПФ и объединение результатов различают БПФ с прореживанием по времени и по частоте.
Рассмотрим пример с делением отсчётов пополам.
Итак, пусть N – чётное число. ДПФ – это по определению:
.
Выделим в последней формуле два слагаемых, соответствующих элементам исходной последовательности с чётными и нечётными номерами:
.
Введём обозначения
и
,
а также вынесем из второй суммы общий
множитель
:
(v)
Две суммы в
последнем выражении представляют собой
ДПФ последовательностей
(отсчёты с чётными номерами) и
(отсчёты с нечётными номерами). Каждое
из этих ДПФ имеет размерность N/2.
Т.о.:
,
(*)
где
и
– ДПФ соответственно последовательностей
отсчётов с чётными и нечётными номерами.
,
.
Т.к. ДПФ размерности
N/2
даёт лишь N/2
спектральных отсчётов, непосредственно
использовать формулу (*) можно только
при
.
Для остальных n
(
)
следует
воспользоваться периодичностью спектра
дискретного сигнала (и, соответственно,
периодичность результатов ДПФ):
,
.
С учётом этого
формула (*) для
представляется в виде:
(vv)
Изобразим графически процесс вычисления 8-точечного ДПФ путём разбиения его на два 4-точечных ДПФ.
Прокомментируем блоки, выполняющие на рисунке объединения результатов двух 4х точечных ДПФ.
Каждый такой блок
имеет два входных и два выходных сигнала.
Один из выходных сигналов умножается
на комплексную экспоненту
,
после чего суммируемся со вторым
выходным сигналом и вычитается из него,
образуя таким образом два выходных
сигнала. Это соответствует реализации
формул (v)
и (vv).
Данная операция получила название
“бабочки ДПФ” (butterfly).
Она расшифровывается так.
Оценим количество
операций, необходимое для вычисления
ДПФ таким способом. Каждое из 2-х ДПФ
половиной размерности требует
операций. Кроме того, при вычислении
окончательных результатов каждый
спектральный коэффициент
(пока ещё не отсчёт) умножается на
комплексный экспоненциальный множитель.
Это добавляет ещё
операций. Итого получается
,
что почти вдвое меньше, чем при вычислении
ДПФ прямым способом.
Если
тоже является чётным числом (т.е. если
),
можно продолжить описанную процедуру,
выразив результат ДПФ через четыре ДПФ
размерности
.
Это позволяет ещё больше сократить
число требуемых вычислительных операций.
Наибольшая степень
ускорения вычислений достигается при
,
хотя в принципе, если N
не является простым числом, его можно
делить на любое число частей. Когда
,
деление последовательностей пополам
можно продолжать до тех пор, пока не
получатся две элементарные
последовательности, ДПФ которых
рассчитывается вообще без использования
операций умножения (достаточно вычислить
сумму и разность двух отсчётов). Число
требуемых при этом пар операций
“умножения-сложения” оценивается как
.
Т.о., вычислительные затраты по сравнению
с непосредственным вычислением ДПФ
уменьшаются в
.
Это число называется коэффициентом
ускорения вычислений (к.у.в.). При больших
N
это отношение весьма велико (например,
при
,
т.е. достигается более, чем 100-кратное
ускорение).
ДПФ с прореживанием по частоте.
(Англ. decimation in frequency, DIF).
Разделим исходную
последовательность
на две следующие друг за другом половины
(как ив предыдущем случае N
должно быть чётным числом):
.
Из второй суммы можно выделить множитель:
.
Это знакочередующая единица в зависимости от номера вычисляемого спектрального отсчёта n, поэтому дальше рассмотрим чётные и нечётные n по отдельности.
После выделения
множителя
комплексные экспоненты в обеих суммах
становится одинаковыми, поэтому выносим
их за скобки, объединяя две суммы:
,
.
Фигурирующие
здесь суммы – это ДПФ суммы и разности
половин исходной последовательности.
При этом разность перед вычислением
ДПФ умножается на комплексные экспоненты
.
Каждое из используемых ДПФ имеет
размерность
.
Итак, при прореживании по частоте
вычисления организуются следующим
образом:
1. Из исходной
последовательности
длиной N
получаем две последовательности
и
длиной
согласно следующим формулам:
,
.
2. ДПФ последовательности даёт спектральные отсчёты с чётными номерами, ДПФ последовательности – с нечётными:
,
.
Все уже сказано о возможности деления на число частей, отличное от двух, и о коэффициенте ускорения вычислений остаётся справедливым и для ДПФ с частотной децимацией.
Процесс вычисления 8-точечного ДПФ путём разбиения его на два 4-х точечных выглядит так.
Для получения обратного БПФ достаточно во всех формулах знак в показателях комплексных экспонент поменять и добавить на выходе или на входе деления на 2, а в общем случае – на используемый коэффициент прореживания.
Лекция №6.
Выводы.
1. Наибольший коэффициент ускорения вычислений по алгоритму БПФ достигается при длине исходной последовательности, равной целой степени двойки.
2. Если длина вектора исходной последовательности – простое число вычисления может быть выполнено только по прямой формуле ДПФ.
3. БПФ не является приближенным алгоритмом: при отсутствии вычислительных погрешностей он даст точно такой же результат, что и по формуле ДПФ. Ускорение достигается только за счёт оптимальной организации вычислений.
Взаимосвязь ДПФ и фильтрации (ДПФ как фильтр).
ДПФ как дискретная фильтрация.
ДПФ можно трактовать как обработку сигнала фильтром с соответствующей импульсной характеристикой (ИХ).
Эту ИХ можно
получить, если заметить, что
при целых n.
С учётом этого формула прямого ДПФ
запишется как:
.
Это есть дискретная
сверка: т.е. N-й
отсчёт результата обработки входного
сигнала
фильтром с ИХ такого вида:
Разумеется для каждого частотного отсчёта ИХ своя, чтобы это подчеркнуть в h мы использовали индекс n.
Найдем частотную
характеристику такого фильтра. Для
этого сначала получим функцию передачи
.
.
Воспользуемся формулой суммы конечной геометрической прогрессии:
.
.
Чтобы рассчитать
частотную характеристику используем
подстановку
.
.
После несложных тригонометрических преобразований можно получить:
.
где
.
Функция вида
называется функцией Дирихле:
,
где n – целое положительное число.
.
Так выглядит АЧХ
одного из каналов (3го)
ДПФ при
(частотная ось в номерах каналов).
Т.е. АЧХ фильтра,
реализуемого при вычислении ДПФ – это
гребёнка (набор) фильтров. Их АЧХ имеет
уровень первого бокового лепестка
примерно
.
При
это выражение стремится к
.
Алгоритм Герцеля (Hoertzel).
Функция передачи отдельного канала ДПФ соответствует рекурсивному фильтру N-го порядка:
.
Алгоритм вычисления отдельного отсчёта можно организовать только с помощью такого фильтра. Бывают случаи, когда нужно вычислить только некоторые отсчёты ДПФ, тогда такой алгоритм предпочтительнее, чем БПФ, принципиально ориентированное на получение срезу всех отсчётов БПФ.
Для практического
использования последнюю формулу
модифицируют. Прежде всего, заметим,
что слагаемое
в числителе выполняет лишь роль
“ограничителя”.
Оно “заставляет” ИХ закончится,
достигнув длины N
отсчётов. Такой фильтр вычисляет отсчёт
мгновенного спектра сигнала, всё время
используя последние N
отсчётов входного сигнала. Если же мы
будем обрабатывать сигнал порциями по
N
отсчётов, перед каждым поступлением
новой порции обнуляя фильтр, то конечность
длины ИХ не будет иметь значения. В этом
случае можно удалить слагаемое
из числителя
,
получив рекурсивный фильтр 1-го порядка:
.
Далее, коэффициент
обратной связи в рекурсивной ветви
данного фильтра является комплексным.
Поэтому при обработке каждого входного
отсчёта нужно выполнить четыре
вещественных умножения и столько же
сложений. Число операций можно сократить,
умножив числитель и знаменатель
на
и перейдя к фильтру второго порядка:
.
Рекурсивная часть
фильтра теперь имеет чисто вещественные
коэффициенты, к тому же коэффициент
при
равен единице. Поэтому при обработке
каждого отсчёта в рекурсивной ветви
фильтра надо выполнить всего одно
вещественное умножение и два вещественных
сложения.
В
рекурсивной ветви фильтра комплексные
коэффициенты исчезли, но появились в
числители (не рекурсивной ветви). Однако
если организовать фильтр в канонической
форме (пока мы не говорили, что это
такое), то расчёты будут выполняться
сначала в рекурсивной ветви, а потом в
не рекурсивной. И тогда комплексное
умножение и сложение в числителе будет
выполняться только один раз – для
получения окончательного результата.
Полученная структура фильтра показана на рисунке. Данный способ вычисления отдельного спектрального отсчёта называется алгоритмом Герцеля (Hoertzel).
Для вычисления
одного спектрального отсчёта по
алгоритму Герцеля требуется
вещественных умножений и
вещественных сложений. Если число
вычисляемых спектральных отсчётов
невелико, то этот алгоритм оказывается
эффективнее, чем расчёт всех спектральных
отсчётов с помощью БПФ.
Дискретная фильтрация с помощью ДПФ.
Мы представим ДПФ как фильтр. Займёмся обратным процессом - реализацией фильтрации с помощью ДПФ.
Рассматривая ДПФ, мы видим, что перемножение дискретных спектров соответствует круговой свёртке периодических последовательностей. В то же время сигнал на выходе фильтра – это линейная свёртка последовательностей конечной длины – входного сигнала и фильтра. Однако и с помощью ДПФ можно получить линейную дискретную свёртку. Для этого необходимо добавить к каждой из сворачиваемых последовательностей нулевые отсчёты в количестве, не меньше, чем длина второй последовательности минус единица. После этого длины последовательностей должны уровняться и стать равными как минимум длине линейной свёртки исходных последовательностей.
Круговая свёртка таких дополненных нулями последовательностей совпадает с линейной свёрткой исходных последовательностей.
Пример. Сворачиваемые
последовательности размера
.
,
.
Исходные сигналы.
1 |
2 |
4 |
8 |
2 |
3 |
4 |
5 |
Линейная свёртка:
|
1 |
2 |
4 |
8 |
|||
5 |
4 |
3 |
2 |
|
|||
|
1 |
2 |
4 |
8 |
|
5 |
4 |
3 |
2 |
|
|
|
1 |
2 |
4 |
8 |
5 |
4 |
3 |
2 |
|
1 |
2 |
4 |
8 |
5 |
4 |
3 |
2 |
1 |
2 |
4 |
8 |
|
|
5 |
4 |
3 |
2 |
1 |
2 |
4 |
8 |
|
|
|
5 |
4 |
3 |
2 |
|
1 |
2 |
4 |
8 |
|
||
|
5 |
4 |
3 |
2 |
||
2 |
7 |
18 |
41 |
50 |
52 |
40 |
Круговая свёртка.
1 |
2 |
4 |
8 |
2 |
5 |
4 |
3 |
1 |
2 |
4 |
8 |
3 |
2 |
5 |
4 |
1 |
2 |
4 |
8 |
4 |
3 |
2 |
5 |
1 |
2 |
4 |
8 |
5 |
4 |
3 |
2 |
52 |
59 |
58 |
41 |
Видно, что результаты ничем не похожи, и собственно, так и должно быть.
Дополним каждую последовательность тремя нулями, вычислим круговую свёртку ещё раз.
Исходные сигналы.
1 |
2 |
4 |
8 |
0 |
0 |
0 |
5 |
4 |
3 |
2 |
0 |
0 |
0 |
Вычисление круговой свёртки.
1 |
2 |
4 |
8 |
0 |
0 |
0 |
2 |
0 |
0 |
0 |
5 |
4 |
3 |
1 |
2 |
4 |
8 |
0 |
0 |
0 |
3 |
2 |
0 |
0 |
0 |
5 |
4 |
1 |
2 |
4 |
8 |
0 |
0 |
0 |
4 |
3 |
2 |
0 |
0 |
0 |
5 |
1 |
2 |
4 |
8 |
0 |
0 |
0 |
5 |
4 |
3 |
2 |
0 |
0 |
0 |
1 |
2 |
4 |
8 |
0 |
0 |
0 |
0 |
5 |
4 |
3 |
2 |
0 |
0 |
1 |
2 |
4 |
8 |
0 |
0 |
0 |
0 |
0 |
5 |
4 |
3 |
2 |
0 |
1 |
2 |
4 |
8 |
0 |
0 |
0 |
0 |
0 |
0 |
5 |
4 |
3 |
2 |
2 |
7 |
18 |
41 |
50 |
52 |
40 |
Результат совпадает с полученным ранее.
Поскольку ДПФ линейной свёртки есть произведение ДПФ сворачиваемых последовательностей, мы получаем такой алгоритм фильтрации в частотной области, который называется быстрая свёртка.
1. Последовательность отсчётов входного сигнала и импульсная характеристика фильтра (две сворачиваемые последовательности) дополняются нулями так, чтобы длины последовательностей стали равными и не меньше, чем сумма длин исходных последовательностей минус единица.
2. Вычисляются ДПФ дополненных нулями последовательностей.
3. Вычисленные ДПФ поэлементно перемножаются.
4. Вычисляется обратное ДПФ от результата перемножения.
Количество нулей, добавляемых к последовательностям, может быть больше минимально необходимого – в этом случае лишь появится соответствующее количество нулей (крайних) в вычисленной свёртке. Т.о. мы всегда можем выбрать количество нулей таким, чтобы длина получившихся последовательностей была равна степени двойки, и использовать БПФ. Именно в этом случае алгоритм, изложенный ранее, называется быстрый свёрткой. Он даёт существенную вычислительную экономию.
Секционированные свёртки.
Недостатком описанного метода фильтрации является то, что входной сигнал обрабатывается сразу целиком, т.е. отсчёты можно начать вычислять только после окончания поступления входного сигнала. На практике длительность входного сигнала, как правило, значительно больше, чем длина ИХ. Например, если мы фильтруем речевой сигнал, он поступает, пока не закончится произнесение сообщения. Необходимо обрабатывать сигналы каким-то образом по мере их поступления.
Возможным решением проблемы является блочная или секционированная свёртка. При этом поступающий входной сигнал делится на секции или блоки заданного размера. Эти блоки подвергаются фильтрации в частотной области по отдельности. При выборе размера блока необходимо реализовать компромисс между числом операций (чем больше длина секции, тем эффективнее использование БПФ) и задержкой в получении выходного сигнала (чем меньше размер секции, тем меньше задержка).
С идеей объединения результатов секционирования свёртки разберёмся на примере.
Пусть
;
.
Сначала вычислим линейную свёртку напрямую.
.
Разобьем входной сигнал на две секции по 6 отсчётов и обработаем по отдельности:
.
Мы видим, что
первые шесть отсчётов
и последние 6 отсчётов
верные. Оставшиеся два отсчёта у могут
быть получены (7й
и 8й)
суммированием двух последних элементов
и двух первых элементов
.
Итак, секционная фильтрация осуществляется так:
1. Входной сигнал разбивается на блоки длиной N отсчётов.
2. Каждый блок
фильтруется независимо. Длина выходного
сигнала составляет (
)
отсчётов, где M
– длина ИХ фильтра. Для повышения
эффективности обычно фильтрация
осуществляется в частотной области с
применением БПФ.
3. Блоки или секции
объединяются, при этом их крайние (
)
отсчётов перекрываются и суммируются.
Данный алгоритм секционированной свёртки называется перекрытие с суммированием (overlap-add). Другой возможный вариант – перекрытие с накоплением (overlap-save). Входной сигнал разбивается на блоки, перекрывающиеся по краям на ( ) отсчёт. У выходного сигнала для каждого блока отбрасываются крайние хвосты длиной по ( ) отсчётов с каждого края. Затем секции объединяются без перекрытия.
Лекция №7