| |
|
  |
 FFT на асм для ARM7TDMI (AT91SAM7xx) |
|
|
|
|
 Nov 23 2012, 08:27
|
山伏
Группа: Свой
Сообщений: 1 827
Регистрация: 3-08-06
Из: Kyyiv
Пользователь №: 19 294
|
Цитата(hd44780 @ Nov 23 2012, 08:58)  Я не знаю, что такое "дельта-импульс"  . И где его взять ... Вы очень много в жизни потеряли. А он заслуживает траты на него времени. Ну хотя-бы потому, что ударь Вы им по любой линейной системе и Вы сразу все узнаете о ней - получите ее АЧХ/ФЧХ или т.н."передаточную функцию" если говорить о временнОй области. Цитата(hd44780 @ Nov 23 2012, 08:58) В гугле поискал - внятного ничего не нашёл. ...0, 0, 0, 0, 0, 1, 0, 0, 0... - это она(функция)... ну или он(импульс)...  Цитата(hd44780 @ Nov 23 2012, 08:58) Про "готовые полуспектры" не понял . Вы хотите сказать, что после Вашего алгоритма не нужно эти корни брать? Нет, уже все хорошо...
--------------------
Нас помнят пока мы мешаем другим...
//--------------------------------------------------------
Хороший блатной - мертвый...
//--------------------------------------------------------
Нет старик, это те дроиды которых я ищу...
|
|
|
|
|
|
|
|
Nov 24 2012, 07:58
|
Профессионал
Группа: Свой
Сообщений: 1 202
Регистрация: 26-08-05
Из: Донецк, ДНР
Пользователь №: 7 980
|
Отписаться вечером не смог .. Вчера проверил - на плавающей точке работает, но тормозит - 1-2 кадра в секунду ... Начал смотреть целочисленный вариант FFT: Код int fn_aT_ditNbrRadix2FFT_int
(
const unsigned int cuic_N,
int ai_inRe[], // i input
int ai_inIm[], // q input
int ai_buffRe[], // i buff
int ai_buffIm[], // q buff
int **pp_outRe, // i output
int **pp_outIm, // q output
int **pp_inRe, // i input
int **pp_inIm, // q input
const int acc_trigT[],
const unsigned int cui_jumpT
) Нашёл пример её вызова: Код fn_aT_ditNbrRadix2FFT_int
(
ui_N,
i, /* i input */
q, /* q input */
bi, /* i buff */
bq, /* q buff */
&p1, /* i output */
&p2, /* q output */
&p3, /* i input */
&p4, /* q input */
cai_TTransf,
STEP(ui_N, dfi_3N4)
); i, q - вход вроде. А что такое - bi, bq? 2-й вход, по аналогии с плавающим вариантом? А p1, p2, p3, p4 - результаты по 2-м входным буферам? Так что-ли?
--------------------
Чтобы возить такого пассажира, необходим лимузин другого класса.
(с) Мария Эдуарда
|
|
|
|
|
|
|
|
Nov 24 2012, 16:24
|
Профессионал
Группа: Свой
Сообщений: 1 202
Регистрация: 26-08-05
Из: Донецк, ДНР
Пользователь №: 7 980
|
DRUID3, чёт я не понял . Написал так: Код int iData [FFT_N], qData [FFT_N];
int bi [FFT_N], bq [FFT_N];
...........................
void DoFFT (byte *pInData, byte iWindow)
{
int i;
int *p1 = 0, *p2 = 0, *p3 = 0, *p4 = 0;
// 1. Заполнение входных массивов
for ( i = 0; i < FFT_N; i ++ )
{
iData[i] = (int)pInData[i];
qData[i] = 0;
} // for
// 2. FFT
fn_aT_ditNbrRadix2FFT_int
(
FFT_N,
iData, // i input
qData, // q input
bi, // i buff
bq, // q buff
&p1, // i output
&p2, // q output
&p3, // i input
&p4, // q input
cai_TTransf, STEP(FFT_N, dfi_3N4)
);
// 3. Расчёт спектра
for ( i = 0; i < FFT_N; i ++ )
{
int Re2, Im2;
Re2 = p1 [ i ];
Im2 = p2 [ i ];
Spectrum [i] = (int) sqrt ( (float)Re2*(float)Re2 + (float)Im2*(float)Im2 );
} // for
} // DoFFT Работает это быстро  , вот только рисует 1-2 полоски где-то вначале .... Я в функции fn_aT_ditNbrRadix2FFT_int использование &p3 и &p4 закрыл, им там что-то присваивается, далее не используется. Потом вообще думаю их убрать ... Хотя и как было не работало.
--------------------
Чтобы возить такого пассажира, необходим лимузин другого класса.
(с) Мария Эдуарда
|
|
|
|
|
|
|
|
 Nov 24 2012, 17:24
