#Lectures #DSP

Быстрое преобразование Фурье — конспект лекции

25 Mar

by Vladimir Leonidov 1

Введение

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

Вспомним выражение для ДПФ:

(1) $\begin{equation*} X[m] = \sum\limits_{n=0}^{N-1}x[n]\cdot e^{-j2\pi n m/N} \end{equation*}$

Из него видно, что для расчёта $N$ -точечного ДПФ требуется $N^2$ вычислений с комплексными числами. Получается ресурсозатратно, особенно если мы имеем дело с большим количеством отсчётов.

Учёным не нравилось делать много вычислений, и несколько человек (существуют разные мнения, кто был первым) предложили алгоритмы, позволяющие вычислять ДПФ, производя меньшее количество математических операций. Называются такие алгоритмы Быстрое преобразование Фурье (БПФ).

Хочу сразу отметить, что БПФ возвращает абсолютно такие же результаты, что и ДПФ, использует в основе ту же формулу, только позволяет достичь этих результатов гораздо быстрее.

В рамках данной лекции рассмотрим два самых распространённых алгоритма БПФ — с прореживанием по времени и с прореживанием по частоте.

Итак, выше мы уже отметили, что для $N$ -точечного ДПФ требуется $N^2$ операций комплексного умножения и сложения. А что, если разделить исходный сигнал на два длиной $N/2$ отсчётов каждый? Тогда получается, что для вычисления ДПФ от одной из половин потребуется $(N/2)^2 = N^2/4$ математических операций. Или $N^2/2$ для двух половин (если не считать процедуру разделения и объединения). Получается, эффективность вычислений возрастает в два раза! Но это ещё не предел. Мы можем делить наш сигнал до тех пор, пока не получим набор из сигналов, состоящих из двух отсчётов:

Теперь давайте немного модифицируем выражение (1). Пусть $W_N^{n m} = e^{-j2\pi n m/N}$ , тогда:

(2) $\begin{equation*} X[m] = \sum\limits_{n=0}^{N-1}x[n]\cdot W_N^{n m} \end{equation*}$

На примере сигнала из $N=8$ отсчётов рассчитаем несколько значений $W_N^{n m}$ :

$\begin{align*} W_8^0 &= e^{0} = 1\\ W_8^1 &= e^{-j2\pi/8} = e^{-j\pi/4} = \frac{1-j}{\sqrt{2}} = a\\ W_8^2 &= a^2 = \left( \frac{1-j}{\sqrt{2}} \right)^2 = \frac{1-2j+j^2}{2} = -j\\ W_8^3 &= a^3 = a^2 \cdot a = -a^*\\ W_8^4 &= a^4 = (a^2)^2 = -1 \\ W_8^5 &= a^5 = a^4\cdot a = -a\\ W_8^6 &= a^6 = a^4\cdot a^2 = j\\ W_8^7 &= a^7 = a^4\cdot a^3 = a^*\\ W_8^8 &= a^8 = a^4 \cdot a^4 = 1 \end{align*}$

Можно заметить, что первая половина значений $W_N^{n m}$ отличается от второй половины только знаком:

$\begin{align*} W_8^0 &= -W_8^4\\ W_8^1 &= -W_8^5\\ W_8^2 &= -W_8^6\\ W_8^3 &= -W_8^7 \end{align*}$

Если произведение $n m$ выходит за диапазон от 0 до 7, мы также увидим повторяемость этих значений. Например, при $n=6$ , а $m=7$ :

$\begin{align*} W_8^{42} = a^{42} = (a^8)^5 \cdot a^2 = a^2 = -j \end{align*}$

Следовательно, в процессе ДПФ вычисление одних и тех же значений повторяется несколько раз. Наша задача — оптимизировать процесс так, чтобы избежать выполнение повторяющихся операций.

Алгоритм БПФ с прореживанием по времени

Давайте разделим сигнал на отсчёты с чётными и нечётными номерами, как показано на рисунке:

Эта операция называется прореживание по времени. Математически это можно записать следующим образом:

