|
|
  |
FFT на асм для ARM7TDMI (AT91SAM7xx) |
|
|
|
Nov 21 2012, 15:52
|

山伏
    
Группа: Свой
Сообщений: 1 827
Регистрация: 3-08-06
Из: Kyyiv
Пользователь №: 19 294

|
Цитата(hd44780 @ Nov 21 2012, 17:01)  DRUID3, а чем отличаются fn_aT_ditNbrRadix2FFT_int и fn_aT_ditNbrRadix2ReFFT_int и какую мне лучше использовать? ...ну у Вас же прямое преобразование потому первую. Отличаются они знаком комплексной экспоненты и масштабирующим коэффициентом 1/N. Я его внес в каскады в виде битовых сдвигов прямого преобразования - так спектр в масштабе сигнала получается. Т.е. если синус амплитудой 16384, то и палка от него будет такой же. Кстати Вы пишете, что картинка Вас не радует на винамп не похожа. Так в винампе она в логарифмическом масштабе. А в прямом там действительно только основные гармоники будут. P.S.: ...там мое тоже не ахти по оптимальности. Нет битреверса потому использованы 2-а дополнительных массива. Но для того чтобы "сэкономить место" используется факт, что синус это косинус + pi/2, а это дополнительные сложения. Плюс там еще "прыжки по таблице тригонометрии" - это чтобы "маленькие" FFT использовали таблицу для "больших" - и это жрет ресурс. Можно сделать битреверс и "вернуться к классике" - это я просто так экспериментировал. В сети полно FFT степени 2 с битреверсом, а такой формы нет нигде что-то  .
--------------------
Нас помнят пока мы мешаем другим... //-------------------------------------------------------- Хороший блатной - мертвый... //-------------------------------------------------------- Нет старик, это те дроиды которых я ищу...
|
|
|
|
|
Nov 21 2012, 17:59
|

Профессионал
    
Группа: Свой
Сообщений: 1 202
Регистрация: 26-08-05
Из: Донецк, ДНР
Пользователь №: 7 980

|
Цитата(DRUID3 @ Nov 21 2012, 17:52)  ну у Вас же прямое преобразование потому первую. Спасибо. Я ещё смотрю. Будут вопросы, спрошу. Цитата(DRUID3 @ Nov 21 2012, 17:52)  Кстати Вы пишете, что картинка Вас не радует на винамп не похожа. Так в винампе она в логарифмическом масштабе. А в прямом там действительно только основные гармоники будут. Да, хотелось бы поближе к винампу  ... Там что, надо еще по каждой полоске 10log10 посчитать и взять модуль, тк. dB всегда отрицательные? А отношение к чему? Написал туда сейчас тупо Spectrum [i] = (int)(10.0 * log10 ( Re2 + Im2 ) ); вообще ничего не рисует Цитата(DRUID3 @ Nov 21 2012, 17:52)  В сети полно FFT степени 2 с битреверсом, а такой формы нет нигде что-то  . Да отож ...
--------------------
Чтобы возить такого пассажира, необходим лимузин другого класса. (с) Мария Эдуарда
|
|
|
|
|
Nov 21 2012, 19:13
|

Профессионал
    
Группа: Свой
Сообщений: 1 202
Регистрация: 26-08-05
Из: Донецк, ДНР
Пользователь №: 7 980

|
Вот сейчас посмотрел Чана - http://elm-chan.org/works/akilcd/report_e.html , логарифмов там нету, а картинка с виду нормальная ... Неужели квадратный корень так влияет  ? Вот проверил, корень добавил - Spectrum [i] = (int) sqrt ( (float)Re2 + (float)Im2 ); Вообще ничего не рисует .... Наверное всё-таки мой алгоритм кривой
--------------------
Чтобы возить такого пассажира, необходим лимузин другого класса. (с) Мария Эдуарда
|
|
|
|
|
Nov 22 2012, 03:32
|

山伏
    
Группа: Свой
Сообщений: 1 827
Регистрация: 3-08-06
Из: Kyyiv
Пользователь №: 19 294

