Програмування БПФ

Скориставшись наведеними вище ідеями, можна отримати рекурсією-

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

Лістинг 242 Повільна реалізація швидкого перетворення

Фурє

static void SlowFftRecursion(complex&ltdouble&gt *samples, int length, int start, int skip, complex&ltdouble&gt *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&ltlength/2 j++) {

/ / Множення другої частини на коефіцієнт

/ / Фазового зсуву

complex&ltdouble&gt t = result[j+length/2] * polar(10,-2*PI*j/

length)

/ / Вирахування

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

/ / Додавання

result[j] += t

}

}

void SlowForwardFft(complex&ltdouble&gt *samples, int length, complex&ltdouble&gt*result) {

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

}

Щоб перевірити цю функцію, необхідно порівняти її вихідні дані з даними реалізації ДПФ, наведеної раніше

Є кілька способів прискорення роботи даної програми Найбільш очевид-

ний усунути повторювані виклики polar у внутрішньому циклі (polar (1, x) еквівалентноeix Наступний крок повністю прибрати рекурсію, однак для цього доведеться попрацювати

Якщо ви уважно вивчите рекурсію в SlowFftRecursion, то побачите, що рекурсії нижнього рівня (1-точкові ШПФ) просто копіюють дані в певному порядку, в той час як основна робота вже зроблена Щоб позбутися цієї рекурсії, спочатку доведеться придумати інший спосіб реалізації переупорядочивания

Після того як переупорядочивание буде завершено, ви можете емулювати рекурсивні виклики в дещо іншому порядку Спочатку виконуємо 2-точкове ШПФ для послідовних пар, потім обчислюємо 4-точкові БПФ для кожної групи з чотирьох відліків і тд Кінцева нерекурсівние реалізація буде виглядати приблизно так:

for(int halfSize=1 halfSize &lt length halfSize *= 2) {

for(int fftStep = 0,fftStep &lt halfSize fftStep++) {

for(int i = fftStep i &lt length i += 2*halfSize) {

complex&ltdouble&gt t = currentPhaseShift *

samples[i+halfSize]

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

samples[i] += t

}

Результатом всієї цієї роботи є наступний модуль, який експортує всього чотири функції, причому дві з них є повільними версіями, включеними тільки для порівняння Найшвидша версія виконує всю роботу на місці розташування даних і, отже, вимагає створення одного масиву

Лістинг 244 Програма ffth

#include &ltcomplex&gt

/ / Довжина повинна бути ступенем двійки

void ForwardFft(complex&ltdouble&gt *, int length)

void InverseFft(complex&ltdouble&gt *, int length)

/ / Повільні версії, експортовані для порівняння

/ / І експериментів

void ForwardDft(complex &ltdouble&gt *samples,

int length, complex&ltdouble&gt*result)

void SlowForwardFft(complex&ltdouble&gt *samples, int length, complex&ltdouble&gt*result)

Якщо ви зрозуміли все вищесказане, то для вас цей код напрочуд простий Єдина додаткова оптимізація проведена мною для реалізації способу обчислення коефіцієнтів зсуву по фазі Так як кожен наступний коефіцієнт є всього лише ступенем попереднього, мені необхідно використовувати тільки одне множення у внутрішньому циклі для обчислення наступного коефіцієнта зсуву по фазі

Лістинг 245 Програма fftcpp

#include &quotffth&quot

#include &ltcmath&gt

#include &ltiostream&gt

const double PI = 314159265358979323846264338

Rearrange(samples,length)

for(int halfSize=1 halfSize &lt length halfSize *= 2) { complex&ltdouble&gt phaseShiftStep = polar(10,-PI/halfSize) complex&ltdouble&gt currentPhaseShift(1,0)

for(int fftStep = 0 fftStep &lt halfSize fftStep++) {

for(int i=fftStep i &lt length i += 2*halfSize) {

complex&ltdouble&gt t = currentPhaseShift *

samples[i+halfSize]

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

samples[i] += t

}

currentPhaseShift *= phaseShiftStep

}

}

}

/ / Використовуючи ряд простих властивостей комплексних чисел,

/ / Ми можемо обчислити зворотне БПФ шляхом транспонування

/ / Відліків до і після обчислення прямого ШПФ

/ / Заодно тут вони домножаются на 1 / N

void InverseFft(complex&ltdouble&gt *samples, int length ) {

for(int i=0 i&ltlength i++)

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

for(int i = 0 i&ltlength i ++ )

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

static_cast&ltdouble&gt(length)

}

Перегрупування відліків складне завдання Для обчислення кожного БПФ необхідно спочатку взяти всі парні елементи, а потім всі непарні Цей поділ потім повторюється для кожної половини

Складність полягає в аналізі довічного значення індексу масиву для кожного елемента Згадаймо, що парне число містить 0 в молодшому бите і що парні елементи будуть розміщуватися в першому половині масиву, для якої індекс містить 0 в старшому біті Перший крок перестановки полягає у переміщенні молодшого біта індексу масиву в старший біт

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

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

Лістинг 246 Перестановка відліків для швидкої реалізації ШПФ

static void Rearrange(complex&ltdouble&gt *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&ltlength limit &lt&lt= 1, bit&gt&gt=1 )

for(int i=0i&ltlimiti++)

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

/ / Розміщуємо 0 в ті елементи, які залишаються

/ / Без змін

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

/ / Щоб кожен обмін відбувався тільки один раз

for(int i=0 i&ltlength i++) {

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

else rearrange[ rearrange[i] ] = 0

}

rearrangeSize = length

}

/ / Використовуємо таблицю перестановки для перестановки елементів

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

complex &ltdouble&gt t

for (int i=0 i&ltlength 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-точкове ШПФ І, нарешті, замість арифметики з плаваючою точкою ви можете використовувати целочисленную арифметику або арифметику з фіксованою точкою

Експерименти з БПФ

Наступні дві невеликі програми призначені для дослідження ШПФ Кожна читає текстовий файл, що містить числа, здійснює відповідне перетворення і генерує вихідний файл аналогічного формату

Лістинг 247 Програма fftmaincpp

#include &quotffth&quot

#include &ltiostream&gt

#include &ltiomanip&gt

#include &ltcstdlib&gt

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

int length = 8

if (argc &gt 1) length = atoi(argv[1])

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

/ / Довжина повинна бути ступенем 2

cerr &lt&lt &quotLength must be a power of two\n"

exit(1)

}

complex&ltdouble&gt *input = new complex&ltdouble&gt[length]

for(int i=0 (i &lt length)&amp&amp(cineof()) i++)

cin &gt&gt input[i]

complex&ltdouble&gt *output = new complex&ltdouble&gt[length]

for(int i=0i&ltlengthi++) output[i] = input[i] ForwardFft(output,length)

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

#if 0

complex&ltdouble&gt *compare = new complex&ltdouble&gt[length]

ForwardDft(input,length,compare)

double maxerr=0 int maxindex=0

for(int i=0i&ltlengthi++) {

/ / Обчислюємо відносну помилку

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

abs(compare[i])

if(error &gt maxerr) {

maxindex=i maxerr=error

}

}

/ / Максимальна відносна помилка

cerr &lt&lt &quotMaximum relative error:&quot &lt&lt maxerr &lt&lt &quot\n"

#endif

for(int i=0i&ltlengthi++)

cout &lt&lt setprecision(20) &lt&lt output[i] &lt&lt &quot\n"

return 0

}

