|
|
Писал кому-то.
БПФ - быстрый алгоритм вычисления ДПФ последовательности конечной длины.
Или периодической последовательности на одном ее периоде.
Обычно считается, что длина последовательности при БПФ - степень 2.
Хотя это не обязательно.
ДПФ: X(jW)=SUMMA(x(n)*exp(-j*W*n)). SUMMA - по n.
W - "цифровая частота" W=2*pi*f/Fd. Fd-частота дискретизации, f - обычная частота.
Ну для периодической последовательности все просто.
Ее спектр дискретен и дискреты кратны величине, обратно пропорциональной
длительности периода последовательности (ряд Фурье).
Если длительность периода в отсчетах =N, то его же длительность в единицах времени =
N/Fd. Т.е. частоты f кратны величине 1/(N/Fd)=Fd/N и, соответственно их можно записать в виде:
f=k*Fd/N. k - целые числа. Соответственно формула для X[jW] примет вид:
X(j*2*pi*k/N)=SUMMA(x(n)*exp(-j*2*pi*k*n/N)).
Упрощенная запись: X(k)= SUMMA(x(n)*exp(-j*2*pi*k*n/N)).
Собственно это исходная формула для БПФ.
Видно, что выражение периодично, если можно так сказать, и по n и по k с одинаковым периодом=N.
Что и понятно - исходная последовательность, во-первых, дискретна, что дает периодичность спектра,
т.е. периодичность по k, во-вторых периодична по времени.
Та же самая формула будет и для произвольной, не обязательно периодичной,
последовательности конечной длины N.
В этом случае можно провести аналогию с теоремой Котельникова.
Т. Котельникова, как известно, говорит, что сигнал, имеющий конечную ширину спектра F,
м.б. восстановлен по своим выборкам, следующим с интервалом не реже 1/2*F.
А для комплексных, точнее аналитических сигналов, с интервалом не реже 1/F.
Ни аналитических, ни комплексных сигналов в природе нет, но это другой вопрос.
Теперь в качестве сигнала возьмем его частотный спектр, а в качестве "спектра"
возьмем сам сигнал (последовательность). Т.е. заменим частоту на время и время на частоту.
Применим т. Котельникова к полученному.
Поскольку сигнал (последовательность) конечной длины, то ее спектр, а он-то аналитический,
может быть точно восстановлен по своим выборкам, следующим с интервалом не реже 1/T.
T - длительность последовательности во временных единицах и равна N/Fd.
Естественный вывод из этого - достаточно вычислять отсчеты спектра в точках,
кратных Fd/N.
Т.е. опять f=k*Fd/N. Подставив в формулу для X(jW), получим уже известное выражение для X(k).
Но в нем, как уже говорилось, периодичность по n и k. И период равен N.
Т.е. Фурье над непериодической последовательностью конечной длины N полностью эквивалентно
Фурье над другой, периодической последовательностью.
И эта другая последовательность получается периодическим, с периодом N, продолжением исходной.
Теперь собственно БПФ. Их много. Остановлюсь на одном.
Длина последовательности, N, равна какой-то степени 2-ки.
Исходно. X(k)= SUMMA(x(n)*exp(-j*2*pi*k*n/N)).
Сумма по всем n, от 0 до N-1.
Для расчета только одного значения X нужно выполнить N операций.
Каждая операция - комплексное умножение с комплексным сложением.
Соответственно для вычисления всех значений X понадобится N^2 таких операций.
Разделим на 2 суммы - одна по четным n, другая по нечетным.
Получим:
X(k)=SUMMA(x(2*n)*exp(-j*2*pi*k*2*n/N)) + SUMMA(x(2*n+1)*exp(-j*2*pi*k*(2*n+1)/N))
Обе суммы тоже по n , но от 0 до (N/2)-1.
Похимичим над выражением.
X(k)=SUMMA(x(2*n)*exp(-j*2*pi*k*n/(N/2))) + SUMMA(x(2*n+1)*exp(-j*2*pi*k*/(N/2))*exp(-j*2*pi*k/N))
Вынесем,как не зависящую от n, exp(-j*2*pi*k/N) за знак второй суммы.
Обозначим xe(n)=x(2*n) и xo(n)=x(2*n+1).
Подставим, получим:
X(k)=SUMMA(xe(n)*exp(-j*2*pi*k*n/(N/2))) +exp(-j*2*pi*k/N)* SUMMA(xo(n)*exp(-j*2*pi*k*/(N/2)))
Что такое каждая из сумм? Вид каждой из них идентичен преобразованию Фурье
(см. исходную формулу для X(k)), но над последовательностями длины N/2.
Преобразованиями Фурье в 2 раза меньшей длины они и являются.
Соответственно каждая из сумм периодична с периодом N/2.
Для упрощения записи обозначим первую сумму как XE(k), а вторую как XE(k).
Получим X(k)=XE(k) +exp(-j*2*pi*k/N)*XO(k). Нас интересуют значения k от 0
до N-1, в то же время и XE(k) и XO(k) периодичны с периодом N/2.
Возьмем k от 0 до (N/2)-1. Выражение для X(k) при таких k только что приведено.
Найдем выражение для X(k+N/2). В силу периодичности XE,XO получим
X(k+N/2)=XE(k)+exp(-j*2*pi*(k+N/2)/N)*XO(k)=XE(k)-exp(-j*2*pi*k/N)*XO(k).
exp(-j*2*pi*k/N) обозначим как WN(k).
Получим пару выражений:
X(k)=XE(k)+WN(k)*XO(k).
X(k+N/2)=XE(k)-WN(k)*XO(k).
Эта пара называется бабочкой с прореживанием по времени по основанию 2.
А число бабочек равно N/2.
А сам метод деления на четные, нечетные - метод прореживания по времени,
по английски DIT (decimation in time).
XE и XO - ПФ длиной N/2. Соответственно для вычисления всех XE и всех XO
понадобится 2*((N/2)^2)=0.5*(N^2) операций. Плюс N операций сшивки XE и XO
(бабочки) для получения X. Величина 0.5*(N^2) существенно больше
N даже при небольших N. Т.е. практически двойная экономия при простом разделении
исходной последовательности на 2 и прямом, без каких-либо хитростей, вычислении
ПФ этих новых последовательностей.
Но и к XE(k) и XO(k) можно применить аналогичную процедуру.
Т.е. разделить xe(n) на 2 последовательности, xee и xeo, одна из четных, другая из
нечетных отсчетов. Аналогично xo(n) на xoe и xoo.
xee=xe(2*l),xeo=xe(2*l+1),xoe=xo(2*l),xoo=xo(2*l+1). l от 0 до (N/4)-1.
Каждая из этих последовательностей имеет длину N/4.
И получим выражения для XE,XO, аналогичные выражению для X.
XE[m]=XEE[m]+WN2[m]*XEO[m];
XE[m+N/4]=XEE[m]-WN2[m]*XEO[m];
XO[m]=XOE[m]+WN2[m]*XOO[m];
XO[m+N/4]=XOE[m]-WN2[m]*XOO[m];
XEE,XEO,XOE,XOO - преобразования Фурье длиной N/4 последовательностей xee,xeo,xoe,xoo
соответственно. m, как и l, от 0 до (N/4)-1.
WN2[m]=exp(-j*2*pi*k/(N/2))= exp(-j*2*pi*k*2/N))
Пара выражений для XE, как и пара для XO, та же бабочка.
Всего нужно N/4 бабочек для XE и столько же для XO.
Аналогично можно вычислить и XEE,XEO,XOE,XOO
Т.е. продолжать до тех пор, пока длины временных последовательностей не станут
равными 1. А число самих последовательностей будет равно N/2.
Вычисление X через XE,XO - это одна ступень БПФ.
Вычисление XE,XO через XEE,XEO,XOE,XOO - еще одна.
Вычисление XEE,XEO,XOE,XOO - еще одна. Ну и т.д.
Поскольку каждый раз временные последовательности делятся на две каждая, то общее число ступеней будет равно
log2(N). И на каждой ступени N/2 бабочек.
Т.е. общее число бабочек=0.5*N*log2(N).
Еще один распространенный метод - прореживания по частоте (DIF).
Последовательность тоже разбивается на 2, но по другому.
В первую последовательность входит первая половина исходной, во вторую - вторая.
А выражения будут для X(k) и X(2*k+1). Бабочка будет другой. Тоже по основанию 2, но DIF.
Общее число бабочек такое же, как для DIT.
Каких либо преимуществ у DIT перед DIF и наоборот нет. Но на каких-то процессорах
может оказаться удобнее DIT, а на других DIF.
Еще немного. Исходную последовательность можно делить не на 2,
а на 4 последовательности , например:
x0(n)=x(4*n), x1(n)=x(4*n+1),x2(n)=x(4*n+2),x3(n)=x(4*n+3).
Далее, по аналогии можно проделать преобразования и получить БПФ
по основанию 4 с прореживанием по времени. Выражения будут для X(k),X(k+N/4),
X(k+N/2),X(k+3*N/4). И бабочка будет состоять из 4 выражений. Так называемая бабочка по основанию 4
с прореживанием по времени, в данном случае.
БПФ по основанию 4 обычно процентов на 20 быстрее, чем по основанию 2.
Но N должно быть степенью 4.
Однако, если N степень 2, но не степень 4, БПФ по основанию 4 тоже можно применить.
Например так. Исходную последовательность делим на 2, xe (четные), xo (нечетные).
Длина и xe и xo уже степень 4. Сл-но их спектры XE,XO м.б. вычислены через БПФ по основанию 4,
а сшивку их для получения X - бабочками по основанию 2. В результате будет Фурье по смешанному
основанию, 2 и 4. Аналогично можно рассмотреть БПФ по основанию 4, но DIF.
Можно ввести и другие основания. 3,5,....
Все вышенаписанное предполагает, что длина последовательности, N, разлагается на сомножители.
Т.е. N - непростое число.
Для простых N тоже существуют быстрые алгоритмы.
Из приведенной выше бабочки видно, что как конечный, так и промежуточные результаты - комплексные.
Ну и, чтобы получить хорошую структуру, предполагают, что и входные данные тоже комплексные.
Даже если они действительные их считают комплексными, но с нулевой мнимой частью.
Однако можно воспользоваться тем, что для действительных сигналов X(k)=X'(-k). Здесь апостроф -
комплексное сопряжение. Для дискретного сигнала, -k эквивалентно
N-k, т.е. частоты от 0 до N/2 - положительные, а от (N/2)+1 до N-1 - отрицательные.
И для действительного сигнала отсчеты спектра на частотах от 1 до (N/2)-1 имеют комплексно-сопряженные
отсчеты на частотах от N-1 до (N/2)+1 соответственно. Отсчеты на частотах 0 и N/2 - особые,
не имеющие комплексно-сопряженных.
Этим можно воспользоваться при вычислении спектра действительного сигнала.
Исходную действительную последовательность x(n) длиной N, n от 0 до N-1, представляют как комплексную
последовательность сл. образом длиной N/2:
y(m)=x(2*m)+j*x(2*m+1)=xe(m)+j*xo(m).
m от 0 до (N/2)-1. Т.е. действит. часть комплексной последовательности - четные
отсчеты исходной, мнимая - нечетные.
И уже над y(m) взять БПФ, но, разумеется, половинной длины. Затем, воспользовавшись
комплексной сопряженностью, из спектра Y(k) можно выделить спектры XE(k),XO(k).
А затем воспользоавшись DIT бабочкой по основанию 2 сшить эти 2 спектра и получить
спектр X(k). Такой способ дает почти двухкратную экономию по сравнению с представлением
действит. сигнала, как комплексного той же длины, но с 0 мнимой частью.
E-mail: 
info@telesys.ru