|
У Чана логарифмов нету. У него чистое FFT. Корень квадратный навряд-ли так влияет. Может что с FFT. А Вы ударьте по своему FFT дельта-импульсом. Что будет на выходе? Цитата(hd44780 @ Nov 21 2012, 19:59)  Да отож ... ...мы не идем проторенным путем...  Код /*=======================================================================*/ int fn_aT_ditNbrRadix2FFT_int ( const unsigned int const cuic_N, int ai_inRe[], int ai_inIm[], int ai_buffRe[], int ai_buffIm[], int **pp_outRe, int **pp_outIm, int **pp_inRe, int **pp_inIm, const int const acc_trigT[], const unsigned int cui_jumpT ) ...там так cui_jumpT - это прыжки по таблице тригонометрии. если у Вас в проекте все FFT одного порядка смело выбрасывайте все математические действия с этой переменной и ее саму. Указатели на входной/выходной массивы потому как в зависимости от порядка FFT(количества каскадов) результат будет либо во входных массивах либо в буфферных - такой алгоритм. Во внешнем алгоритме(realFFT за 2-а прохода) обявите просто еще и эти массивы. Код ui_TargS = cui_jumpT*ui_cWcount*ui_distance; /* шаг для sin(t) */ ui_TargC = ui_TargS + cui_jumpT*(cuic_N>>2); /* шаг для cos(t) */ если синус и косинус просто без изврата поместить в 2-е раздельные таблицы, то вычисление индекса массива тригонометрии упростится до: Код ui_TargT = ui_cWcount*ui_distance; /* шаг для sin(t)/cos(t) */ Вот где-то так...
--------------------
Нас помнят пока мы мешаем другим... //-------------------------------------------------------- Хороший блатной - мертвый... //-------------------------------------------------------- Нет старик, это те дроиды которых я ищу...
|
|
|
|
|
Nov 22 2012, 09:52
|

Профессионал
    
Группа: Свой
Сообщений: 1 202
Регистрация: 26-08-05
Из: Донецк, ДНР
Пользователь №: 7 980

