Алгоритмы вычисления ДПФ. Быстрое преобразование Фурье.

Цифровая обработка сигналов; лекция 17 апреля 2017 г. МФТИ 10. Алгоритмы вычисления ДПФ Дискретное во времени преобразование Фурье (ДВПФ) и ДПФ обеспечивают успешный анализ и разработку систем с дискретным временем. В п. 8 мы выяснили, что N коэффициентов ДПФ равны отмасштабированным отсчётам X (n), 0  n  N  1, последовательности x(k ), 0  k  N  1 , X (n) / N ДВПФ X ( ) , взятым с шагом   1/ N . В этом параграфе обсуждаются конкретные методы вычисления ДПФ. Особое внимание будет уделено высокоэффективным алгоритмам цифрового вычисления, которые обычно называют быстрым преобразованием Фурье (БПФ). Алгоритмы БПФ наиболее хорошо работают, когда нужно определить все N коэффициентов ДПФ. Однако практически часто достаточно знать только часть из N отсчётов ДПФ. В этом случае более эффективными могут оказаться другие алгоритмы. Примером является алгоритм Герцеля. Будут также рассмотрены алгоритмы скользящего ДПФ-анализа. 10.1. Быстрое преобразование Фурье Быстрое преобразование Фурье (БПФ) представляет собой эффективный метод вычисления ДПФ. Его эффективность заключается в существенном уменьшении числа операций умножения и суммирования, затрачиваемых для получения всех N коэффициентов ДПФ, которые запишем в виде Equation Section 10 X ( n)  N 1  x (k ) W k 0 где WNn k  exp( j nk N , 2 n k ) – дискретные экспоненциальные функции, N (10.1) n  0, 1, 2, , N 1 . Здесь и далее масштабирующий множитель 1 / N в прямом преобразовании для простоты опущен. Прямое вычисление всех N коэффициентов ДПФ по (10.1) требует N 2 операций типа «комплексное умножение плюс сложение». Основная идея БПФ состоит в том, чтобы разбить исходную N-точечную последовательность на две более короткие последовательности, из ДПФ которых можно получить ДПФ исходной N-точечной последовательности. Так, например, если N чётное, то исходная N-точечная последовательность разбивается на две ( N / 2)-точечные последовательности. Для вычисления искомого Nточечного ДПФ потребуется ( N / 2)2  2 операций типа «комплексное умножение плюс сложение», т. е. вдвое меньше по сравнению с прямым вычислением. Такое уменьшение размерности ДПФ вдвое называется итерацией. Эту операцию можно повторить, вычисляя вместо ( N / 2)-точечного ДПФ два ( N / 4)-точечных ДПФ, сокращая тем самым объём вычислений ещё в два раза. Процесс уменьшения размера ДПФ продолжается до тех пор, пока не останутся только двухточечные ДПФ. Как будет показано далее, алгоритм БПФ с основанием 2 затрачивает на вычисление искомого N-точечного ДПФ N log 2 N операций, так что выигрыш в числе операций составляет N / log 2 N и при больших N может достигать нескольких сотен. Алгоритм БПФ основан на периодичности ядра преобразования WNn k . Алгоритм БПФ с составным основанием Прямой метод вычисления всех коэффициентов в соответствии с (10.1) требует N 2 операций комплексного умножения и суммирования. Это число может оказаться очень большим. Возможность сокращения числа операций основывается на представлении одномерного ДПФ в виде многомерного. Для этого необходимо, чтобы длина массива являлась составным числом N  N1  N2  ...  N p . (10.2) Разберем эту возможность на пример двух сомножителей: N  N1 N 2 . 1 (10.3) Цифровая обработка сигналов; лекция 17 апреля 2017 г. МФТИ При этом входной массив из N отсчетов разбивается на N 2 блоков Рис. 10.1. Образование двумерного массива из одномерного по N1 элементов в каждом (рис. 10.1). Расположив блоки один под другим, получаем двумерный массив. Нетрудно убедиться, что одномерный номер k может быть представлен в виде k  N1 k2  k1 . (10.4) Первое слагаемое соответствует целому числу  k 2  блоков, предшествующих номеру k , а второе слагаемое определяет номер элемента в блоке, содержащем номер k . ДПФ этого двумерного массива также будет иметь вид двумерного массива с переменными n1 и n2 : n1  0, 1, 2, ..., N1  1; n2  0, 1, 2, ..., N 2  1. Пусть двумерный массив имеет вид (рис. 10.2) Рис. 10.2 Одномерный номер n может быть представлен в виде n  N 2 n 1  n2 . (10.5) Для базисной функции ДПФ WNnk с учетом (10.4) и (10.5) можем записать WNk n  WN N1 k2  k1   N 2 n1  n2   WNN k2 n1 WNN1 k2 n 2 WNN2 k1 n1 WNk1 n 2 . Учитывая, что k n 2 WNN k2 n1  1; WNN1k2 n 2  WN22 k n 1 ; WNN2 k1 n1  WN11 получаем k n 1 WNk n  WN11 k n 2 WNk1 n 2 WN22 Подставим это выражение в (10.1), тогда 2 . , Цифровая обработка сигналов; лекция 17 апреля 2017 г. МФТИ X  n 1, n 2   N1 1  W k n 1 1 k 1 0 N1 WNk 1n 2 N2 1  x k k 2 0 k n 2 2 , k 1  WN22 (10.6) . N2-точечные ДПФ вектор поворота N1-точечные ДПФ В соответствии с (10.6) можно выделить следующие этапы вычисления коэффициентов ДПФ: k n 1) сначала вычисляются N2-точечные ДПФ с ядром WN22 2 по столбцам матрицы  x  k2 , k1  : Y  k1 , n2   N2 1  x  k , k  WNk n 2 2 k 2 0 1 2 2 ; число таких ДПФ соответствует числу столбцов матрицы  x  k2 , k1  и равно N1 ; 2) n k 1 умножение на поворачивающие множители WN 2 , в результате образуется новый массив Z  k1 , n2   Y  k1 , n2  WNn 2k 1 ; число умножений равно N1  N 2 ; 3) n k 1 вычисляются ДПФ с ядром WN11 по столбцам матрицы  Z  k1 , n2  : X  n 1, n 2   N1 1  Z k k 1 0 , n 2  WN 12 1 . n k 1 Полное число операций в сумме по всем трем этапам будет M  N1   N2   N1  N2  N2   N1   N1  N2   N1  N2  1 , 2 2 т. е. меньше, чем N 2   N1    N2  при прямом методе вычислений. 2 2 Таким образом, исходное N-точечное ДПФ оказалось сведённым к ДПФ, производимым над уменьшенными массивами N 2 и N1 . Пример. 10.1. Пусть N  35  5  7. БПФ-алгоритм выполняется за M  5  7   5  7  1  455 операций вместо N 2  1225 при прямом вычислении по (10.1). Если N1 и N 2 сами являются составными, то каждое из ДПФ N 1 и ДПФ N 2 может быть выполнено с применением рассмотренного алгоритма. Выигрыш в числе операций при этом еще более возрастает. В общем случае, если выполняется (10.2), такая процедура позволяет уменьшить общее число операций при вычислении ДПФ до величины p N N i 1 p i вместо N 2  N  N i . i 1 В заключение на рис. 10.3 приведен граф БПФ для N  6  2  3 ( N1  2, N 2  3 ). Рис. 10.3. Граф БПФ для составного N  6 3 Цифровая обработка сигналов; лекция 17 апреля 2017 г. МФТИ Умножение на поворачивающие множители обозначено стрелками, рядом с которыми записаны значения n k коэффициентов W6 2 1 . Звёздочкой обозначены ячейки памяти. Алгоритмы БПФ с основанием 2 Рассмотрим случай, когда N является степенью двойки. Пусть в (10.2) N1  2, N 2  N / 2, тогда на первой итерации N-точечное ДПФ представляется через двухточечные и N / 2-точечные ДПФ: X  n 1, n 2   1  W n k k 0 1 2 1 WNn 2k 1 1 ( N / 2) 1  k 2 0 WNn /22k 2 x  k 2 , k 1  . (10.7) Эта же схема вычислений может быть использована на следующей итерации для получения каждого из N / 2-точечных ДПФ. В результате перейдем к N / 4-точечным ДПФ и т. д., пока не останутся только двухточечные ДПФ. Соответствующий граф вычислений для N  8 приведен на рис. 10.4. Рис. 10.4. Граф алгоритма БПФ с прореживанием по времени для N  8 В данном графе умножение на вектор поворота выполняется до двухточечного ДПФ. Такой алгоритм получил название алгоритма БПФ с прореживанием по времени. Его базовая операция «бабочка» приведена на рис. 10.5. Само двухточечное ДПФ не содержит умножений: все умножения сводятся к умножениям на поворачивающие множители WNm , расположенные у стрелочек. Рис. 10.5. Базовая операция «бабочка» алгоритма БПФ с прореживанием по времени На каждой итерации имеется N / 2 базовых операций алгоритма БПФ, каждая из которых состоит из умножения на вектор поворота и двухточечного ДПФ. Еще одной особенностью графа на рис. 10.4 является то, что входные отсчеты расположены в так называемом разрядно-инверсном порядке, а выходные – в естественном. 4 Цифровая обработка сигналов; лекция 17 апреля 2017 г. МФТИ Номер Двоичное представление Инверсия разрядов Разрядно-инв. порядок 000 000 1 001 100 4 2 010 010 2 3 011 110 6 4 100 001 1 5 101 101 5 6 110 011 3 7 111 111 7 Пусть теперь в (10.2) N1  N / 2, X  n 1, n 2   N 2  2. Тогда N / 2 1  WNn /2k k 0 1 1 WNn 2k 1 1 1  x  k , k W2n k k 0 2 2 1 2 . (10.8) 2 повор. двухточечное. ДПФ множит. N / 2-точечные ДПФ Здесь вначале выполняются двухточечные ДПФ над входным массивом, затем N / 2-точечные. На следующей итерации аналогично представляются N / 2-точечные ДПФ. Соответствующий данной процедуре граф вычислений для N  8 представлен на рис. 10.6. Рис. 10.6. Граф алгоритма БПФ с прореживанием по частоте для N  8 Особенностью этого графа является естественный порядок данных на входе и разрядно-инверсный порядок на выходе. В этом алгоритме в отличие от предыдущего самые сложные умножения на поворачивающие множители производятся в начале графа. 5 Цифровая обработка сигналов; лекция 17 апреля 2017 г. МФТИ В базовой операции алгоритма (рис.10.7) сначала выполняется двухточечное ДПФ, а затем умножение на вектор поворота. Такой алгоритм носит название алгоритма БПФ с прореживанием по частоте. Рис. 10.7. Базовая операция «бабочка» алгоритма БПФ с прореживанием по частоте Граф алгоритма БПФ с прореживанием по частоте можно получить из рис. 10.4, если входные и выходные отсчёты поменять местами, направления стрелок изменить на обратное и читать справа налево. В алгоритме БПФ по основанию 2 на каждой итерации, по крайней мере, половина поворачивающих множителей – единицы, а число итераций равно log 2 N , поэтому общее число комплексных умножений не превышает M   N 2  log 2 N . На каждой итерации в алгоритме БПФ по основанию 2 выполняется N / 2 двухточечных ДПФ или N сложений. Всего в алгоритме БПФ по основанию 2 затрачивается A  N log 2 N комплексных сложений. Можно считать, что при вычислении всех N коэффициентов ДПФ требуется около N log 2 N вычислительных операций типа «комплексное умножение плюс сложение» вместо N 2 при прямом вычислении по (10.1). При больших N выигрыш в числе операций может быть значительным. Например, при N  210  1024 этот выигрыш составляет N /log 2 N  100. Вычисление обратного ДПФ может быть проведено с использованием любого из описанных выше алгоритмов БПФ: x k     N 1  X  n  WN nk    X   n  WNnk  .  n0 n  0  N 1 (10.9) Выражение в квадратных скобках представляет собой ДПФ последовательности X   n  , комплексносопряженной с X  n  , и может быть вычислено с использованием одного из описанных выше алгоритмов. Искомая последовательность x  k  получается, таким образом, комплексным сопряжением этого ДПФ. Если последовательность x  k  действительная, то комплексного сопряжения ДПФ не требуется. В (10.9) попрежнему WNnk  e  j 2 nk / N . При вычислениях на каждой базовой операции БПФ по основанию 2, включающей вычисление одного двухточечного ДПФ и умножение на поворачивающий множитель, выходные результаты можно размещать в те же ячейки памяти, где находились исходные данные. Такой алгоритм БПФ называется алгоритмом с замещением. Графы, представленные на рис. 10.4 и рис. 10.6, соответствуют алгоритмам с замещением. Математическая операция перехода из одномерного массива данных к двумерному является основой всех алгоритмов БПФ. Существует большое разнообразие алгоритмов БПФ, отличающихся порядком следования входных, промежуточных и выходных отсчетов, организацией вычислений, регулярностью структуры и т. д. Известна, например, другая, часто используемая модификация алгоритмов БПФ, отличная от алгоритмов с замещением. Она получила название алгоритмов с постоянной структурой. Граф этого алгоритма для N  8 с нормальным порядком входных и разрядно-инверсным выходных отсчетов приведен на рис. 10.8. Отличительной особенностью этого графа является то, что все итерации имеют одинаковую структуру, различаются лишь поворачивающие множители между итерациями. 6 Цифровая обработка сигналов; лекция 17 апреля 2017 г. МФТИ Рис. 10.8. Граф алгоритма БПФ с постоянной структурой с нормальным порядком входных и разрядно- инверсным порядком выходных отсчетов для N  8 Этот алгоритм имеет и другую модификацию, при которой порядок входных отсчётов разрядноинверсный, а выходных – нормальный (рис. 10.9). Рис.10.9. Граф алгоритма БПФ с постоянной структурой с разрядно-инверсным порядком входных и нормальным выходных отсчётов для N  8 Алгоритм БПФ по основанию 2 был предложен Кули и Тьюки в 1965 году и дал огромный импульс развитию цифровых методов обработки сигналов. Однако алгоритмы БПФ в наиболее общем виде (10.6) были получены известным математиком Гауссом (1777–1855) и опубликованы в 1865 году. Примеры решения задач Пример 4. Найти последовательность x(k) по ее z- преобразованию   2 1  z 1 1  2 z 1  z 2 X  z    z 1 3 1 1 2  1 1  1 1 z  z 1  z  1  z 2 2  2  Для нахождения последовательности x(k) применим метод простейших дробей. Оба полюса имеют кратность 1. Поэтому X(z) можно представить в виде A1 A2 X  z   B0   1 1  0,5 z 1  z 1 Константу B0 ищем делением в столбик  z 2  2 z 1  1  0,5 z 2  1,5 z 1  1 2 z 2  3 z 1  2 5 z 1  1 Степень числителя и знаменателя совпадают, поэтому потребовался лишь один шаг деления при вычислении остатка. X(z) можно переписать как  X  z  2  1  5 z 1 1  0,5 z 1 1  z 1   7  Цифровая обработка сигналов; лекция 17 апреля 2017 г. МФТИ Теперь находим коэффициенты А1 и А2  1  5 z 1   A1  2   1  0,5 z 1 1  z 1     2  1  5 z 1 A2   2   1  0,5 z 1 1  z 1       1  1  z 1    9   2    z 1    1  z 1       8   z 1 Следовательно, X  z  2  9 8   1 1  1  z 1 1  z   2    Поскольку область сходимости определена неравенством |z| >1, находим z 2   2 1 k  1 1  k k  1 2 1  z 1 2 1  k  1  z 1 Благодаря линейности z- преобразования получаем 9 x  k   2 1 k   k u  k   8u  k  . 2 Пример 5. Передаточная функция фильтра a0 . H z  1 1  b1 z  b2 z 2 Найти соответствующее разностное уравнение. Изобразить блок-схему фильтра. Решение. y  k   b1 y  k  1  b2 y  k  2   a0 x  k  x(k) a0 y(k) ∑ z−1 ∑ y(k−1) b1 z−1 y(k−2) b2 8 Цифровая обработка сигналов; лекция 17 апреля 2017 г. МФТИ Пример 6. Разностное уравнение фильтра 1 1 y  k  2    y  k  1  y  k   2 x  k  , 6 6 Найти отклик на сигнал x(k) = σ(k). По теореме опережающего сдвига y  0   0, y 1  1 . y  k  1  z Y  z   y  0    z  Y  z  ; y  k  2   z  zY  z   y 1   z  zY  z   1 . Исходное уравнение в терминах z- преобразования преобразуется к виду 1 1 z  zY  z   1  zY  z   Y  z   2 X  z  . 6 6 Так как 1 z X  z   , 1 1 z  z  1 (а) то уравнение (а) будет 1 1 2z z 2Y  z   z  zY  z   Y  z   . 6 6  z  1 Поэтому z2  z z2  z   1 1  1  2 1   z  1  z  z    z  1  z   z   6 6 2  3   . z 1  z 2  . z 1 z 2  1  1  z 1    6 6   Y  z    Находим 1 1 k1  3  вычет в точке z=1; k2  1,8  вычет в точке z= ; k3  0, 2  вычет в точке z=  . 2 3 Разложим Y  z  на простые дроби Y  z  3 1,8 0, 2   1 1 1 1 z 1  z 1 1  z 1 2 3 С учетом того, что ak  1 1  az 1 получаем k к 1  1 y  k   3  1,8    0, 2    . 2  3 9 Цифровая обработка сигналов; лекция 17 апреля 2017 г. МФТИ Задачи к лекции 17 апреля 1. КИХ-фильтр задаётся отсчётами импульсной характеристики h(0)  h(4)  2 h(1)  h(3)  1 h(2)  0,5 h( k )  0 для других k Для рекурсивного и нерекурсивного способов реализации :  записать разностные уравнения,  изобразить блок-схемы фильтра,  определить передаточные функции,  найти фазовую задержку. 2. Сигнал x  t   cos w0t дискретизуется с шагом Δt и из полученной последовательности вырезается отрезок x(k) длиной в N отсчетов, т.е. k = 0, 1, 2…N-1. Для проведения спектрального анализа вычисляется N-точечное ДПФ N 1 X  n   x  k  e j 2 nk N , п  0, 1, 2 N –1. k 0 а) Для фиксированных w0, N, n0 подобрать Δt так, чтобы X  n0   0, X  N  n0   0, а остальные Х(п) равны нулю. б) Единственное ли это Δt? 3. Рассмотрим последовательность x  k  с периодом N=10 ••• ••• 0 1 2 3 4 5 6 7 8 9 10 k а) Найти и изобразить по модулю ДВПФ этой периодической последовательности. б) Изобразить модули и фазы коэффициентов ДПФ. 10