(3) $\begin{equation*} X[m] = \sum\limits_{n=0}^{N/2-1}x[2n]\cdot W_N^{2nm} + \sum\limits_{n=0}^{N/2-1}x[2n+1]\cdot W_N^{(2n+1)m} \end{equation*}$

Отметим, что:

(4) $\begin{equation*} W_N^{2n m} = e^{-j2\pi 2n m/N} = e^{\frac{-j2\pi n m}{N/2}} = W_{\frac{N}{2}}^{nm} \end{equation*}$

(5) $\begin{equation*} W_N^{(2n+1)m} = W_N^{2nm} \cdot W_N^m = W_{\frac{N}{2}}^{nm} \cdot W_N^m \end{equation*}$

Тогда, выражение (3) принимает вид:

(6) $\begin{equation*} X[m] = \sum\limits_{n=0}^{N/2-1}x[2n]\cdot W_{\frac{N}{2}}^{nm} + W_N^m\sum\limits_{n=0}^{N/2-1}x[2n+1]\cdot W_{\frac{N}{2}}^{nm} \end{equation*}$

Это был расчёт ДПФ для первых $N/2$ отсчётов. Теперь запишем подобное выражение для второй половины отсчётов:

(7) $\begin{equation*} \begin{aligned} X\left[m+\frac{N}{2}\right] = \sum\limits_{n=0}^{N/2-1}x[2n]\cdot W_N^{2n(m+\frac{N}{2})} +\\ +\sum\limits_{n=0}^{N/2-1}x[2n+1]\cdot W_N^{(2n+1)(m+\frac{N}{2})} \end{aligned} \end{equation*}$

Давайте теперь по очереди разбираться со всеми множителями:

(8) $\begin{equation*} \begin{aligned} W_N^{2n(m+\frac{N}{2})} = W_N^{2nm} \cdot W_N^{2n\frac{N}{2}} = W_{\frac{N}{2}}^{n m} \cdot e^{\frac{-j2\pi 2n\frac{N}{2}}{N}} = \\ = W_{\frac{N}{2}}^{n m} \cdot e^{-j2\pi n} \end{aligned} \end{equation*}$

Т.к. $e^{-j2\pi n} = 1$ , то выражение (8) принимает вид:

(9) $\begin{equation*} W_N^{2n(m+\frac{N}{2})} = W_{\frac{N}{2}}^{n m} \end{equation*}$

Рассмотрим следующий множитель:

(10) $\begin{equation*} W_N^{(2n+1)(m+\frac{N}{2})} = W_N^{2n m} \cdot W_N^{2nN/2} \cdot W_N^m \cdot W_N^{N/2} \end{equation*}$

Т.к.

(11) $\begin{equation*} W_N^{N/2} = e^{\frac{-j2\pi N}{2N}} = e^{-j\pi} = -1, \end{equation*}$

то выражение (10) можно упростить:

(12) $\begin{equation*} W_N^{(2n+1)(m+\frac{N}{2})} = -W_{\frac{N}{2}}^{n m} \cdot W_N^m \end{equation*}$

Подставим результаты из (9) и (12) в уравнение (7) и получим:

(13) $\begin{equation*} \begin{aligned} X\left[m+\frac{N}{2}\right] = \sum\limits_{n=0}^{N/2-1}x[2n]\cdot W_{\frac{N}{2}}^{nm} - \\ - W_N^m\sum\limits_{n=0}^{N/2-1}x[2n+1]\cdot W_{\frac{N}{2}}^{n m} \end{aligned} \end{equation*}$

Для простоты восприятия сделаем две замены:

(14) $\begin{equation*} G[m] = \sum\limits_{n=0}^{N/2-1}x[2n]\cdot W_{\frac{N}{2}}^{n m} \end{equation*}$

(15) $\begin{equation*} H[m] = \sum\limits_{n=0}^{N/2-1}x[2n+1]\cdot W_{\frac{N}{2}}^{n m} \end{equation*}$

Тогда, с учётом (14) и (15), сгруппируем уравнения (6) и (13):