|
山伏
Группа: Свой
Сообщений: 1 827
Регистрация: 3-08-06
Из: Kyyiv
Пользователь №: 19 294
|
Цитата(hd44780 @ Nov 24 2012, 18:24) Хотя и как было не работало. Это БПФ рабочее, даже шедевральное в своем роде . Тест я зачем привел? Но такую специфическую форму я написал просто что-бы выстебнуться. Степени 2 с битреверсом в интернете полно БПФ. Хотя целочисленное я нашел только одно. Но вот во всех учебниках пишут - а возможна форма БПФ и без битреверса. Ну я и попытался такую написать. Получилось. Без буферных массивов тут не получится. Их можно было бы внутри объявить, но тут прикол такой - в зависимости от количества точек конечный результат будет то во входных массивах то в буферных. Потому указатели нужны что-бы знать знать откуда читать а куда писать. Внутренние финты с указателями просто что-бы отказаться от if...else тормозящей конвеер и/или сбрасывающей кэш. У Вас же ARM7? Но Вы забейте на все это если лень разбираЦЦо. Вставьте туда (в 2RealFFT) обычное БПФ с битреверсом(оно там уже вставлено). Просто поперезначайте все float на int и битовые сдвиги не забудьте.
--------------------
Нас помнят пока мы мешаем другим...
//--------------------------------------------------------
Хороший блатной - мертвый...
//--------------------------------------------------------
Нет старик, это те дроиды которых я ищу...
|
|
|
|
|
|
|
|
Nov 25 2012, 11:35
|
Профессионал
Группа: Свой
Сообщений: 1 202
Регистрация: 26-08-05
Из: Донецк, ДНР
Пользователь №: 7 980
|
Цитата(DRUID3 @ Nov 24 2012, 19:24) У Вас же ARM7? Да. ARM7TDMI. Цитата(DRUID3 @ Nov 24 2012, 19:24) Вставьте туда (в 2RealFFT) обычное БПФ с битреверсом(оно там уже вставлено). Просто поперезначайте все float на int и битовые сдвиги не забудьте. Я сейчас сделал так (по аналогии с плавающей точкой): Код int a_inReal1[FFT_N/2], a_inReal2[FFT_N/2];
int a_outRe1[FFT_N/2], a_outIm1[FFT_N/2];
int a_outRe2[FFT_N/2], a_outIm2[FFT_N/2];
.....................
// Выполнение FFT
// pInData - входные данные (выборка FFT) - FFT_N
// оконная функция: 0 - нет; 1-Hanning; 2 - Hamming
void DoFFT (byte *pInData, byte iWindow)
{
int i;
// 1. Заполнение входных массивов
for ( i = 0; i < FFT_N/2; i ++ )
{
a_inReal1[i] = pInData [ i ];
a_inReal2[i] = pInData [ i + FFT_N/2 ];
} // for
// 2. FFT
fn_a_2RealFFT_int ( FFT_N/2, a_inReal1, a_inReal2,
a_outRe1, a_outIm1, a_outRe2, a_outIm2,
cai_TTransf, STEP ( FFT_N/2, dfi_3N4 ) );
// memset ( Spectrum, 0, FFT_N);
// 3. Расчёт спектра
for ( i = 0; i < FFT_N/2; i ++ )
{
long Re2, Im2;
long long s;
Re2 = a_outRe1 [ i ];
Im2 = a_outIm1 [ i ];
s = Re2*Re2 + Im2*Im2;
Spectrum [i] = sqrt_int ( s );
Re2 = a_outRe2 [ i ];
Im2 = a_outIm2 [ i ];
s = Re2*Re2 + Im2*Im2;
Spectrum [ FFT_N/2 + i ] = sqrt_int ( s );
} // for
} // DoFFT
// Целочисленный sqrt
unsigned long sqrt_int ( long long l )
{
unsigned long div = 1, rslt = 0;
while (l > 0)
{
l -= div;
div += 2;
rslt += (l < 0) ? 0 : 1;
} // while
return rslt;
} // sqrt_int Тригонометрию подключил из IntTrigArray.h. sqrt_int - нашёл на просторах ... fn_a_2RealFFT_int - переделанный из плавающей точки. Переделки только здесь: Код /*----------------------------------------------------------------------------*/
/* умножение 2-х комплексных чисел - бабочка -4-*/
/*----------------------------------------------------------------------------*/
temp_Re = (((acc_trigT[i_TargC]*a_Re[i_top])>>14) + ((acc_trigT[i_TargS]*a_Im[i_top])>>14))>>1;
temp_Im = (((acc_trigT[i_TargC]*a_Im[i_top])>>14) - ((acc_trigT[i_TargS]*a_Re[i_top])>>14))>>1;
// temp_Re = acc_trigT[i_TargC]*a_Re[i_top] - acc_trigT[i_TargS]*a_Im[i_top]*(cc_dir);
// temp_Im = acc_trigT[i_TargC]*a_Im[i_top] + acc_trigT[i_TargS]*a_Re[i_top]*(cc_dir); Закомментарено - оригинал из плавающей точки. Заменил на аналогичный кусок из fn_aT_ditNbrRadix2FFT_int ... Результат - тормоза при расчёта спектра. Даже нельзя понять, правильн ли рисует Замучился уже
--------------------
Чтобы возить такого пассажира, необходим лимузин другого класса.
(с) Мария Эдуарда
|
|
|
|
|
|
|
|
Nov 25 2012, 11:59
|
山伏
Группа: Свой
Сообщений: 1 827
Регистрация: 3-08-06
Из: Kyyiv
Пользователь №: 19 294
|
Цитата(hd44780 @ Nov 25 2012, 13:35) Замучился уже ...но с другой стороны, а что из того что само идет к нам в руки мы ценим? Цитата(hd44780 @ Nov 25 2012, 13:35) Результат - тормоза при расчёта спектра. Даже нельзя понять, правильн ли рисует А ударить дельтой таки не судьба? В идеале оно выдаст все спектральные бины равны "1"(32768 для int). Код (l < 0) ? 0 : 1 Тоже кэш будет сбрасывать и конвеер тормозить(в 2-3 раза). Лучше ln(sqrt(x))/ln(10) в виде полиномов получить и взять 2-3 первых члена, причем множить "по Горнеру". Но Вы не спешите с этим раз точно не знаете корректно ли заработало FFT.
--------------------
Нас помнят пока мы мешаем другим...
//--------------------------------------------------------
Хороший блатной - мертвый...
//--------------------------------------------------------
Нет старик, это те дроиды которых я ищу...
|
|
|
|
|
|
|
|
Nov 30 2012, 10:56
|
Профессионал
Группа: Свой
Сообщений: 1 202
Регистрация: 26-08-05
Из: Донецк, ДНР
Пользователь №: 7 980
|
Это снова я . Окончательно переехал на целочисленный вариант, как писал выше. Также гонял, проверял, смотрел цифры через DBGU - COM-порт - терминал. В итоге спектр теперь по размеру 256 байт - ровно половина от исходной выборки (раньше там было до хрена нулей - следы от неиспользуемых половинок выходных буферов). Также сделал целочисленный квадратный корень. В итоге длительность работы - FFT+спектр - 18-19 мс на 512 выборке. Хорошо это или плохо - не знаю. Оконной функции пока нет. С отрисовкой проблемы. Точнее не с отрисовкой, а с расчётом высоты полосок. Вот лог отработки типичного блока (98 полосок): Код val=6179 --> Line: (0,48)-(0,0);
val=296 --> Line: (1,48)-(1,0);
val=426 --> Line: (2,48)-(2,0);
val=186 --> Line: (3,48)-(3,0);
val=101 --> Line: (4,48)-(4,0);
val=382 --> Line: (5,48)-(5,0);
val=86 --> Line: (6,48)-(6,0);
val=73 --> Line: (7,48)-(7,0);
val=45 --> Line: (8,48)-(8,5);
val=20 --> Line: (9,48)-(9,30);
val=9 --> Line: (10,48)-(10,41);
val=7 --> Line: (11,48)-(11,43);
val=29 --> Line: (12,48)-(12,21);
val=18 --> Line: (13,48)-(13,32);
val=18 --> Line: (14,48)-(14,32);
val=12 --> Line: (15,48)-(15,38);
val=18 --> Line: (16,48)-(16,32);
val=23 --> Line: (17,48)-(17,27);
val=8 --> Line: (18,48)-(18,42);
val=13 --> Line: (19,48)-(19,37);
val=19 --> Line: (20,48)-(20,31);
val=13 --> Line: (21,48)-(21,37);
val=13 --> Line: (22,48)-(22,37);
val=10 --> Line: (23,48)-(23,40);
val=11 --> Line: (24,48)-(24,39);
val=22 --> Line: (25,48)-(25,28);
val=11 --> Line: (26,48)-(26,39);
val=29 --> Line: (27,48)-(27,21);
val=2 --> Line: (28,48)-(28,48);
val=7 --> Line: (29,48)-(29,43);
val=9 --> Line: (30,48)-(30,41);
val=5 --> Line: (31,48)-(31,45);
val=2 --> Line: (32,48)-(32,48);
val=16 --> Line: (33,48)-(33,34);
val=8 --> Line: (34,48)-(34,42);
val=17 --> Line: (35,48)-(35,33);
val=15 --> Line: (36,48)-(36,35);
val=7 --> Line: (37,48)-(37,43);
val=7 --> Line: (38,48)-(38,43);
val=10 --> Line: (39,48)-(39,40);
val=9 --> Line: (40,48)-(40,41);
val=12 --> Line: (41,48)-(41,38);
val=4 --> Line: (42,48)-(42,46);
val=13 --> Line: (43,48)-(43,37);
val=8 --> Line: (44,48)-(44,42);
val=15 --> Line: (45,48)-(45,35);
val=16 --> Line: (46,48)-(46,34);
val=11 --> Line: (47,48)-(47,39);
val=13 --> Line: (48,48)-(48,37);
val=6141 --> Line: (49,48)-(49,0);
val=498 --> Line: (50,48)-(50,0);
val=528 --> Line: (51,48)-(51,0);
val=165 --> Line: (52,48)-(52,0);
val=192 --> Line: (53,48)-(53,0);
val=152 --> Line: (54,48)-(54,0);
val=60 --> Line: (55,48)-(55,0);
val=15 --> Line: (56,48)-(56,35);
val=32 --> Line: (57,48)-(57,18);
val=53 --> Line: (58,48)-(58,0);
val=10 --> Line: (59,48)-(59,40);
val=25 --> Line: (60,48)-(60,25);
val=17 --> Line: (61,48)-(61,33);
val=12 --> Line: (62,48)-(62,38);
val=20 --> Line: (63,48)-(63,30);
val=21 --> Line: (64,48)-(64,29);
val=17 --> Line: (65,48)-(65,33);
val=26 --> Line: (66,48)-(66,24);
val=13 --> Line: (67,48)-(67,37);
val=14 --> Line: (68,48)-(68,36);
val=10 --> Line: (69,48)-(69,40);
val=8 --> Line: (70,48)-(70,42);
val=7 --> Line: (71,48)-(71,43);
val=25 --> Line: (72,48)-(72,25);
val=14 --> Line: (73,48)-(73,36);
val=4 --> Line: (74,48)-(74,46);
val=20 --> Line: (75,48)-(75,30);
val=3 --> Line: (76,48)-(76,47);
val=6 --> Line: (77,48)-(77,44);
val=2 --> Line: (78,48)-(78,48);
val=9 --> Line: (79,48)-(79,41);
val=15 --> Line: (80,48)-(80,35);
val=8 --> Line: (81,48)-(81,42);
val=6 --> Line: (82,48)-(82,44);
val=8 --> Line: (83,48)-(83,42);
val=10 --> Line: (84,48)-(84,40);
val=15 --> Line: (85,48)-(85,35);
val=1 --> Line: (86,48)-(86,49);
val=12 --> Line: (87,48)-(87,38);
val=8 --> Line: (88,48)-(88,42);
val=5 --> Line: (89,48)-(89,45);
val=29 --> Line: (90,48)-(90,21);
val=18 --> Line: (91,48)-(91,32);
val=6 --> Line: (92,48)-(92,44);
val=5 --> Line: (93,48)-(93,45);
val=14 --> Line: (94,48)-(94,36);
val=9 --> Line: (95,48)-(95,41);
val=0 --> Line: (96,48)-(96,50);
val=9 --> Line: (97,48)-(97,41);
Draw in FFT buffer: duration 367 ms val - значение из спектра. Я рисую добро в окне 100x50, значит мне надо смасштабировать/отсечение всё под высоту окна 50 точек. Но лезут эти дурацкие "иголки" типа val=6179, которые портят всё. Я понимаю, что м.б. дело в отсутствии оконной функции, но всё равно как-то странно. В чём м.б. дело? PS. До дельта-импульсов пока не добрался. Спасибо.
--------------------
Чтобы возить такого пассажира, необходим лимузин другого класса.
(с) Мария Эдуарда
|
|
|
|
|
|
|
|
Nov 30 2012, 17:07
|
Профессионал
Группа: Свой
Сообщений: 1 202
Регистрация: 26-08-05
Из: Донецк, ДНР
Пользователь №: 7 980
|
Посмотрел на своё изделие под музыку. Реакция хорошая и очень адекватная  . Осталось только с этими всплесками разобраться. Кстати заметил, что они идут только начале каждом половинки спектра. Сейчас ещё оконную фунцию сделаю, гляну. Сдаётся мне, что её надо раздельно накладывать на каждую половинку исходного спектра.
--------------------
Чтобы возить такого пассажира, необходим лимузин другого класса.
(с) Мария Эдуарда
|
|
|
|
|
|
|
|
Nov 30 2012, 21:07
|
山伏
Группа: Свой
Сообщений: 1 827
Регистрация: 3-08-06
Из: Kyyiv
Пользователь №: 19 294
|
Цитата(hd44780 @ Nov 30 2012, 19:07) Посмотрел на своё изделие под музыку. Реакция хорошая и очень адекватная Ну пока это "показометр"... Цитата(hd44780 @ Nov 30 2012, 19:07) Осталось только с этими всплесками разобраться. Это не всплески, это спектр такой. Если у Вас АЦП 10 бит то будут палки ~1000 + где-то неверное масштабирование. Что-бы не было длинных палок его нужно еще и логарифмировать. Цитата(hd44780 @ Nov 30 2012, 19:07) Сдаётся мне, что её надо раздельно накладывать на каждую половинку исходного спектра. Оконная функция "развязывает" частотные бины между собой. Вам она не нужна! При такой длине выборки и с int в качестве данных + логарифмирование - нет смысла. P.S.: 18-19 мс на 512 выборке - пока это очень много что-то.
--------------------
Нас помнят пока мы мешаем другим...
//--------------------------------------------------------
Хороший блатной - мертвый...
//--------------------------------------------------------
Нет старик, это те дроиды которых я ищу...
|
|
|
|
|
|
|
|
Dec 2 2012, 18:42
|
Профессионал
Группа: Свой
Сообщений: 1 202
Регистрация: 26-08-05
Из: Донецк, ДНР
Пользователь №: 7 980
|
Цитата(DRUID3 @ Nov 30 2012, 23:07) Ну пока это "показометр"... Я его и делаю . На большее я пока не претендую .... Слишком мало знаю. Но стараюсь, чтобы хоть правдоподобно это выглядело... Сейчас меня только эти 2 палки угнетают - одна вначале (0-й элемент), 2-я ровно посередине - в логе видно. И их длина практически не меняется. Даже если на входе ничего нет .. Пропускать их при рендеринге что ли  ... Да и сам по себе FFT, на мой взгляд, некий показометр . В том смысле, что какие частоты он анализирует - одному Богу известно, к тому же это прямо зависит от размера выборки. Цитата(DRUID3 @ Nov 30 2012, 23:07) Это не всплески, это спектр такой. Если у Вас АЦП 10 бит то будут палки ~1000 + где-то неверное масштабирование. Что-бы не было длинных палок его нужно еще и логарифмировать. Да, АЦП 10 бит. Но в буфер помещаются 8-битные (2 младших бита я отбрасываю сдвигом вправо). Может в этом дело? Цитата(DRUID3 @ Nov 30 2012, 23:07) Оконная функция "развязывает" частотные бины между собой. Вам она не нужна! При такой длине выборки и с int в качестве данных + логарифмирование - нет смысла. Спасибо. Тогда и делать её не буду . Про логарифмирование вопрос - как вообще его делать? Я попытался при расчёте спектра писать Spectrum [i] = 20*(int)log10 ( (float) sqrt_int ( s ) ); // типа dB вместо Spectrum [i] = sqrt_int ( s ); Но это явно фигня - что-то еле заметное внизу болтается . Если 20 оттуда убрать, там вообще ничего не видно. Почему - понятно. Цитата(DRUID3 @ Nov 30 2012, 23:07) P.S.: 18-19 мс на 512 выборке - пока это очень много что-то. Посмотрю ещё ... Может там вывод в UART подмешался. Там же буфера аппаратного нету, всё в лоб .... Один раз я такую "гаву" словил уже .
--------------------
Чтобы возить такого пассажира, необходим лимузин другого класса.
(с) Мария Эдуарда
|
|
|
|
|
|
|
|
1 чел. читают эту тему (гостей: 1, скрытых пользователей: 0)
Пользователей: 0
|
|
|
