Главная » Программирование звука » Программирование БПФ

0

Воспользовавшись   приведенными   выше   идеями,   можно   получить   рекурсив-

ную реализацию преобразования, которая будет похожа на приведенную ниже.

Листинг 24.2. Медленная реализация быстрого преобразования

Фурье

static void SlowFftRecursion(complex<double> *samples, int length, int start, int skip, complex<double> *result) {

if (length =1) {    // Простое одноточечное  БПФ.

*result = samples[start];

return;

}

// Вычисление двух  БПФ "половинного"

// размера.

SlowFftRecursion(samples,length/2,start,skip*2 , result); SlowFftRecursion(samples,length/

2,start+skip,skip*2,result+length/2);

// Вычисление  суммы  и  разности пар.

for(int j=0; j<length/2; j++) {

// Умножение второй  части  на   коэффициент

// фазового сдвига.

complex<double> t = result[j+length/2] * polar(1.0,-2*PI*j/

length);

// Вычитание.

result[j + length/2] = result[j] t;

// Добавление.

result[j] += t;

}

}

void SlowForwardFft(complex<double> *samples, int length, complex<double>*result) {

SlowFftRecursion(samples,length,0,1,result);

}

Чтобы   проверить   эту   функцию,   необходимо   сравнить   ее   выходные   данные с данными реализации ДПФ, приведенной ранее.

Есть несколько способов ускорения работы данной программы. Наиболее очевид-