(16) $\begin{equation*} \begin{aligned} X[m] = G[m]+W_N^m \cdot H[m]\\ X\left[m+\frac{N}{2}\right] = G[m]-W_N^m \cdot H[m] \end{aligned} \end{equation*}$

Полученное выражение (16) называется графом “бабочка”. Почему? Это станет ясно, если посмотреть на его графическое представление:

Данное выражение лежит в основе БПФ с прореживанием по времени.

Таким образом, берём анализируемый сигнал и продолжаем разбивать его на сигналы, состоящие из чётных и нечётных индексов до тех пор, пока не получим набор сигналов из двух отсчётов. На рисунке ниже показан пример такого разбиения для сигнала из восьми отсчётов.

Алгоритм БПФ с прореживанием по времени для N=8

Для вычисления ДПФ такого сигнала требуется 3 стадии разбиения. На первом этапе получаем набор из четырёх ДПФ от сигналов из двух отсчётов:

$\begin{align*} X_{00}[0] = x[0] + W_2^0 x[4]\\ X_{00}[1] = x[0] - W_2^0 x[4]\\\\ X_{01}[0] = x[2] + W_2^0 x[6]\\ X_{01}[1] = x[2] - W_2^0 x[6]\\\\ X_{02}[0] = x[1] + W_2^0 x[5]\\ X_{02}[1] = x[1] - W_2^0 x[5]\\\\ X_{03}[0] = x[3] + W_2^0 x[7]\\ X_{03}[1] = x[3] - W_2^0 x[7] \end{align*}$

Далее, на основе полученных результатов формируется два ДПФ для сигналов, состоящих из четырёх отсчётов:

$\begin{align*} X_{10}[0] = X_{00}[0] + W_4^0 X_{01}[0]\\ X_{10}[1] = X_{00}[1] + W_4^1 X_{01}[1]\\ X_{10}[2] = X_{00}[0] - W_4^0 X_{01}[0]\\ X_{10}[3] = X_{00}[1] - W_4^1 X_{01}[1]\\\\ X_{11}[0] = X_{02}[0] + W_4^0 X_{03}[0]\\ X_{11}[1] = X_{02}[1] + W_4^1 X_{03}[1]\\ X_{11}[2] = X_{02}[0] - W_4^0 X_{03}[0]\\ X_{11}[3] = X_{02}[1] - W_4^1 X_{03}[1] \end{align*}$

И на последнем этапе формируются итоговые результаты восьмиточечного ДПФ:

$\begin{align*} X[0] = X_{10}[0] + W_8^0 X_{11}[0]\\ X[1] = X_{10}[1] + W_8^1 X_{11}[1]\\ X[2] = X_{10}[2] + W_8^2 X_{11}[2]\\ X[3] = X_{10}[3] + W_8^3 X_{11}[3]\\ X[4] = X_{10}[0] - W_8^0 X_{11}[0]\\ X[5] = X_{10}[1] - W_8^1 X_{11}[1]\\ X[6] = X_{10}[2] - W_8^2 X_{11}[2]\\ X[7] = X_{10}[3] - W_8^3 X_{11}[3] \end{align*}$

ДПФ рассчитано!

Двоично-инверсная перестановка

Если количество отсчётов исходного сигнала $N=2^n$ , где $n$ — целое положительное число, то его разбиение на сигналы, состоящие из чётных и нечётных индексов вплоть до последнего уровня можно сделать легко и за одну итерацию с помощью двоично-инверсной перестановки.

Для этого записываем индексы всех отсчётов в двоичной системе счисления, при этом должны быть записаны все $n$ бит индекса, включая ведущие нули. Затем зеркально отображаем код каждого из этих двоичных чисел и записываем полученные результаты обратно в десятичную систему счисления. Готово! Если расставить элементы исходного массива в соответствии с полученными индексами, мы получим подряд идущие пары отсчётов для двухточечного ДПФ. Пример, как это работает для сигнала из восьми отсчётов, представлен в таблице:

Номер до перестановки	Двоичное	Двоично-инверсная перестановка	Десятичное
0	000	000	0
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

