Алгоритм БПФ составной длины

Содержание

Введение

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

Пример расчета для N = 6

N = 6

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

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

Алгоритмы БПФ по основанию 2 как частный случай алгоритма составной длины

Список литературы

DSPL-2.0 — свободная библиотека алгоритмов цифровой обработки сигналов

Распространяется под лицензией LGPL v3

Страница проекта на SourceForge

Введение

В предыдущих разделах мы рассмотрели алгоритмы БПФ по основанию два с прореживанием по времени и с прореживанием по частоте. Данные алгоритм очень эффективны, но они имеют существенное ограничение: длина входных данных и выходного вектора должна быть целой степенью двойки N = 2^t .

В данном разделе мы рассмотрим алгоритмы БПФ для произвольной составной длины $N = M \cdot L$ , где и — произвольные положительные целые. Данный подход позволит построить эффективные процедуры дискретного преобразования Фурье практически для всех длин, а не только для N = 2^t . Кроме того, мы покажем, что рассмотренные ранее алгоритмы БПФ по основанию два являются частным случаем алгоритма БПФ составной длины.

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

Пусть длина ДПФ равна произведению $N =M \cdot L$ , где и — целые положительные числа.

Выражение для ДПФ сигнала s(n) , $n = 0\ldots N-1$ можно записать в виде:

equation 1

(1)

где $W_N = \exp(-j 2\pi /N)$ — основание поворотного коэффициента.

Для составной длины $N = M\cdot L$ можно разделить индексы и следующим образом:

equation 2

(2)

Тогда относительно индексов

,

,

и

выражение ДПФ (1) разделится на вложенные суммы:

equation 3

(3)

Раскроем скобки в показателе поворотных коэффициентов:

equation 4

(4)

Учтем, что $N = M \cdot L$ :

equation 5

(5)

тогда получаем выражение для ДПФ составной длины:

equation 6

(6)

Разделим входной сигнал длины $M\cdot L$ на сегменты x_m(p)

x_m(p)

, $m = 0\ldots M-1$ , $p = 0\ldots L-1$ получим внутреннюю сумму в виде

-точечного ДПФ:

equation 7

(7)

Таких

-точечных ДПФ будет

штук для каждого индекса $m = 0\ldots M-1$ .

После этого необходимо каждый из -точечных ДПФ y_m(q) умножить на соответствующие поворотные коэффициенты $W_N^{mq}$ и сформировать сигналов z_q(m)

equation 8

(8)

длительности

отсчетов. Тогда внешняя сумма представляет собой

-точечное ДПФ сигнала z_q(m)

z_q(m)

, которое возвращает один сегмент выходного спектра

equation 9

(9)

Всего таких сегментов будет

для каждого $q = 0 \ldots M-1$ .

Можно подвести итог: для расчета ДПФ составной длины $N = M\cdot L$ требуется ДПФ размера точек, умножений на поворотные коэффициенты и ДПФ размера точек.

Пример расчета для N = 6

N = 6

Рассмотрим пример для N=6 . Тогда M = 2 , L = 3 . Вектор входного сигнала имеет вид:

equation 10

(10)

Индексы $m = 0,\,1$ , $p = 0,\,1,\,2$ и согласно (6) получаем два набора сигналов x_m(p) в виде векторов $\mathbf{x}_0$ и $\mathbf{x}_1$ для индексов m=0 и m = 1 соответственно:

equation 11

(11)

Рассчитаем 2 ДПФ размера 3 точки для $\mathbf{x}_0$ и $\mathbf{x}_1$ . Получим два вектора $\mathbf{y}_0$ и $\mathbf{y}_1$ , хранящие значения y_m(q)

y_m(q)

:

equation 12

(12)

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

y_m(q)

в

z_q(m)

. Тогда получим 3 вектора $\mathbf{z}_q$ , $q =0,\,1,\,2$ :

equation 13

(13)

После этого применяем для каждого $\mathbf{z}_q$ , $q =0,\,1,\,2$ ДПФ размера 2 точки и получаем векторы $\mathbf{v}_q$ , хранящие значения V_q(r)

V_q(r)

:

equation 14

(14)

equation 15

(15)

Тогда окончательно спектр исходного вектора сигнала $\mathbf{s}$ равен:

equation 16

(16)

На этом расчет 6-точечного ДПФ окончен.

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

Для практической реализации можно интерпретировать алгоритм БПФ составной длины как двумерное преобразование, как это показано на рисунке 1 для приведенного выше примера для N = 6 .

ДПФ составной длины как двумерное преобразование при , ,

Рисунок 1. ДПФ составной длины как двумерное преобразование при N = 6

N = 6

,

M = 2

,

L = 3

Шаг1. Исходный вектор размерности $[6 \times 1]$ Преобразуется в матрицу $\mathbf{A}$ размерности $[2 \times 3]$ , содержащую 2 строки и 3 столбца.

Шаг 2. Транспонируем матрицу $\mathbf{A}$ и получаем матрицу $\mathbf{B} =\mathbf{A}^T$ размерности $[3 \times 2]$ .