ный устранить повторяющиеся вызовы polar во внутреннем цикле (polar (1, x) эквивалентно  eix.  Следующий  шаг   полностью  убрать  рекурсию,  однако  для  этого придется поработать.

Если  вы  внимательно  изучите  рекурсию  в  SlowFftRecursion,  то  увидите, что  рекурсии  нижнего  уровня  (1-точечные  БПФ)  просто  копируют  данные  в  определенном  порядке,  в  то  время  как  основная  работа  уже  проделана.  Чтобы  избавиться  от  этой  рекурсии,  сначала  придется  придумать  другой  способ  реализации переупорядочивания.

После  того  как  переупорядочивание  будет  завершено,  вы  можете  эмулировать рекурсивные  вызовы  в  несколько  ином  порядке.  Сначала  выполняем  2-точечное БПФ  для  последовательных  пар,  затем  вычисляем  4-точечные  БПФ  для  каждой группы  из  четырех  отсчетов  и  т.д.  Конечная  нерекурсивная  реализация  будет  выглядеть приблизительно так:

for(int halfSize=1; halfSize < length; halfSize *= 2) {

for(int fftStep = 0,fftStep < halfSize; fftStep++) {

for(int i = fftStep; i < length; i += 2*halfSize) {

complex<double> t = currentPhaseShift *

samples[i+halfSize];

samples[i+halfSize] = samples[i] t;

samples[i] += t;

}

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

Листинг 24.4. Программа fft.h

#include <complex>

// Длина   должна быть  степенью  двойки.

void ForwardFft(complex<double> *, int length);

void InverseFft(complex<double> *, int length);

// Медленные  версии,  экспортированные для сравнения

// и  экспериментов.

void ForwardDft(complex <double> *samples,

int length, complex<double>*result);

void SlowForwardFft(complex<double> *samples, int length, complex<double>*result);

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

Листинг 24.5. Программа fft.cpp

#include "fft.h"

#include <cmath>

#include <iostream>

const double PI = 3.14159265358979323846264338;

Rearrange(samples,length);

for(int halfSize=1; halfSize < length; halfSize *= 2) { complex<double> phaseShiftStep = polar(1.0,-PI/halfSize); complex<double> currentPhaseShift(1,0);

for(int fftStep = 0; fftStep < halfSize; fftStep++) {

for(int i=fftStep; i < length; i += 2*halfSize) {

complex<double> t = currentPhaseShift *

samples[i+halfSize];

samples[i+halfSize] = samples[i] t;

samples[i] += t;

}

currentPhaseShift *= phaseShiftStep;

}

}

}

// Используя ряд  простых свойств  комплексных чисел,

// мы можем  вычислить обратное  БПФ путем транспонирования

// отсчетов до   и  после  вычисления прямого БПФ.

// Заодно здесь  они   домножаются  на   1/N.

void InverseFft(complex<double> *samples, int length ) {

for(int i=0; i<length; i++)

samples[i] = conj(samples[i]); ForwardFft(samples,length);

for(int i = 0; i<length; i ++ )

samples[i] = conj(samples[i]) /

static_cast<double>(length);

}

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

Сложность  заключается  в  анализе  двоичного  значения  индекса  массива  для каждого элемента. Вспомним, что четное число содержит 0 в младшем бите и что четные  элементы  будут  размещаться  в  первой  половине  массива,  для  которой  индекс содержит 0 в старшем бите. Первый шаг перестановки заключается в перемещении младшего бита индекса массива в старший бит.

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

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

Листинг 24.6. Перестановка отсчетов для быстрой реализации БПФ

static void Rearrange(complex<double> *samples, int length) {

static int rearrangeSize = 0;   // Размер    таблицы перестановки.

static int *rearrange = 0;

if (rearrangeSize != length) {

// Производим перерасчет  каждый   раз при   изменении  размера.

if (rearrange) delete [] rearrange;

rearrange = new int[length];

// Заполняем  назначение   каждой величины.

rearrange[0] = 0;

for(int limit=1, bit=length/2; limit<length; limit <<= 1, bit>>=1 )

for(int i=0;i<limit;i++)

rearrange[i+limit] = rearrange[i] + bit;

// Помещаем   0 в те элементы,  которые остаются

// без изменений.

// Также обнуляем один из  элементов  так,

// чтобы  каждый   обмен происходил  только  один раз.

for(int i=0; i<length; i++) {

if (rearrange[i] == i) rearrange[i] = 0;

else rearrange[ rearrange[i] ] = 0;

}

rearrangeSize = length;

}

// Используем таблицу перестановки  для  перестановки  элементов.

// Нулевые индексы просто  пропускаются.

complex <double> t;

for (int i=0; i<length; i++)

if (rearrange[i]) {   // Этот элемент  надо  переставлять?

t = samples[i];    // Да, переставляем.

samples[i] = samples[ rearrange[i] ];

samples[ rearrange[i] ] = t;

}

}

?

Скорость

При разработке сложных систем обработки звука основной проблемой является скорость. Например, если вам нужно использовать 1024-точечные БПФ с качеством звука,  соответствующим  компакт-диску,  то  для  обеспечения  непрерывного  потока данных  необходимо  вычислять  каждое  БПФ  менее  чем  за  6  миллисекунд.  Этому условию  достаточно  сложно  следовать.  Подпрограмма  ForwardFft,  показанная мною  ранее,  требует,  например,  36  миллисекунд  для  вычисления  1024-точечного БПФ  на  процессоре  Intel 486/66.  (Для  сравнения,  подпрограмма  SlowForwardFft работает  163  миллисекунды  на  486/66,  в  то  время  как  подпрограмма  ForwardDft  почти  70  секунд!)  200-мегагерцовый  PowerPC  603e  способен  выполнить  такое  БПФ  за  6  миллисекунд,  но  при  этом  не  остается  времени  на  дополнительную обработку, которая также может потребоваться.

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

короткие   БПФ.   Два   512-точечных   БПФ   отнимают   меньше   времени,   чем   одно

1024-точечное  БПФ.  И,  наконец,  вместо  арифметики  с  плавающей  точкой  вы  можете  использовать  целочисленную  арифметику  или  арифметику  с  фиксированной точкой.

Эксперименты с БПФ

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

Листинг 24.7. Программа fftmain.cpp

#include "fft.h"

#include <iostream>

#include <iomanip>

#include <cstdlib>

int main(int argc, char **argv) {

int length = 8;

if (argc > 1) length = atoi(argv[1]);

if (((~length+1)&length) != length) {

// Длина   должна быть  степенью  2!

cerr << "Length must be a power of two!\n";

exit(1);

}

complex<double> *input = new complex<double>[length];

for(int i=0; (i < length)&&(!cin.eof()); i++)

cin >> input[i];

complex<double> *output = new complex<double>[length];

for(int i=0;i<length;i++) output[i] = input[i]; ForwardFft(output,length);

// Раскомментировать  следующий раздел  для проверки?

#if 0

complex<double> *compare = new complex<double>[length];

ForwardDft(input,length,compare);

double maxerr=0; int maxindex=0;

for(int i=0;i<length;i++) {

// Вычисляем относительную   ошибку

double error = abs(output[i] compare[i]) /

abs(compare[i]);

if(error > maxerr) {

maxindex=i; maxerr=error;

}

}

// Максимальная относительная  ошибка.

cerr << "Maximum relative error:" << maxerr << "\n";

#endif

for(int i=0;i<length;i++)

cout << setprecision(20) << output[i] << "\n";

return 0;

}

Эти   две   программы   используют   стандартную   библиотеку   C++   для   чтения комплексных   чисел   из   стандартного   потока   ввода.   Стандартная   библиотека   C++ работает  с  комплексными  числами  двух  форматов.  Если  первый  символ   круглая  скобка,  ожидается  ввод  заключенной  в  скобки  пары  чисел,  представляющих действительную  и  мнимую  части.  B  противном  случае  считывается  одно  число с  плавающей  точкой,  представляющее  действительную  часть,  а  мнимая  часть  полагается  равной  нулю.  Например,  и  (1,  0),  и  1,0  будут  восприняты  как  комплексное число 1 + 0i.

Листинг 24.8. Программа ifftmain.cpp

#include "fft.h"

#include <iostream>

#include <iomanip>

#include <cstdlib>

int main(int argc, char **argv) {

int length = 8;

if (argc > 1) length = atoi(argv[1]);

if (((~length+1)&length) != length) {

// Длина   должна быть степенью  числа  2.

cerr << "Length must be a power of two!\n";

exit(1);

}

complex<double> *a = new complex<double>[length];

for(int i=0; (i < length)&&(!cin.eof()); i++)

cin >> a[i]; InverseFft(a,length); for(int i=0;i<length;i++)

cout << setprecision(20) << a[i] << "\n";

return 0;

}

Источник: Кинтцель Т.  Руководство программиста по работе со звуком = A Programmer’s Guide to Sound: Пер. с англ. М.: ДМК Пресс, 2000. 432 с, ил. (Серия «Для программистов»).

По теме:

  • Комментарии