|
DRUID3 Я взял пока функцию, работающую с плавающей точкой. Вот, что у меня получилось: Код #define FFT_N 512
void DoFFT (byte *pInData, byte iWindow) { float a_inReal1[FFT_N], a_inReal2[FFT_N]; float a_outRe1[FFT_N], a_outIm1[FFT_N]; float a_outRe2[FFT_N], a_outIm2[FFT_N]; int i; float temp;
// 1. Заполнение входных массивов for (i=0; i<FFT_N; i++) { temp = (float)pInData[i]; a_inReal1[i] = a_inReal2[i] = temp; } // for // 2. FFT fn_a_2RealFFT ( FFT_N, a_inReal1, a_inReal2, a_outRe1, a_outIm1, a_outRe2, a_outIm2, a_TTransf, STEP(FFT_N, 1536) );
// 3. Расчёт спектра for ( i = 0; i < FFT_N; i ++ ) { float Re2, Im2; Re2 = a_outRe1 [ i ]; Im2 = a_outIm1 [ i ];
// Squareing Spectrum [i] = (int) sqrt ( Re2 + Im2 ); } // for } // DoFFT Я сомневаюсь, правильно ли я написал первый аргумент макроса STEP. У Вас там 8 == размеру массивов, т.е. FFT_N. Оконные функции пока выкинул вообще (чтоб не путались пока под ногами  ), потом добавлю. Теперь мне надо определить спектр. Но я не понял, использовать пару a_outRe1/a_outIm1 или a_outRe2/a_outIm2? Как я понял из алгоритма и описания в конце, что, т.к. я возвожу всё добро в квадрат, то всё равно. Я прав? Спасибо.
--------------------
Чтобы возить такого пассажира, необходим лимузин другого класса. (с) Мария Эдуарда
|
|
|
|
|
Nov 22 2012, 10:39
|

山伏
    
Группа: Свой
Сообщений: 1 827
Регистрация: 3-08-06
Из: Kyyiv
Пользователь №: 19 294

|
Цитата(hd44780 @ Nov 22 2012, 11:52)  DRUID3 Я взял пока функцию, работающую с плавающей точкой. Вот, что у меня получилось: Код // 1. Заполнение входных массивов for (i=0; i<FFT_N; i++) { temp = (float)pInData[i]; a_inReal1[i] = a_inReal2[i] = temp; } // for Неверно. a_inReal1[] это один буфер, а a_inReal2[] следующий за ним... Иначе в чем ускорение то? Секрет в том что на основе БПФ на 512 комплексных отсчета можно посчитать 2-а БПФ на столько же отсчетов одновременно!!! Цитата(hd44780 @ Nov 22 2012, 11:52)  Я сомневаюсь, правильно ли я написал первый аргумент макроса STEP. У Вас там 8 == размеру массивов, т.е. FFT_N. ...кажется верно... Цитата(hd44780 @ Nov 22 2012, 11:52)  Теперь мне надо определить спектр. Но я не понял, использовать пару a_outRe1/a_outIm1 или a_outRe2/a_outIm2? Как я понял из алгоритма и описания в конце, что, т.к. я возвожу всё добро в квадрат, то всё равно. Я прав?
Спасибо. Это 2-а абсолютно независимых спектра от первого и от второго буфера соответственно.
--------------------
Нас помнят пока мы мешаем другим... //-------------------------------------------------------- Хороший блатной - мертвый... //-------------------------------------------------------- Нет старик, это те дроиды которых я ищу...
|
|
|
|
|
Nov 22 2012, 11:39
|

Профессионал
    
Группа: Свой
Сообщений: 1 202
Регистрация: 26-08-05
Из: Донецк, ДНР
Пользователь №: 7 980

|
Цитата(DRUID3 @ Nov 22 2012, 12:39)  Неверно. a_inReal1[] это один буфер, а a_inReal2[] следующий за ним... Иначе в чем ускорение то? Секрет в том что на основе БПФ на 512 комплексных отсчета можно посчитать 2-а БПФ на столько же отсчетов одновременно!!! Т.е. я могу представить свою 512-байтовую выборку как 2 256-байтовых, так что ли? Типа такого: Код float a_inReal1[FFT_N/2], a_inReal2[FFT_N/2]; float a_outRe1[FFT_N/2], a_outIm1[FFT_N/2]; float a_outRe2[FFT_N/2], a_outIm2[FFT_N/2]; int i;
// 1. Заполнение входных массивов for (i=0; i<FFT_N/2; i++) { a_inReal1[i] = (float)pInData[i]; a_inReal2[i] = (float)pInData[i+FFT_N/2]; } // for Цитата(DRUID3 @ Nov 22 2012, 12:39)  Это 2-а абсолютно независимых спектра от первого и от второго буфера соответственно. В свете первого замечания понятно.
--------------------
Чтобы возить такого пассажира, необходим лимузин другого класса. (с) Мария Эдуарда
|
|
|
|
|
Nov 22 2012, 12:35
|

Профессионал
    
Группа: Свой
Сообщений: 1 202
Регистрация: 26-08-05
Из: Донецк, ДНР
Пользователь №: 7 980

|
Понял. Тут ещё один вопрос.... Для входного окна размером FFT_N обычно я получаю спектр размером FFT_N/2, т.к. вторая половина симметрична первой относительно середины массива. По Вашему алгоритму также? Если я просчитаю спектр: Код int Spectrum [FFT_N];
for ( i = 0; i < FFT_N/2; i ++ ) { float Re2, Im2; Re2 = a_outRe1 [ i ]; Im2 = a_outIm1 [ i ]; Spectrum [i] = (int) sqrt ( Re2 + Im2 ); Re2 = a_outRe1 [ i + FFT_N/2 ]; Im2 = a_outIm1 [ i + FFT_N/2 ]; Spectrum [i + FFT_N/2] = (int) sqrt ( Re2 + Im2 ); } // for То я должен буду, по сути взять 1-ю и 3-ю четверти маcсива Spectrum - последовательно рисовать элементы 0..FFT_N/4-1 и FFT/2 .. FFT/2+FFT_N/4-1? Или ещё как-то?
--------------------
Чтобы возить такого пассажира, необходим лимузин другого класса. (с) Мария Эдуарда
|
|
|
|
|
Nov 22 2012, 14:45
|

Профессионал
    
Группа: Свой
Сообщений: 1 202
Регистрация: 26-08-05
Из: Донецк, ДНР
Пользователь №: 7 980

|
Вроде понятно. В итоге получил такое: Код float a_inReal1[FFT_N/2], a_inReal2[FFT_N/2]; float a_outRe1[FFT_N/2], a_outIm1[FFT_N/2]; float a_outRe2[FFT_N/2], a_outIm2[FFT_N/2]; int i;
// 1. Заполнение входных массивов for (i=0; i<FFT_N/2; i++) { a_inReal1[i] = (float)pInData[i]; a_inReal2[i] = (float)pInData[i+FFT_N/2]; } // for // 2. FFT fn_a_2RealFFT ( FFT_N/2, a_inReal1, a_inReal2, a_outRe1, a_outIm1, a_outRe2, a_outIm2, a_TTransf, STEP(FFT_N/2, 1536) );
// 3. Расчёт спектра for ( i = 0; i < FFT_N/2; i ++ ) { float Re2, Im2; Re2 = a_outRe1 [ i ]; Im2 = a_outIm1 [ i ]; Spectrum [i] = (int) sqrt ( Re2 + Im2 ); Re2 = a_outRe1 [ i + FFT_N/2 ]; Im2 = a_outIm1 [ i + FFT_N/2 ]; Spectrum [i + FFT_N/2] = (int) sqrt ( Re2 + Im2 ); } // for Дома на железяке проверю.
--------------------
Чтобы возить такого пассажира, необходим лимузин другого класса. (с) Мария Эдуарда
|
|
|
|
|
Nov 22 2012, 19:29
|

Профессионал
    
Группа: Свой
Сообщений: 1 202
Регистрация: 26-08-05
Из: Донецк, ДНР
Пользователь №: 7 980

|
Исправил косяк Код // 3. Расчёт спектра for ( i = 0; i < FFT_N/2; i ++ ) { float Re2, Im2; Re2 = a_outRe1 [ i ]; Im2 = a_outIm1 [ i ]; Spectrum [i] = (int) sqrt ( Re2*Re2 + Im2*Im2 ); Re2 = a_outRe1 [ FFT_N/2 + i ]; Im2 = a_outIm1 [ FFT_N/2 + i ]; Spectrum [ FFT_N/2 + i ] = (int) sqrt ( Re2*Re2 + Im2*Im2 ); } // for Заработал. Результат воде ничего, но тормозит .. Единственное - рисуется только левая половина. Наверное с этими частями намутил  Завтра целочисленный вариант попробую прикрутить.
--------------------
Чтобы возить такого пассажира, необходим лимузин другого класса. (с) Мария Эдуарда
|
|
|
|
|
Nov 23 2012, 06:58
|

Профессионал
    
Группа: Свой
Сообщений: 1 202
Регистрация: 26-08-05
Из: Донецк, ДНР
Пользователь №: 7 980

|
Цитата(DRUID3 @ Nov 23 2012, 07:19)  лучше ударить дельта-импульсом и посмотреть на выход - RealFFT это тоже ведь линейная система. Я не знаю, что такое "дельта-импульс"  . И где его взять ... В гугле поискал - внятного ничего не нашёл. Могу придумать только генерить частоты в каком-нибудь CoolEdit-е и ими тестировать. Цитата(DRUID3 @ Nov 23 2012, 07:19)  1) а a_outRe2[], a_outIm2[] куда дели? 2) a_outIm1 [ FFT_N/2 + i ] - это что такое вообще? Зачем, откуда? У Вас уже есть готовые полуспектры, забудьте Вы о математическом формализме и этой сопряженной симметрии. Да то уже в полусне писал  , вот и получилась фигня ... И работает там кстати только левая половина. Наверно поэтому. Изначально я так собирался написать: Код // 3. Расчёт спектра for ( i = 0; i < FFT_N/2; i ++ ) { float Re2, Im2; Re2 = a_outRe1 [ i ]; Im2 = a_outIm1 [ i ]; Spectrum [i] = (int) sqrt ( Re2*Re2 + Im2*Im2 ); Re2 = a_outRe2 [ i ]; Im2 = a_outIm2 [ i ]; Spectrum [ FFT_N/2 + i ] = (int) sqrt ( Re2*Re2 + Im2*Im2 ); } // for Про "готовые полуспектры" не понял  . Вы хотите сказать, что после Вашего алгоритма не нужно эти корни брать? Цитата(DRUID3 @ Nov 23 2012, 07:19)  Вы же хотите так как в Winamp? Совершенно верно
--------------------
Чтобы возить такого пассажира, необходим лимузин другого класса. (с) Мария Эдуарда
|
|
|
|
|
  |
1 чел. читают эту тему (гостей: 1, скрытых пользователей: 0)
Пользователей: 0
|
|
|