Можете сравнить с рисунком выше, действительно совпадает. Именно поэтому, рассматриваемые нами алгоритмы также называются БПФ по основанию два. У них есть одно небольшое ограничение — количество отсчётов анализируемого сигнала должно быть равно степени двойки (например, 16, 32, 64, 128, 256 и т.д.). “Но как же быть, если, скажем, в моём сигнале всего 200 отсчётов?” — спросите вы. Ничего страшного, нужно просто дополнить сигнал нулевыми отсчётами, пока его длина не станет равна ближайшей степени двойки.

Алгоритм БПФ с прореживанием по частоте

Как вы догадались из названия раздела, прореживать можно не только время, но ещё и частоту. Граф “бабочка” для БПФ с прореживанием по частоте показан на рисунке:

Подробно рассматривать данный алгоритм не будем, просто примем его как данность. Пример БПФ с прореживанием по частоте для сигнала, состоящего из 8 отсчётов показан ниже:

Алгоритм БПФ с прореживанием по частоте для N=8

Прореживание (двоично-инверсная перестановка) в данном случае производится уже после вычисления ДПФ, т.е. по частотным отсчётам, отсюда и такое название.

По количеству операций умножения данный алгоритм схож с алгоритмом прореживания по времени, поэтому эффективность их схожая. Какой метод выбрать — решать вам.

Поворотные коэффициенты

В вышеупомянутых “бабочках” есть множитель $W_N^m$ . Он называется поворотным коэффициентом. На нижнем уровне БПФ, когда имеем дело с двумя отсчётами, требуется всего один поворотный коэффициент $W_2^0 = 1$ . Умножение на единицу считается тривиальной операцией, откуда следует, что на первом уровне не требуется ни одной операции умножения, только сложение и вычитание.

На втором уровне у нас получается два поворотных коэффициента $W_4^0=1$ и $W_4^1=-j$ . Умножение на $-j$ также делается просто: нужно поменять местами действительную и мнимую часть, а также инвертировать знак мнимой части. На этом уровне также не требуется операций умножения.

На следующем уровне требуется четыре поворотных коэффициента $W_8^0$ , $W_8^1$ , $W_8^2$ и $W_8^3$ , их значения мы рассчитывали в начале нашей статьи.

Давайте проиллюстрируем на графике все упомянутые выше поворотные коэффициенты (и даже немного больше):

Из графика видно, что на каждом следующем уровне количество поворотных коэффициентов увеличивается в два раза, при этом половина из них совпадает с поворотными коэффициентам предыдущего уровня.

Когда БПФ используется для решения каких-либо задач в устройствах цифровой обработки сигналов, количество отсчётов анализируемого сигнала как правило заранее известно и скорее всего не будет меняться в процессе работы устройства. Поэтому, чтобы не тратить время на вычисление $W_N^m$ , их заранее рассчитанные значения для заданного количества отсчётов и всех уровней разбиения записывают в памяти в виде таблицы констант, которые в дальнейшем будут использоваться для умножения в графах “бабочка”.

Выводы

БПФ — это всего лишь алгоритм эффективного вычисления ДПФ, поэтому результаты вычисления БПФ и ДПФ получаются абсолютно идентичны.
Для вычисления ДПФ сигнала, состоящего из $N$ отсчётов нам требуется $N^2$ вычислений с комплексными числами, в случае с БПФ — $\frac{N}{2}\log_2{N}$ вычислений с комплексными числам.

Чтобы наглядно продемонстрировать эффективность алгоритма, давайте взглянем на таблицу:

Кол-во отсчётов, N	Количество вычислений с комплексными числами		Эффективность
Кол-во отсчётов, N	ДПФ	БПФ	Эффективность
256	65 536	1 024	64:1
512	262 144	2 304	114:1
1024	1 048 576	5 120	205:1
2048	4 194 304	11 264	373:1
4096	16 777 216	24 576	683:1

Особенно заметен прирост эффективности с увеличением количества отсчётов.

Скачать конспект в pdf: FFT Lecture – V.V. Leonidov.pdf

Vladimir Leonidov

Visit Card