Шаг 3. Рассчитываем 3-точечное ДПФ каждого столбца матрицы $\mathbf{B}$ и сохраняем результат в матрицу $\mathbf{Y}$ той же размерности $[3 \times 2]$ .

Шаг 4. Производим поэлементное умножение матрицы $\mathbf{Y}$ и предварительно подготовленной матрицы поворотных коэффициентов $\mathbf{W}$ (знак $\odot$ означает поэлементное умножение). На выходе имеем матрицу $\mathbf{Z}$ размерности $[3 \times 2]$ .

Шаг 5. Транспонируем матрицу $\mathbf{Z}$ и получаем $\mathbf{Q = Z}^T$ размерности $[2 \times 3]$ .

Шаг 6. Рассчитываем 2-точечное ДПФ каждого столбца матрицы $\mathbf{Q}$ и сохраняем результат в матрицу $\mathbf{V}$ той же размерности $[2 \times 3]$ .

Шаг 7. Наконец транспонируем матрицу $\mathbf{V}$ и преобразуем матрицу $\mathbf{P = V}^T$ в вектор 6-точечного ДПФ исходного сигнала.

Данная схема подходит для любого размера $N = M \cdot L$ и требует три операции транспонирования матриц размерности $[M \times L]$ и $[L \times M]$ , комплексных умножений $\mathbf{Y} \odot \mathbf{W}$ , а также ДПФ размера точек и ДПФ размера точек.

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

Следующий скрипт dft_composite.m реализует расчет ДПФ составной длины в соответствии с приведенным алгоритмом. Скрипт можно запустить в пакете GNU Octave или в Matlab.

dft_composite.m


clear all; close all; clc;


M = 2;
L = 3;
N = M*L;


% Заполняю матрицу W [L x M] поворотных к-тов составного FFT 
W = zeros(L, M);
for l = 0:L-1
  for m = 0:M-1
    W(l+1, m+1) = exp( -j * 2 * pi * l * m / N);
  end
end

% Случайный комплексный сигнал 6 отсчетов
x = randn(N,1) + 1i * randn(N,1);

% Перевожу входной сигнал в матрицу [M x L],
A = reshape(x, M, L);

% транспонирую матрицу A и получаю матрицу B [L x M] 
B = transpose(A);

% FFT по столбцам матрицы B 
Y = fft(B);

% Поэлементное умножение матрицы Y и 
% матрицы поворотных к-тов W
Z = Y .* W;

% Транспонирую матрицу Z и получаю матрицу Q размера [M x L]
Q = transpose(Z);

% FFT по столбцам матрицы Q. Получим матрицу V [M x L]
V = fft(Q);

% Транспонируем матрицу V и получаем P [L x M]
P = transpose(V);

%стыкуем столбцы матрицы P друг за другом и получаем результат FFT
X = reshape(P, M*L, 1);

err = max(abs(X - fft(x)));
fprintf('max error = %.3e\n', err);

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

Во всех вышеприведенных примерах для N = 6 мы представляли $N = M\cdot L$ , где M = 2 , L = 3 . Однако мы можем изменить порядок представления, т.е. задать M = 3 , а L = 2 . В этом случае получим равносильный алгоритм, в котором на первом этапе потребуется 3 ДПФ размера 2 точки, а на втором этапе — 2 ДПФ размера 3 точки.

Тогда произведем разложение для сигнала (10) при изменении порядка представления.

Индексы $m = 0,\,1,\,2$ , $p = 0,\,1,$ и согласно (6) получаем три набора сигналов x_m(p) в виде векторов $\mathbf{x}_0$ , $\mathbf{x}_1$ и $\mathbf{x}_2$ для индексов $m=0, \, 1,\,2$ соответственно:

equation 17

(17)

Рассчитаем 3 ДПФ размера 2 точки для $\mathbf{x}_0$ , $\mathbf{x}_1$ и $\mathbf{x}_2$ . Получим три вектора $\mathbf{y}_0$ , $\mathbf{y}_1$ , и $\mathbf{y}_2$ , хранящие значения y_m(q) :

equation 18

(18)

equation 19

(19)

Следующий шаг: умножение на поворотные коэффициенты и переиндексация y_m(q) в z_q(m) . Тогда получим 2 вектора $\mathbf{z}_q$ , $q =0,\,1$ :

equation 20

(20)

После этого применяем для каждого $\mathbf{z}_q$ , $q =0,\,1,$ ДПФ размера 3 точки и получаем векторы $\mathbf{v}_q$ , хранящие значения V_q(r)

V_q(r)

:

equation 21

(21)

тогда окончательно спектр исходного вектора сигнала $\mathbf{s}$ равен:

equation 22

(22)

На этом расчет 6-точечного ДПФ окончен. Алгоритм при изменении порядка разложения на составные множители также можно представить в виде двумерного преобразования, как это показано на рисунке 2.

ДПФ составной длины как двумерное преобразование при , , .

