Быстрое Преобразование Фурье, java

Решил я вплотную подступиться к Цифровой Обработке Сигналов.
Достаточно долгое время не мог понять, как реализовать алгоритм БПФ. Математика была понятна, но просмотрев готовые решения, понять как они работают не удавалось. На какое-то время необходимость в этом отпала, а потом снова появилась.
Ко второму заходу опыта программирования стало уже немного большое, и, наверно, поэтому мне удалось самому написать программу, выполняющую БПФ на java (код без особого труда можно портировать на C++ или C).
Я не буду выводить основные формулы и доказывать, что они верны — с этим справляются на многих сайтах и во многих книжках (по крайней мере лично у меня это не вызвало трудностей).
Основные формулы1

Я просто расскажу о том, как пользуясь этими формулами я писал саму программу, с какими трудностями столкнулся и как их преодолел.
Статья, по сути предназначена для тех кто:
1. Прочитал и понял математику, стоящую за БПФ.
2. Понял алгоритм бит-реверса.
3. Понял почему граф преобразования выглядит именно так

Первое, что я сделал — поставил основную задачу «Разработать программу быстрого преобразования Фурье на java». Затем приступил ко второму этапу: сбору информации. Прочитав доказательства и поняв математику, стоявшую за алгоритмом, я понял, следующие вещи:
1. Алгоритм можно спроектировать пользуясь парадигмой Divide and Conquer.
2. В данной парадигме задачу можно решить рекурсивно или итеративно.
Рекурсивный подход, как мне кажется, в какой-то мере проще в написании и понимании, но приведет к некоторым потерям в скорости из-за вложенных вызовов методов (функций).

Напомню или расскажу тем, кто не в курсе, что парадигма Divide and Conquer (Разделяй и Властвуй) сводится к трем общим этапам.
1. Этап разделения задачи на более простые задачи меньшего размера.
2. Этап решения всех этих задач.
3. Этап объединения результатов всех меньших задач в общий результат.

БПФ же можно разбить на три основных этапа:

1. Перестановка входных отсчетов (значений в массиве). Почему это нужно много рассказывать не буду, но в общем, как я понял, это позволяет сэкономить память, т.к. в таком случае не требуется массив для хранения промежуточных результатов и в целом упрощает написание программы.
2. БПФ основано на том факте, что если N = 2^n (т.е. 1,2,4,...,1024...2^n), то N-точечное ДПФ (т.е. количество отсчетов исходных данных равно числу n) можно разбить на два ДПФ размером вдвое меньше (2^(n-1)). В этом и заключается второй этап парадигмы. Мы разбиваем ДПФ, пока не приходим к базовому случаю — ДПФ от двух отсчетов, которое сводится к сложению и вычитанию двух элементов массива.
3. Объединение результатов каждого из разбиений.

Дальше можно собственно приступать к разработке алгоритма предварительно разбив его на эти отдельные задачи.
1. Разработать метод перестановки входных отсчетов с помощью алгоритма бит-реверса.
Если непонятно почему бит-реверс, пишите в коммантариях.
Собственно алгоритм бит-реверса

i — целое число, в котором надо переставить t-битов.

private static int reverse(int i, int t) {
		// TODO Auto-generated method stub
	    int Shift = t - 1;
	    int LowMask = 1;
	    int HighMask = 1 << Shift;
	    int R;
	    for(R = 0; Shift >= 0; LowMask <<= 1, HighMask >>= 1, Shift -= 2)
	        R |= ((i & LowMask) << Shift) | ((i & HighMask) >> Shift);
	    return R;
	}


А это, собственно, перестановка значений.

for(int i = 1; i < Nmax-1; i++)
		{
		    int J = reverse(i,powerOfTwo);           // reverse ����������� ���� � I � �������� �������
		    if (i >= J)                     // ���������� ��� ��������������
		    	continue;
		    double S = array[i]; array[i] = array[J]; array[J] = S;  // ������������ ��������� xI � xJ
		}


Т.е. здесь мы берем число i, переставляя в нем биты, получаем число J. i и J — это индексы двух элементов массива, которые нужно поменять местами. Ну и меняем их местами. Повторяем эту процедуру пока не пройдемся по всем элеметам массива, т.е. пока i не станет больше или равно Nmax-1.

А дальше начинается самое интересное.