Ці дві програми використовують стандартну бібліотеку C + + для читання комплексних чисел із стандартного потоку введення Стандартна бібліотека C + + працює з комплексними числами двох форматів Якщо перший символ кругла дужка, очікується введення укладеної в дужки пари чисел, що представляють дійсну і уявну частини B іншому випадку зчитується одне число з плаваючою точкою, що представляє дійсну частину, а уявна частина покладається рівною нулю Наприклад, і (1, 0), і 1,0 будуть сприйняті як комплексне число 1 + 0i

Лістинг 248 Програма ifftmaincpp

#include &quotffth&quot

#include &ltiostream&gt

#include &ltiomanip&gt

#include &ltcstdlib&gt

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

int length = 8

if (argc &gt 1) length = atoi(argv[1])

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

/ / Довжина повинна бути ступенем числа 2

cerr &lt&lt &quotLength must be a power of two\n"

exit(1)

}

complex&ltdouble&gt *a = new complex&ltdouble&gt[length]

for(int i=0 (i &lt length)&amp&amp(cineof()) i++)

cin &gt&gt a[i] InverseFft(a,length) for(int i=0i&ltlengthi++)

cout &lt&lt setprecision(20) &lt&lt a[i] &lt&lt &quot\n"

return 0

}

Джерело: Кінтцель Т Керівництво програміста по роботі зі звуком = A Programmers Guide to Sound: Пер з англ М: ДМК Пресс, 2000 432 с, іл (Серія «Для програмістів»)

Схожі статті:


Сподобалася стаття? Ви можете залишити відгук або підписатися на RSS , щоб автоматично отримувати інформацію про нові статтях.

Коментарів поки що немає.

Ваш отзыв

Поділ на параграфи відбувається автоматично, адреса електронної пошти ніколи не буде опублікований, допустимий HTML: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

*

*