Рисунок 2. ДПФ составной длины как двумерное преобразование при N = 6

N = 6

,

M = 3

,

L = 2

.

Можно заметить, что размеры матриц на рисунке 2 являются транспонированными по отношению к матрицам на рисунке 1. C точки зрения вычислительных операций, оба алгоритма эквивалентны.

Алгоритмы БПФ по основанию 2 как частный случай алгоритма составной длины

Ранее мы рассмотрели алгоритмы БПФ по основанию два: с прореживанием по времени и с прореживание по частоте для расчета ДПФ размера N = 2^t .

В данном случае составное число, которое можно представить как $N = M \cdot L$ , где M = 2 , $L = 2^{t-1}$ и тогда применив алгоритм составной длины получим БПФ по основанию 2 с прореживанием по времени, как это показано на рисунке 3 для N = 8 , M = 2 , L = 4 .

Алгоритм БПФ с прореживанием по времени, как частный случай алгоритма составной длины. , ,

Рисунок 3. Алгоритм БПФ с прореживанием по времени, как частный случай алгоритма составной длины. N = 8

N = 8

,

M = 2

,

L = 4

Процедура «прореживание по времени» возвращает два столбца матрицы $\mathbf{B}$ . Далее мы выполняем два ДПФ размера 4 точки над каждым столбцом матрицы $\mathbf{B}$ и получаем матрицу $\mathbf{Y}$ . После этого производим поэлементное умножение $\mathbf{Y}$ на матрицу поворотных коэффициентов $\mathbf{Y}$ . Выходной граф бабочка на структуре алгоритм БПФ с прореживанием по времени описывают 4 ДПФ размера 2 точки по столбцам транспонированной матрицы $\mathbf{Z}^T$ , в результате получаем матрицу переставленных отсчетов спектра V_q(r) = S(rL + q) .

Заметим, что каждое из 4-точечных ДПФ мы также можем рассчитывать с использованием алгоритма составной длины, применив который мы получим полный граф алгоритма с прореживанием по времени.

Если же мы изменим порядок разложения, т.е. для N = 2^t зададим $M = 2^{t-1}$ , L = 2 , то получим алгоритм БПФ по основанию 2 с прореживанием по частоте, как это показано на рисунке 4.

Алгоритм БПФ с прореживанием по частоте, как частный случай алгоритма составной длины. , ,

Рисунок 4. Алгоритм БПФ с прореживанием по частоте, как частный случай алгоритма составной длины. N = 8

N = 8

,

M = 4

,

L = 2

Первый каскад «бабочек» рассчитывает 4 ДПФ размера 2 точки и выводит матрицу $\mathbf{Y}^T$ , которая умножается на транспонированную матрицу поворотных коэффициенты $\mathbf{W}^T$ и получаем транспонированную матрицу $\mathbf{Z}^T$ . После этого рассчитываем два ДПФ размера 4 точки и возвращаем матрицу спектра $\mathbf{V}$ . На последнем этапе транспонируем матрицу $\mathbf{V}$ и получаем результат ДПФ.

Производя расчет 4-точечных ДПФ как ДПФ составной длины получаем алгоритм БПФ с прореживанием по частоте.

Таким образом, мы показали, что рассмотренные алгоритмы по основанию два это частный случай алгоритма составной длины, когда все множители равны двойке. Однако использование данного подхода позволяет получить высокоэффективные алгоритмы для ДПФ произвольных составных длин . Например ДПФ размера N = 100 можно представить как $N = 5 \cdot 20$ , и для расчета потребуется 5 ДПФ размера 20 точек и 20 ДПФ размера 5 точек. Причем каждое 20-точечное ДПФ также можно представить в виде составного алгоритма с множителями 5 и 4.

Таким образом можно заключить, что для эффективной реализации алгоритма составной длины необходимо иметь процедуры коротких ДПФ простых длин (2, 3, 5, 7, 11 и т.д.), тогда можно будет построить эффективные алгоритмы БПФ практически для всех размеров входных векторов. Эффективные алгоритмы простых длин были разработаны Шмуэлем Виноградом и опубликованы в 1978 году. Данные алгоритмы мы рассмотрим в следующем разделе.

Список литературы

[1] James W. Cooley and John W. Tukey An Algorithm for the Machine Calculation of Complex Fourier Series. Mathematics of Computation, Vol. 19, no. 2, April 1965, pp. 297-301.

[2] Оппенгейм А., Шаффер Р. Цифровая обработка сигналов. Москва, Техносфера, 2012. 1048 с. ISBN 978-5-94836-329-5

[3] Bracewell R. The Fourier Transform and Its Applications McGraw-Hills, 1986, 474 c. ISBN 0-07-007-015-6

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

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

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

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

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

Последнее изменение страницы: 12.05.2022 (19:41:47)

Страница создана Latex to HTML translator ver. 5.20.11.14

© Бахурин Сергей 2015 - 2022. Все права защищены.
Копирование материалов сайта без разрешения автора запрещено.