*************************************************************
Статья пока недописана, но публикую, т.к. интересуют комментарии. О чем написать подробнее, что может вызывать вопросы и пр. Любые мнения и помощь приветствуются!
  • 0
  • 12 ноября 2012, 15:43
  • Boris_B

Комментарии (33)

RSS свернуть / развернуть
Предыдущий дубль удали, пока его не прокомментировали.
0
  • avatar
  • Vga
  • 12 ноября 2012, 16:04
нельзя, он в коллективных
0
Его можно утянуть в личку и удалить уже из лички.
0
а, ну да:)
0
Спасибо за подсказку, тут какие-то проблемы были (502 Bad Gateway), поэтому дубль получился.
0
и КАТ поставь.
-1
А может нам стоит написать кат-бота? :))
+1
живых мало? ;)

проще ди допинать, чтобы прилепил автокат. благо есть такой плаг к лайвстриту.
0
Хм, а можно посмотреть. Не знал
0
да я когда-то уже вспоминал о нем. причем автокат умный. если руками не поставить, ставит принудительно после определенного количества символов.
0
Кому-то вообще интересна тема?
Потому что если это никому не нужно, то мне тем более и пост тогда просто удалю =)
0
Этот уже не удалишь) Правов не хватит (да и на тот, возможно, тоже). Так что дописывай. Ща пока особо комментировать нечего.
0
пиши, автор, штука нужная, ща сам как раз собираюсь вкуривать и писать, ибо понадобилось.
можно максиматьно подробно, можно неочень.
и желательно чтобы не завязанно на длинну, ато везде примеры на 4 или 8 точек, а как больше сделать я не пойму никак.
0
  • avatar
  • kest
  • 12 ноября 2012, 18:35
интересует очень
0
Спасибо за комментарии и проявленный интерес. Раз интерес есть, тогда через пару дней можете ждать полную версию. =)
0
Если не секрет, зачем БПФ на яве?
0
БПФ основано на том факте, что если n = 2^n

Формулу поправь.
0
Спасибо за поправку. пока нет времени писать. При установки убунты почему-то полетела винда и теперь я разбираюсь с убунтой. Четыре часа разбирался почему не работает блютуз и как его включить (у меня мышь беспроводная). В общем со временем напряги.

Любые части по пожеланиям товарищей распишу подробнее, если будет необходимо (бит-реверс, почему граф именно такой и пр.)

Позже напишу собственно алгоритм объединения результатов.

БПФ на джаве… БПФ мне нужно не на джаве в итоге, а на DSP реализовывать (надеюсь, с ПЛИС связываться не придется, т.к. знаний в этой области маловато). А на джаве хочу сделать математическую модель тракта цифровой обработки. Мне с неё проще код перевести будет, чем с матлаба. А матлаб вообще плохо знаю.
Кстати, знает ли кто тут хорошие библиотеки ЦОС для C или C++?
0


Пока можете посмотреть как выглядит итог работы программы.
0
Спасибо что затронули интересную тему. Ждем продолжения.
Хорошо алгоритм FFT описан в книжках по ЦОС и DSP-процессорам.
В некоторых DSP даже есть аппаратный модуль выбора адреса для FFT (Bit-reversed addressing). Так что реализация FFT на Java будет сложнее чем для DSP.
Качество и время работы программы очень сильно зависит от реализации. Разрабатывать FFT самому имеет смысл только для обучения, в большинстве случаев, своя реализация будет менее оптимальна, чем написанная и отлаженная кем то до Вас.
0
чем не устроили готовые библиотеки?
0
1. Пишу, конечно, в целях обучения.
Если обрабатывать данные не налету, то можно это делать и на джаве, как мне кажется. Не думаю, что она будет сильно медленнее при том же алгоритме. BigO-то тоже самое будет. Например, можно увеличить скорость бит-реверса, если делать его по маскам, но все алгоритмы бит-реверса будут линейно зависеть от n, если правильно помню.

Если библиотеки куплены или есть уже отлаженные, почему бы не воспользоваться, кто бы спорил =)

Только лично я свободной библиотеки на Си для DSP не нашел. Вроде, каждая компания производитель DSP процессоров выпускает свои библиотеки (оно и логично, потому что там будет много фишек конкретной архитектуры использоваться), но, к примеру, под тот комплект отладочный, что заказали у нас на работе, они платные. У любителя на них, наверно, просто не хватит денег.

Был бы признателен, если бы просветили, где можно взять библиотеки и на что именно от них можно рассчитывать?
0
лично я свободной библиотеки на Си для DSP не нашел
www.fftw.org/, spuc.sourceforge.net/, www.gnu.org/software/gsl/, www.vsipl.org/
0
Бит реверс:

