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

Содержание
Введение
В предыдущем разделе был подробно рассмотрен алгоритм БПФ с прореживанием по времени.
Рассмотрим еще один способ разделения и объединения, который носит название прореживание по частоте.
Алгоритм БПФ с прореживанием по частоте
Пусть имеется $N$ отсчетов входного сигнала $s(n)$, $n = 0 \ldots N-1$, при этом $N$ представляет собой целую степень двойки $N = 2^L$.
В алгоритме БПФ с прореживанием по времени производилось разделение исходного сигнала в соответствии с двоично-инверсной перестановкой. Таким образом, мы получили первую и вторую половины дискретного преобразования Фурье (ДПФ).
В алгоритме с прореживанием по частоте исходный сигнал $s(n)$, $n = 0 \ldots N-1$, делится пополам, т.е. $s_0(m) = s(m)$ и $s_1(m) = s\left(m+\frac{N}{2}\right)$, $m = 0 \ldots \frac{N}{2} - 1$.
Тогда ДПФ сигнала $s(n)$ может быть записано в виде:
$$ S(k) = \sum_{n = 0}^{N - 1} s(n) \cdot W_N^{n\cdot k} =\sum_{m = 0}^{\frac{N}{2} - 1} \left[ s(m) \cdot W_N^{m\cdot k} + s\left( m + \frac{N}{2} \right) \cdot W_N^{\left(m + \frac{N}{2} \right)\cdot k} \right] = $$ $$ =\sum_{m = 0}^{\frac{N}{2} - 1} \left[ s(m) \cdot W_N^{m\cdot k} + s\left( m + \frac{N}{2} \right) \cdot W_N^{m \cdot k} \cdot W_N^{\frac{N}{2} \cdot k}\right] = $$ $$ =\sum_{m = 0}^{\frac{N}{2} - 1} \left[ s(m) + s\left( m + \frac{N}{2} \right) \cdot W_N^{\frac{N}{2} \cdot k}\right]\cdot W_N^{m \cdot k}, \textrm{ }k=0 \ldots N-1. $$
(1)
Рассмотрим выражение (1) для спектральных отсчетов $S(2\cdot p)$, $p = 0 \ldots \frac{N}{2}-1$, с четными индексами:
\begin{equation} \label{fft_dec_in_freqs:eq2} S(2 \cdot p) = \sum_{m = 0}^{\frac{N}{2} - 1} \left[ s(m) + s\left( m + \frac{N}{2} \right) \cdot W_N^{\frac{N}{2} \cdot 2 \cdot p}\right]\cdot W_N^{m \cdot 2\cdot p}, \textrm{ }p = 0 \ldots \frac{N}{2}-1. \end{equation}
(2)
Учтем в выражении (2), что $W_N^{\frac{N}{2} \cdot 2 \cdot p} = W_N^{N \cdot p} = 1$, а также $W_N^{m \cdot 2\cdot p} = W_{\frac{N}{2}}^{m \cdot p}$, тогда получим:
\begin{equation} \label{fft_dec_in_freqs:eq3} S(2 \cdot p) = \sum_{m = 0}^{\frac{N}{2} - 1} \underbrace{\left[ s(m) + s\left( m + \frac{N}{2} \right) \right]}_{x_0(m)}\cdot W_{\frac{N}{2}}^{m \cdot p} = \underbrace{\sum_{m = 0}^{\frac{N}{2} - 1} x_0(m) \cdot W_{\frac{N}{2}}^{m \cdot p}}_{\frac{N}{2}-\textrm{точечное ДПФ}}, \textrm{ } p = 0 \ldots \frac{N}{2}-1. \end{equation}
(3)
Таким образом, спектральные отсчеты $S(2\cdot p)$, $p = 0 \ldots \frac{N}{2}-1$, с четными индексами есть $\frac{N}{2}-$ точечное ДПФ сигнала $x_0(m) = s(m) + s\left( m + \frac{N}{2} \right)$.
Аналогично рассмотрим выражение (1) для спектральных отсчетов $S(2\cdot p+1)$, $p = 0 \ldots \frac{N}{2}-1$, с нечетными индексами:
\begin{equation} \label{fft_dec_in_freqs:eq4} S(2 \cdot p+1) = \sum_{m = 0}^{\frac{N}{2} - 1} \left[ s(m) + s\left( m + \frac{N}{2} \right) \cdot W_N^{\frac{N}{2} \cdot (2 \cdot p+1)}\right]\cdot W_N^{m \cdot (2 \cdot p+1)}, \textrm{ }p = 0 \ldots \frac{N}{2}-1. \end{equation}
(4)
Учтем в выражении (4), что $W_N^{\frac{N}{2} \cdot (2 \cdot p+1)} = W_N^{\frac{N}{2} \cdot 2 \cdot p} \cdot W_N^{\frac{N}{2}} = -1$, а также, что $W_N^{m \cdot (2 \cdot p+1)} =W_N^{m} \cdot W_N^{m \cdot 2 \cdot p} = W_N^{m} \cdot W_{\frac{N}{2}}^{m \cdot p}$. Тогда выражение (4) можно представить в виде:
\begin{equation} \label{fft_dec_in_freqs:eq5} S(2 \cdot p+1) = \sum_{m = 0}^{\frac{N}{2} - 1} \underbrace{\left[ s(m) - s\left( m + \frac{N}{2} \right) \right] \cdot W_N^m}_{x_1(m)} \cdot W_{\frac{N}{2}}^{m \cdot p} = \underbrace{\sum_{m = 0}^{\frac{N}{2} - 1} x_1(m) \cdot W_{\frac{N}{2}}^{m \cdot p}}_{\frac{N}{2}- \textrm{точечное ДПФ}}, \textrm{ }p = 0 \ldots \frac{N}{2}-1. \end{equation}
(5)
Спектральные отсчеты $S(2\cdot p+1)$, $p = 0 \ldots \frac{N}{2}-1$, с нечетными индексами $2 \cdot p+1$ есть $\frac{N}{2}-$ точечное ДПФ сигнала $x_1(m)= \left[s(m) - s\left( m + \frac{N}{2} \right)\right] \cdot W_N^{m}$.
Таким образом, процедура разделения заключается в расчете сигналов половинной длительности $x_0(m)$ и $x_1(m)$, $m = 0 \ldots \frac{N}{2} - 1$. После чего можно заменить $N-$точечное ДПФ двумя $\frac{N}{2}-$точечными ДПФ сигналов $x_0(m)$ и $x_1(m)$.
Граф «бабочка» алгоритма БПФ с прореживанием по частоте
Процедура расчета сигналов половинной длительности
\begin{align} \label{fft_dec_in_freq:eqbutterfly} x_0(m) &= s(m) + s\left( m + \frac{N}{2}\right), & m &= 0 \ldots \frac{N}{2}-1; \nonumber \\ x_1(m) &= \left[ s(m) - s\left( m + \frac{N}{2}\right) \right] \cdot W_N^m, & m &= 0 \ldots \frac{N}{2}-1 \end{align}
(6)
может быть представлена в виде графа «бабочка», как это показано на рисунке 1.
Рисунок 1. Граф «бабочка» алгоритма БПФ с прореживанием по частоте
Отличие графа «бабочка» алгоритма с прореживанием по частоте от аналогичного графа для алгоритма с прореживанием по времени заключается в том, что умножение на поворотный коэффициент $W_N^m$ производится после вычитания ветвей. В графе «бабочка» алгоритма с прореживанием по времени умножение на поворотный коэффициент производилось до вычитания ветвей.
Полный граф алгоритма БПФ с прореживанием по частоте
Мы заменили расчет $N-$точечного ДПФ двумя $\frac{N}{2}-$точечными. В результате граф алгоритма БПФ с прореживанием по частоте для $N=8$ приведен на рисунке 2.
Рисунок 2. Граф алгоритма БПФ с прореживанием по частоте для $N=8$
При этом каждое из $\frac{N}{2}-$точечных ДПФ также может быть рассчитано с использованием алгоритма с прореживанием по частоте. В результате получим полный граф алгоритма БПФ с прореживанием по частоте, как это показано на рисунке 3.
Рисунок 3. Полный граф алгоритма БПФ с прореживанием по частоте для $N=8$
На первом этапе получаем $x_0(m)$ и $x_1(m)$, в соответствии с (6):
\begin{align} x_0(0) & = s(0) + s(4);\nonumber \\ x_0(1) & = s(1) + s(5);\nonumber \\ x_0(2) & = s(2) + s(6);\nonumber \\ x_0(3) & = s(3) + s(7);\nonumber \\ x_1(0) & = [s(0) - s(4)] \cdot W_8^0;\nonumber \\ x_1(1) & = [s(1) - s(5)]\cdot W_8^1 ;\nonumber \\ x_1(2) & = [s(2) - s(6)]\cdot W_8^2 ;\nonumber \\ x_1(3) & = [s(3) - s(7)]\cdot W_8^3 . \end{align}
(7)
Для расчета 4-точечных ДПФ сигналов $x_0(m)$ и $x_1(m)$, $m = 0 \dots 3$, снова используем алгоритм с прореживанием по частоте. Тогда получим сигналы:
\begin{align} x_{00}(0) & = x_0(0) + x_0(2);\nonumber \\ x_{00}(1) & = x_0(1) + x_0(3);\nonumber \\ \\ x_{01}(0) & = [x_0(0) - x_0(2)] \cdot W_4^0;\nonumber \\ x_{01}(1) & = [x_0(1) - x_0(3)] \cdot W_4^1;\nonumber \\ \\ x_{10}(0) & = x_1(0) + x_1(2);\nonumber \\ x_{10}(1) & = x_1(1) + x_1(3);\nonumber \\ \\ x_{11}(0) & = [x_1(0) - x_1(2)] \cdot W_4^0;\nonumber \\ x_{11}(1) & = [x_1(1) - x_1(3)] \cdot W_4^1.\nonumber \end{align}
(8)
После расчета двухточечных ДПФ на последнем уровне разбиения получаем переставленные спектральные отсчеты:
\begin{align} S(0) & = x_{00}(0) + x_{00}(1);\nonumber \\ S(4) & = [x_{00}(0) - x_{00}(1)] \cdot W_2^0;\nonumber \\ \\ S(2) & = x_{01}(0) + x_{01}(1);\nonumber \\ S(6) & = [x_{01}(0) - x_{01}(1)] \cdot W_2^0;\nonumber \\ \\ S(1) & = x_{10}(0) + x_{10}(1);\nonumber \\ S(5) & = [x_{10}(0) - x_{10}(1)] \cdot W_2^0;\nonumber \\ \\ S(3) & = x_{11}(0) + x_{11}(1);\nonumber \\ S(7) & = [x_{11}(0) - x_{11}(1)] \cdot W_2^0,\nonumber \end{align}
(9)
которые мы должны подвергнуть двоично-инверсной перестановке аналогично тому, как производится эта процедура в алгоритме с прореживанием по времени.
Выводы
В данном разделе мы рассмотрели алгоритм БПФ по основанию два с прореживанием по частоте. В следующем разделе мы рассмотрим вычислительную эффективность алгоритмов БПФ по основанию два и некоторые вопросы программной реализации алгоритмов БПФ.

Ваши комментарии, вопросы и пожелания вы можете оставить на форуме или в гостевой книге.
Список литературы
[1] James W. Cooley and John W. Tukey An Algorithm for the Machine Calculation of Complex Fourier Series. Mathematics of Computation, 1965 p. 297–301.

[2] Блейхут Р. Быстрые алгоритмы цифровой обработки сигналов. Москва, Мир, 1989.

[3] Нуссбаумер Г. Быстрое преобразование Фурье и алгоритмы вычисления сверток. Москва, Радио и связь, 1985.

[4] Оппенгейм А., Шафер Р. Цифровая обработка сигналов. Москва, Техносфера, 2006.

[5] Гольденберг Л.М., Матюшкин Б.Д., Поляк М.Н. Цифровая обработка сигналов. Москва, Радио и связь, 1985.

[6] Рабинер Л., Гоулд Б. Теория и применение цифровой обработки сигналов. Москва, Мир, 1978.

[7] Сергиенко А.Б. Цифровая обработка сигналов. Санкт-Петербург, Питер, 2002.

Oбнаружили ошибку в тексте? Выделите ее мышкой и нажмите