unsigned rev (unsigned x)
{
 x = (x & 0x55555555) << 1 | (x >> 1) & 0x55555555;
 x = (x & 0x33333333) << 2 | (x >> 2) & 0x33333333;
 x = (x & 0x0F0F0F0F) << 4 | (x >> 4) & 0x0F0F0F0F;
 x = (x <<24) | ((x & 0xFF00) <<8) |
     ((x >> 8) & 0xFF00) | (x >> 24);
 return x;
}

«Алгоритмические трюки для программистов»
0
Байты тоже можно реверсить половинками, как биты в них:
unsigned rev (unsigned x)
{
 x = (x & 0x55555555) << 1 | (x >> 1) & 0x55555555;
 x = (x & 0x33333333) << 2 | (x >> 2) & 0x33333333;
 x = (x & 0x0F0F0F0F) << 4 | (x >> 4) & 0x0F0F0F0F;
 x = (x & 0x00FF00FF) << 8 | (x >> 8) & 0x00FF00FF;
 x = (x & 0x0000FFFF) << 16 | (x >> 16) & 0x0000FFFF;
 return x;
}
0
Стоит правда отметить, что на этих ваших МК сдвиги неэффективны и может оказаться быстрее сдвигать по одному биту источника в бит переноса и затем задвигать его в приемник. Правда, это только на асме.
0
как там дела со статьей?;)

хочу отметить что реализация БПФ с википедии очень даже неплохо переносится куда либо.
но таки разобраться что и как было бы интересно
0
  • avatar
  • kest
  • 01 января 2013, 18:18
непонятно, при чём тут Ява? Из статьи пока непонятно, при чём тут БПФ, пока только видно бит-реверс. Я когда-то писал статью в журнал про БПФ на SSE, тогда ещё первом, там была описана пара приёмов, которые оптимизируются именно на SSE. А тут какое-то сферическое БПФ в вакууме. Т.е. я не обидеть хочу, а намекнуть, что когда пишут про БПФ относительно железа, то пишут, что это за железо, и что оно может хорошего для БПФ.
Кстати, советую автору обратить внимание на FFTW. Если ему просто прикрутить к произвольному процу — самое оно.
0
Про бит-реверс всё понятно в двух словах: каждая стадия на графе БПФ переписывает данные задом-наперёд относительно их номера в битах. В итоге все данные оказываются на местах, номера которых в битах являются зеркальным отражением входных. Кстати, отсюда следует, что можно переставлять данные на входе, а можно на выходе.
0
Извините, не сдержался…
Особенно убила фраза:
2. Понял алгоритм бит-реверса.

Если вы действительно всё поняли, надо было написать
2. Понял, откуда в БПФ берётся бит-реверс.
0
Как говорится, или крестик снимите, или трусы наденьте.
Я про Divide and conquer.
Надо писать разделяй и властвуй, либо Divide et impera.
0
Кстати, бит-реверс я делал через LFSR
Один сдвигался вправо, другой — влево. Один сдвиг и один XOR. Т.о. было 2 индекса, скажем, j,i, они были всё время зеркальными по построению. Бит реверс на выходе тогда делался t=X[i]; X[i]=X[j]; X[j]=t;
Чтоб ускорить, можно было генерить индексы пачками:
например, i, j покрывают всё от 0 до 7, а у нас 16 сэмплов, делаем так:

j1 = (j<<1) | 0;
j2 = (j<<1) | 1;
i1 = (j<<1) | (0<<4);
i2 = (j<<1) | (1<<4);
<code>
или, например, i, j покрывают всё от 0 до 7, а у нас 32 сэмпла, делаем так:<br />
<code>
j1 = (j<<2) | 0;
j2 = (j<<2) | 1;
j3 = (j<<2) | 2;
j4 = (j<<2) | 3;
i1 = (j<<2) | (0<<5);
i2 = (j<<2) | (2<<5);
i3 = (j<<2) | (1<<5);
i4 = (j<<2) | (3<<5);
<code>
думаю, смысл понятен.
0
внимательный читатель заметит косяк — конешно, LSFR покрывает не от 0 до 7, а от 1 до 7. Но 0й элемент переставлять не нужно, так что всё нормально.
0
Только зарегистрированные и авторизованные пользователи могут оставлять комментарии.