FFT на асм для ARM7TDMI (AT91SAM7xx) Опции

[hd44780]

На Си нашёл, вставил в программу - тормозит оно. Конечно, гораздо лучше, чем на AVR, но всё равно не айс.
Надо сделать аудио-анализатор. Сделать-сделал, осциллограммы рисует великолепно, рендер быстрый для дисплея написал, а с FFT проблемы.
А если ещё и стерео запустить - вообще ступор почти будет...

Нет ли у кого реализации на асме, типа как Чен для AVR-ов?
Сам я ейный ассемблер практически не знаю .
Спасибо.

[_Pasha]

На Си нашёл, вставил в программу - тормозит оно.
Спасибо.

Небось, в плавучке?
Тем более, что смотря какой анализатор - набор согласованных фильтров по фиксированной сетке частот может оказаться лучше.

Сообщение отредактировал _Pasha - Nov 15 2012, 08:22

[hd44780]

Нет, не на плавучке
Выдрал у Геннадия - http://www.cqham.ru/forum/showthread.php?t=9688 100% целочисленная арифметика.
Даже слегка пооптимизировал его .
Там есть макросы
#define FFT_POWER 8
#define FFT_N (1 << FFT_POWER) // 256

У него величина FFT_POWER в доп. цикле определялась. Зачем - непонятно, ведь она уже есть ...

А где почитать про "набор согласованных фильтров по фиксированной сетке частот"? Или что-нибудь "для подражания"
Я в этих делах - почти ноль Это типа набора фильтров на определённые частоты?

Вот, нашёл сие - http://www.renan.org/ARM/doc/Apps16pdf.pdf
но, похоже, только общие слова

[Xenia]

Выдрал у Геннадия - http://www.cqham.ru/forum/showthread.php?t=9688 100% целочисленная арифметика.

Уточните, откуда вы его там выдрали? А то я в упор не вижу в той ссылке чего-либо про FFT.

[hd44780]

В конце есть архив исходниками разными - TC1-Oct-30-2012.zip‎ Там есть файлы fft.c и прочее.
Окна Хэмминга и Ханнинга в комплекте .
Добро обнаружил, в общем-то случайно, когда с дисплеем на ILI9320 разбирался.

Я вложил сюда, чтобы там не шарить.

Даже для какого-то ARM9 нашёл на асме - http://www.elsevierdirect.com/v2/companion...N=9781558608740

Вот бабочка для ARM7TDMI - http://www.platan.ru/shem/pdf/bpf.pdf
Остальное обещает на http://www.platan.ru/shem/ , но там Error 404

Прикрепленные файлы

 FFT.ZIP ( 6.31 килобайт ) Кол-во скачиваний: 26

 

[_Pasha]

Нет, не на плавучке

А где почитать про "набор согласованных фильтров по фиксированной сетке частот"? Или что-нибудь "для подражания"
Я в этих делах - почти ноль Это типа набора фильтров на определённые частоты?

Вот, с кортехом-м4 уже не шутошное действо, эта плавучка.
Книжка:
затертых годов... чётта с разбегу не нахожу. Страница 67-74
Да, и поделитесь размером окна FFT, чтоли ...

Сообщение отредактировал _Pasha - Nov 15 2012, 12:10

[hd44780]

Спасибо, почитаю.
Размер окна сейчас - 256 байт. Хоть осциллограмма красивее на 8кБ

Как я выше приводил:
#define FFT_POWER 8
#define FFT_N (1 << FFT_POWER) // 256

[Genadi Zawidowsk...]

Я затащил к себе в проект (на будущее использование) это из конкурсной работы AR1731 одного программиста на приз журнала Circuit Cellar.
Оригинальные файлы здесь. Там не в реальном времени, просто демонстрация принципа работы.

Сообщение отредактировал Genadi Zawidowski - Nov 15 2012, 12:28

[_Pasha]

А заинлайнить BitShift() и Butterfly() слабо?

[hd44780]

Genadi Zawidowski, я так и подумал, т.к. оно у Вас фактически нигде не используются. Даже пытался искать по AR1731 - не нашёл. Спасибо.
_Pasha, попробую вечером дома. потом отпишусь.
Там ещё входной разъём почти отвалился, запаять надо.

[hd44780]

Разъём припаял, inline вставил.
Закрутился чуть быстрее, но всё равно тормоза.
Но даже не это суть. Там полосок 10 вначале прыгают, дальше тишина .... Остальные - редко-редко ... На входе обычная песня (попса ) с компа, частот там валом всяких.
Наверное там сам алгоритм "неправильный" .

[ViKo]

A. N. Sloss ... ARM System Developer's Guide - лежит в закромах.
Начиная со стр 303.

[_Pasha]

Наверное там сам алгоритм "неправильный" .

Дык можно ж свипом 0-20кГц пощщупать..

[hd44780]

Viko, нашёл, скачал...

_Pasha, видимо придётся... А не знаете какую-нибудь прогу для генерации частот? Где-то натыкался, не помню ...

[Genadi Zawidowsk...]

прогу для генерации частот? Где-то натыкался, не помню ...

В том же моём проекте встроена генерация многочастотного сигнала - я с её помощью и проверял алгоритм (то как я его переносил к себе).

Получение очередного семпла с АЦП - GetADCVoltage() - самое начало файла hardware.c - там я даже шум пытался замешать.

[ReAl]

А заинлайнить BitShift() и Butterfly() слабо?

Когда там запилят уже fixed point в gcc? В С11 уже ж вроде попало. Или в 4.8 уже есть?
Сюда

Код

typedef int _Complex cmplx_t;

//cmplx_t a = 2 + 3I;
//cmplx_t b = -1 + 2I;

void foo(cmplx_t *px, cmplx_t *py, cmplx_t w, int len)
{
    while(len--) {
        cmplx_t temp = *px + *py * w;
        *py = *px - *py * w;
        *px = temp;
        ++px;
        ++py;
    }
}

arm-none-eabi-gcc -Os -S -mcpu=cortex-m3 -mthumb cmplx-int.c

CODE

foo:
push {r4, r5, r6, r7, r8, lr}
ldr r4, [sp, #24]
b .L2
.L3:
ldmdb r1, {r5, r8}
mul r7, r2, r5
mls r7, r3, r8, r7
mul r8, r2, r8
mla r5, r3, r5, r8
ldr ip, [r0, #-8]
ldr r6, [r0, #-4]
rsb r8, r7, ip
str r8, [r1, #-8]
add r7, ip, r7
rsb r8, r5, r6
adds r5, r6, r5
str r8, [r1, #-4]
str r7, [r0, #-8]
str r5, [r0, #-4]
.L2:
subs r4, r4, #1
adds r0, r0, #8
adds r1, r1, #8
cmp r4, #-1
bne .L3
pop {r4, r5, r6, r7, r8, pc}

ещё fixed point встроенный, так и вообще хорошо.
(правда это Cortex-M3, тут ещё mla/mls хорошо пошли)



Закат солнца вручную тоже нормально выглядит:

Код

#define FIXED_FACTOR (1<<12)

typedef int _Complex cmplx_t;

void foo(cmplx_t *px, cmplx_t *py, cmplx_t w, int len)
{
    while(len--) {
        cmplx_t temp = *py * w / FIXED_FACTOR;
        *px = *px + temp;
        *py = *px - temp;
        ++px;
        ++py;
    }
}

CODE

foo:
push {r4, r5, r6, r7, r8, lr}
ldr r4, [sp, #24]
b .L2
.L5:
ldmdb r1, {r5, r8}
mul r6, r2, r5
mls r6, r3, r8, r6
mul r8, r2, r8
mla r5, r3, r5, r8
cmp r6, #0
itt lt
addlt r6, r6, #4064
addlt r6, r6, #31
cmp r5, #0
ldr ip, [r0, #-8]
ldr r7, [r0, #-4]
itt lt
addlt r5, r5, #4064
addlt r5, r5, #31
add r6, ip, r6, asr #12
add r5, r7, r5, asr #12
str r6, [r0, #-8]
str r5, [r0, #-4]
str ip, [r1, #-8]
str r7, [r1, #-4]
.L2:
subs r4, r4, #1
adds r0, r0, #8
adds r1, r1, #8
cmp r4, #-1
bne .L5
pop {r4, r5, r6, r7, r8, pc}

(ой, ну по нормальному там не ++px; ++py; а раскрыв крыльев бабочки передавать надо)

[_Pasha]

Вроде ж давно есть fixed...

Да, пока не забыл - бит реверс по приснопамятному красивому трюку
Здесь
Т.е. вместо тяжеленького BitShift пишем хард-кодед вариант

Код

b = ((b * 0x0802LU & 0x22110LU) | (b * 0x8020LU & 0x88440LU)) * 0x10101LU >> 16;

Где там ошибка - сходу вроде нету. Мож, в представлении Q11 косяк? Не вникал, короче.

[hd44780]

По-моему, BitShift это мелочи , я сейчас пытаюсь основной цикл перекрутить ...

Код

  // Main FFT loop
  for ( stc = 0; stc < NumOfStages; stc ++ )
  {
    NumOfButterflys = (1<<stc);
    NumOfBlocks = FFT_N>>(stc+1);

    for ( blc=0; blc < NumOfBlocks; blc ++ )
    {
     base = (1<<(stc+1))*blc;
     for ( bfc=0; bfc < NumOfButterflys; bfc ++ )
     {
      Butterfly ( X+base+bfc, X+base+bfc+NumOfButterflys, NumOfBlocks*bfc );
     } // for
    } // for
  } // for

Butterfly - самая тяжёлая вещь ...

http://www.embeddedsignals.com/ARM.htm - есть всё, но Cortex-M3

[_Pasha]

Butterfly - самая тяжёлая вещь ...

Глянул - точно. Надо избавиться от complex.c - сразу инлайнами сделать операции над комплексными числами. Потом - сделать unroll для хотя бы внутреннего цикла. Оптимизация -О2 или -О3

[hd44780]

inline сделал, слегка видоизменил структуру вызовов кода - пошустрей пошёл, но недостаточно - сильно заметен рендеринг, хотя тот же рендер на осциллограмме летает.

_Pasha, с unroll непонятно - как я могу его сделать, ведь число оборотов циклов заранее неизвестно?
Возможно я не понял, что Вы имели в виду . Поясните пожалуйста.

[_Pasha]

_Pasha, с unroll непонятно - как я могу его сделать, ведь число оборотов циклов заранее неизвестно?
Возможно я не понял, что Вы имели в виду . Поясните пожалуйста.

Вы про unroll loops? Это - можно -O2 сделать - оно автоматом включится, там где выгоднее.
Про перестановку бит - я тут более артист разговорного жанра... надо еще посчитать, что лучше. Но может лучше - табличкой, для 256 отсчетов-то?

[hd44780]

Вы про unroll loops? Это - можно -O2 сделать - оно автоматом включится, там где выгоднее.

Спасибо, попробую.

может лучше - табличкой, для 256 отсчетов-то?

с этим вообще пока не разбирался почти. Гляну ещё.

Вот, сделал замеры времени одного прохода FFT через DBGU и COM-порт:

---> DoFFT start
Create complex data: 1 ms
Window function 0: 0 ms
FFT: 9 ms
Determine spektrum: 0 ms
--> DoFFT end

Window function 0 - прямоугольное окно, т.е. вообще там ничего нет.
FFT - скачет 8-9 ms.
Determine spektrum - всегда 0 - значит, меньше 1 ms.

Эти цифры вообще нормальные? Может там собака в другом месте зарыта ..
У меня сейчас отображает всё вместе - осциллограмму, FFT, ещё уровень рисует - просто максимум по каждому блоку данных АЦП. Еще и часиси рисует .
Щас выкину всё, оставлю только FFT.

[hd44780]

_Pasha, спасибо.

Со скоростью вроде управился, работает довольно быстро.
Осталось, в общем-то главное - работает он как-то ненормально, рисует что-то странное..
Видимо, придётся перебирать разные реализации FFT

[Genadi Zawidowsk...]

Там очень проблемно с переполнением - я у себя переделал оригинальный проект (с потерей скорости - на 64 бит арифметику) для простоты.

[hd44780]

Спасибо.
Вечерком проверю вариант с заменой int на long.
Геннадий, а вы не копались, в каком месте алгоритма эти переполнения лезут? Главный цикл, последуэщий расчёт спектра, ещё что?

[hd44780]

Проверил. Особых тормозов от long int не заметил.
Полосок заметно больше стало, но всё равно WinAmp круче.

Сейчас другие реализации алгоритма буду пробовать ...

[esaulenka]

А какой должна быть разница от замены int на long в случае ARM-а ? :-)

[hd44780]

Медленнее. То же самое как на 8-битном проце складывать-вычитать 16-битные числа.
Ну и т.д. т.п.

[adnega]

Полосок заметно больше стало, но всё равно WinAmp круче.

Может, он в децибелах кажет?!

[DRUID3]

На Си нашёл, вставил в программу - тормозит оно. Конечно, гораздо лучше, чем на AVR, но всё равно не айс.
Надо сделать аудио-анализатор. Сделать-сделал, осциллограммы рисует великолепно, рендер быстрый для дисплея написал, а с FFT проблемы.

Вам нужен RealFFT и аппроксимация формулы для нахождения магнитуды (Р.Лайонс страница 475). По-поводу asm я где-то статью когда-то чью-то бросал сюда на форум. По поводу мегаоптимизаций - откажитесь от битреверса, примените RADIX4 и на асме вручную разверните внутренние циклы (не каскады бабочек). Вряд ли есть куда двигаться дальше в этом вопросе...

кстати раз уж пошли аппроксимации то набор IIR + блочное суммирование под частоту обновления экрана вполне претендует на спасение отца русской демократии. Ну а если уж совсем упростить - полосовые фильтры можно строить на основе фильтров скользящего среднего (можно ведь и в дифференцирующем включении их реализовать). Но результат будет... ммм... ну главное что он таки будет.

[hd44780]

Может, он в децибелах кажет?!

хрен его знает... Вообще вряд ли.
Я где-то год назад родственный анализатор на C#, .NET 2.0 (VS.NET 2005) написал.
Алгоритм откуда-то выдрал готовый. Как вы понимаете, всё в лоб, плавающая точка в любых количествах и всё прочее . Логарифмов вроде нету. Результат визуально близок к винампу.
Могу выложить, если кого интересует.

Кстати, как децибелы по таблице просчитать, я уже придумал, осталось реализовать и проверить .

DRUID3, я RealFFT для 80x86 только нашёл - неактуально .

Асма я этого не знаю. Учить потихоньку начал, но когда это будет ....
Нашёл реализацию FFT на нём - вложение. Для какого проца - не знаю, это с сидюка к какой-то буржуйской книжке. К какой - не знаю. На форумах англоязычных нашёл.
Сам определить я пока не способен. Если кто ориентируется - гляньте плиз.
Меня смущают имена функций типа fft_16_arm9m ....

Сейчас курю это - http://code.google.com/p/falab/
Ещё лежит kissFFt, руки пока не дошли ....

Спасибо.

Прикрепленные файлы

 ch08_dsp.zip ( 58.64 килобайт ) Кол-во скачиваний: 12

 

[esaulenka]

Медленнее. То же самое как на 8-битном проце складывать-вычитать 16-битные числа.
Ну и т.д. т.п.

Задам вопрос по-другому.
Какой разрядности у Вас int и long ?

И компилятор какой?

[hd44780]

Компилятор - IAR
int - 32 (вроде бы). Для 64 написал long long.

Вообще у IAR help какой-то тупой, даже размер типов данных я там не нашёл .
Хотя в том же CvAVR и прочих это с пол-пинка находится .....

[_Артём_]

Вообще у IAR help какой-то тупой, даже размер типов данных я там не нашёл .
Хотя в том же CvAVR и прочих это с пол-пинка находится .....

Смотрите EWARM_DevelopmentGuide.ENU.pdf
Part 2. Reference information
Data representation
Basic data types

[hd44780]

Смотрите EWARM_DevelopmentGuide.ENU.pdf

спасибо. Так и есть,long long - 64 бита.
Вечером ещё разок проверю на железе ..
Рассобачился я на этих __int32 и __int64 в VS.NET . Зато как удобно ...

[esaulenka]

Рассобачился я на этих __int32 и __int64 в VS.NET . Зато как удобно ...

Что, и stdint.h нет?


Хех, а я уж решил, что свершилось чудо, и поведение int и long (ещё вчера упоминался long int, не long long) на 32-битниках как-то различается :-)

[hd44780]

long long поставил.
Объём кода подскочил где-то с 33кил до 38 кил
Полосок больше стало, но все равно редкие какие-то ...

[DRUID3]

DRUID3, я RealFFT для 80x86 только нашёл - неактуально .

Прикрепленные файлы

 aT_RealFFTrFFT_GCC_TEST.zip ( 29.4 килобайт ) Кол-во скачиваний: 18
 aT_ditNbrRadix2FFT_int_c_GCC_TEST.zip ( 25.02 килобайт ) Кол-во скачиваний: 18

 

[hd44780]

DRUID3, спасибо, буду смотреть. Попытаюсь к иару прикрутить.

Кста. Вот тут - http://electronix.ru/forum/index.php?showt...ost&p=57271 вчера поздно вечером что-то нарыл, добавил тупо в иаровский проект - ругается, синтаксис на такой.

[DRUID3]

DRUID3, спасибо, буду смотреть. Попытаюсь к иару прикрутить.

...там это... для "реал" они во флоат, но никто не мешает принять синус равным 32768 на время военного конфликта в том числе и в тригонометрических формулах. Но вообще Вам сподручнее будет не само realFFT брать, а функцию 2-а реал за одно комплексное двойной длинны(есть там и такая). Кстати, а сколько у Вас точек? Для правильного масштабирования бабочек я там еще свое intFFT приложил. fixed 0.15 не применял, и так хорошо. С asm не спешите заморачиваться. Вы как вычисляете логарифм магнитуды(или просто магнитуду)бинов в частотной области? Там эти вычисления отъедят намного больше MIPS чем FFT на C.

[hd44780]

Точек сейчас - 512.
Про логарифм магнетуды, извините не понял. Наверное терию вопроса плохо знаю . Поясните пожалуйста.

Ваш код (целочисленный) сейчас читаю.

Как я понял, у Вас там 2 функции - int fn_aT_ditNbrRadix2FFT_int и fn_aT_ditNbrRadix2ReFFT_int - где какой (имею в виду "за одно комплексное двойной длины")?

Спасибо.

[DRUID3]

Про логарифм магнетуды, извините не понял. Наверное терию вопроса плохо знаю . Поясните пожалуйста.

Ну вот Вы график потом отрисовываете. Это график чего? У Вас после 512 точек действительной последовательности будет 512 точоек комплексного спектра. Но рисуете то Вы физический а не комплексный спектр. Т.е. длину вектора образованного каждым комплексным бином - корень из квадрата реал плюс квадрат имеджинари (sqrt(I*I+Q*Q)). Ну и такой спектр еще очень любят потом логарифмировать для качественной визуальной оценки. Корень квадратный и логарифмирование запросто сожрут все эти попытки оптимизации FFT. Тем более не понятно где Вы их взяли для int. Сами писали?

Точек сейчас - 512.

Ваш код (целочисленный) сейчас читаю.

Как я понял, у Вас там 2 функции - int fn_aT_ditNbrRadix2FFT_int и fn_aT_ditNbrRadix2ReFFT_int - где какой (имею в виду "за одно комплексное двойной длины")?

Спасибо.

целочисленный - это просто 2-а комплексных БПФ. Прямое и обратное. Отличаются знаком комплексной экспоненты и внесением множителя 1/n для обратного БПФ в одиночный битовый сдвиг при каждом проходе каскада.

А нужно еще завернуть такое БПФ для получения БПФ вещественной последовательности. Это другой архив. Там ряд функций. БПФ прямое и обратное для действительной последовательности с отбрасыванием зеркально-сопряженной симметрии и без. Но Вам нужна функция fn_a_2RealFFT() делающая 2 БПФ действительного ряда чисел за раз внутри которой комплексное БПФ такой же длинны. Переделайте эту функцию для int(тривиально) и внутрь вставте либо мою(из того другого архива) целочисленную комплексную БПФ либо чью-то другую.

[hd44780]

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

Понял, о чём Вы .
Спектр не логарифмический, логарифмов вообще нигде нету ... И без корней:

#define FMUL(a,b,q) (((long)(a)*()>>(q))
enum { FFT_Q=11 }; //exponent
#define FFT_POWER 9
#define FFT_N (1 << FFT_POWER) // 512 - размер выборки

...........
for ( i = 0; i < FFT_N; i ++ )
{
// Scaling by FFT_N/2=128
fix re = FFTResult[i].Re >> (FFT_POWER / 2);
fix im = FFTResult[i].Im >> (FFT_POWER / 2);
fix Re2, Im2;

// Squareing
Re2 = FMUL ( re, re, FFT_Q );
Im2 = FMUL ( im, im, FFT_Q );
Spectrum [i] = Re2 + Im2;
} // for

Писал не сам, выдрал, по сути, отсюда - http://electronix.ru/forum/index.php?showt...t&p=1111653 и там выше по теме ..
Насколько это корректно, не знаю . Передирал, как зелёный студент ....

целочисленный - это просто 2-а комплексных БПФ. Прямое и обратное. Отличаются знаком комплексной экспоненты и внесением множителя 1/n для обратного БПФ в одиночный битовый сдвиг при каждом проходе каскада.

А нужно еще завернуть такое БПФ для получения БПФ вещественной последовательности. Это другой архив. Там ряд функций. БПФ прямое и обратное для действительной последовательности с отбрасыванием зеркально-сопряженной симметрии и без. Но Вам нужна функция fn_a_2RealFFT() делающая 2 БПФ действительного ряда чисел за раз внутри которой комплексное БПФ такой же длинны. Переделайте эту функцию для int(тривиально) и внутрь вставте либо мою(из того другого архива) целочисленную комплексную БПФ либо чью-то другую.

Щас буду смотреть ...

[DRUID3]

Понял, о чём Вы .
Спектр не логарифмический, логарифмов вообще нигде нету ... И без корней:

Понял. Все уже упрощено дальше некуда. разве что маску наложить.

---> DoFFT start
Create complex data: 1 ms
Window function 0: 0 ms
FFT: 9 ms
Determine spektrum: 0 ms
--> DoFFT end

Window function 0 - прямоугольное окно, т.е. вообще там ничего нет.
FFT - скачет 8-9 ms.
Determine spektrum - всегда 0 - значит, меньше 1 ms.

Эти цифры вообще нормальные? Может там собака в другом месте зарыта ..

А какая у Вас тактовая? В одном из приведенных Вами же документов - 40 MHz ~ 3ms.

Тут у меня такая мысль. Вот этот тип данных complex. Это же структура с двумя полями, так? Т.е. перед каждым обращением к полю нужно вначале перейти на указатель на структуру. А это MIPSы. Инлайнь там что или нет. Попробуйте с обыкновенными масивами I и Q.

[hd44780]

А какая у Вас тактовая? В одном из приведенных Вами же документов - 40 MHz ~ 3ms.

48 MHz. Для USB. Хоть ещё ни разу его не использовал .

Тут у меня такая мысль. Вот этот тип данных complex. Это же структура с двумя полями, так? Т.е. перед каждым обращением к полю нужно вначале перейти на указатель на структуру. А это MIPSы. Инлайнь там что или нет. Попробуйте с обыкновенными масивами I и Q.

попробую. Спасибо.
Я сейчас Ваш смотрю.
Мой вариант по скорости уже приемлемо крутитися, но результаты какие-то странные показывает - от силы 10 полосок из 128, которые я отрисовываю на дисплее. Могу фотку дисплея кинуть вечером. Может "симптоматика" что-то подскажет.

[hd44780]

DRUID3, а чем отличаются fn_aT_ditNbrRadix2FFT_int и fn_aT_ditNbrRadix2ReFFT_int и какую мне лучше использовать?

[DRUID3]

DRUID3, а чем отличаются fn_aT_ditNbrRadix2FFT_int и fn_aT_ditNbrRadix2ReFFT_int и какую мне лучше использовать?

...ну у Вас же прямое преобразование потому первую. Отличаются они знаком комплексной экспоненты и масштабирующим коэффициентом 1/N. Я его внес в каскады в виде битовых сдвигов прямого преобразования - так спектр в масштабе сигнала получается. Т.е. если синус амплитудой 16384, то и палка от него будет такой же.

Кстати Вы пишете, что картинка Вас не радует на винамп не похожа. Так в винампе она в логарифмическом масштабе. А в прямом там действительно только основные гармоники будут.

P.S.: ...там мое тоже не ахти по оптимальности. Нет битреверса потому использованы 2-а дополнительных массива. Но для того чтобы "сэкономить место" используется факт, что синус это косинус + pi/2, а это дополнительные сложения. Плюс там еще "прыжки по таблице тригонометрии" - это чтобы "маленькие" FFT использовали таблицу для "больших" - и это жрет ресурс. Можно сделать битреверс и "вернуться к классике" - это я просто так экспериментировал. В сети полно FFT степени 2 с битреверсом, а такой формы нет нигде что-то .

[hd44780]

ну у Вас же прямое преобразование потому первую.

Спасибо. Я ещё смотрю. Будут вопросы, спрошу.

Кстати Вы пишете, что картинка Вас не радует на винамп не похожа. Так в винампе она в логарифмическом масштабе. А в прямом там действительно только основные гармоники будут.

Да, хотелось бы поближе к винампу ...
Там что, надо еще по каждой полоске 10log10 посчитать и взять модуль, тк. dB всегда отрицательные? А отношение к чему?

Написал туда сейчас тупо
Spectrum [i] = (int)(10.0 * log10 ( Re2 + Im2 ) );

вообще ничего не рисует

В сети полно FFT степени 2 с битреверсом, а такой формы нет нигде что-то .

Да отож ...

[hd44780]

Вот сейчас посмотрел Чана - http://elm-chan.org/works/akilcd/report_e.html , логарифмов там нету, а картинка с виду нормальная ...
Неужели квадратный корень так влияет ?

Вот проверил, корень добавил -
Spectrum [i] = (int) sqrt ( (float)Re2 + (float)Im2 );

Вообще ничего не рисует ....

Наверное всё-таки мой алгоритм кривой

[DRUID3]

У Чана логарифмов нету. У него чистое FFT. Корень квадратный навряд-ли так влияет. Может что с FFT.

А Вы ударьте по своему FFT дельта-импульсом. Что будет на выходе?

Да отож ...

...мы не идем проторенным путем...

Код

/*=======================================================================*/
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) */

Вот где-то так...

[hd44780]

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?
Как я понял из алгоритма и описания в конце, что, т.к. я возвожу всё добро в квадрат, то всё равно. Я прав?

Спасибо.

[DRUID3]

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-а БПФ на столько же отсчетов одновременно!!!

Я сомневаюсь, правильно ли я написал первый аргумент макроса STEP. У Вас там 8 == размеру массивов, т.е. FFT_N.

...кажется верно...

Теперь мне надо определить спектр.
Но я не понял, использовать пару a_outRe1/a_outIm1 или a_outRe2/a_outIm2?
Как я понял из алгоритма и описания в конце, что, т.к. я возвожу всё добро в квадрат, то всё равно. Я прав?

Спасибо.

Это 2-а абсолютно независимых спектра от первого и от второго буфера соответственно.

[hd44780]

Неверно. 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

Это 2-а абсолютно независимых спектра от первого и от второго буфера соответственно.

В свете первого замечания понятно.

[DRUID3]

Т.е. я могу представить свою 512-байтовую выборку как 2 256-байтовых, так что ли?

Да, причем внутреннее(комплексное) БПФ будет на 256 точек. В этом то и ускорение.

[hd44780]

Понял.

Тут ещё один вопрос....
Для входного окна размером 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?

Или ещё как-то?

[DRUID3]

Код

ui_N = cuic_RN>>1;

...пусть cuic_RN = 512 длинна двух Real подпоследовательностей. 512 - будет и complexFFT. А спектры по 256. Да I и Q спектры Real последовательностей зеркально сопряжены и имеют по 512 точек. Но зачем они нам? Информацию несут только их половины.

[hd44780]

Вроде понятно.
В итоге получил такое:

Код

  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

Дома на железяке проверю.

[hd44780]

Исправил косяк

Код

  // 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

Заработал. Результат воде ничего, но тормозит .. Единственное - рисуется только левая половина. Наверное с этими частями намутил
Завтра целочисленный вариант попробую прикрутить.

[DRUID3]

...Результат воде ничего...

Сомневаюсь я однако. И повторю свой совет - чем тестить какими-то гармониками наугад лучше ударить дельта-импульсом и посмотреть на выход - RealFFT это тоже ведь линейная система.

1) а a_outRe2[], a_outIm2[] куда дели?

2) a_outIm1 [ FFT_N/2 + i ] - это что такое вообще? Зачем, откуда? У Вас уже есть готовые полуспектры, забудьте Вы о математическом формализме и этой сопряженной симметрии. Вы же хотите так как в Winamp?

[hd44780]

лучше ударить дельта-импульсом и посмотреть на выход - RealFFT это тоже ведь линейная система.

Я не знаю, что такое "дельта-импульс" . И где его взять ... В гугле поискал - внятного ничего не нашёл.
Могу придумать только генерить частоты в каком-нибудь CoolEdit-е и ими тестировать.

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

Про "готовые полуспектры" не понял . Вы хотите сказать, что после Вашего алгоритма не нужно эти корни брать?

Вы же хотите так как в Winamp?

Совершенно верно

[DRUID3]

Я не знаю, что такое "дельта-импульс" . И где его взять ...

Вы очень много в жизни потеряли. А он заслуживает траты на него времени. Ну хотя-бы потому, что ударь Вы им по любой линейной системе и Вы сразу все узнаете о ней - получите ее АЧХ/ФЧХ или т.н."передаточную функцию" если говорить о временнОй области.

В гугле поискал - внятного ничего не нашёл.

...0, 0, 0, 0, 0, 1, 0, 0, 0... - это она(функция)... ну или он(импульс)...

Про "готовые полуспектры" не понял . Вы хотите сказать, что после Вашего алгоритма не нужно эти корни брать?

Нет, уже все хорошо...

[hd44780]

Отписаться вечером не смог ..

Вчера проверил - на плавающей точке работает, но тормозит - 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-м входным буферам? Так что-ли?

[DRUID3]

...забейте... посмотрите как в том примере с int сделан сдвиг(в каких местах) и сделайте точно так же в том с плавающей...

[hd44780]

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 закрыл, им там что-то присваивается, далее не используется.
Потом вообще думаю их убрать ...
Хотя и как было не работало.

[DRUID3]

Хотя и как было не работало.

Это БПФ рабочее, даже шедевральное в своем роде . Тест я зачем привел? Но такую специфическую форму я написал просто что-бы выстебнуться. Степени 2 с битреверсом в интернете полно БПФ. Хотя целочисленное я нашел только одно. Но вот во всех учебниках пишут - а возможна форма БПФ и без битреверса. Ну я и попытался такую написать. Получилось.
Без буферных массивов тут не получится. Их можно было бы внутри объявить, но тут прикол такой - в зависимости от количества точек конечный результат будет то во входных массивах то в буферных. Потому указатели нужны что-бы знать знать откуда читать а куда писать. Внутренние финты с указателями просто что-бы отказаться от if...else тормозящей конвеер и/или сбрасывающей кэш. У Вас же ARM7?

Но Вы забейте на все это если лень разбираЦЦо. Вставьте туда (в 2RealFFT) обычное БПФ с битреверсом(оно там уже вставлено). Просто поперезначайте все float на int и битовые сдвиги не забудьте.

[hd44780]

У Вас же ARM7?

Да. ARM7TDMI.

Вставьте туда (в 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 ...

Результат - тормоза при расчёта спектра. Даже нельзя понять, правильн ли рисует
Замучился уже

[DRUID3]

Замучился уже

...но с другой стороны, а что из того что само идет к нам в руки мы ценим?

Результат - тормоза при расчёта спектра. Даже нельзя понять, правильн ли рисует

А ударить дельтой таки не судьба? В идеале оно выдаст все спектральные бины равны "1"(32768 для int).

Код

(l < 0) ? 0 : 1

Тоже кэш будет сбрасывать и конвеер тормозить(в 2-3 раза). Лучше ln(sqrt(x))/ln(10) в виде полиномов получить и взять 2-3 первых члена, причем множить "по Горнеру". Но Вы не спешите с этим раз точно не знаете корректно ли заработало FFT.

[hd44780]

Да рано мне до дельты ...
Нашёл в интернете алгоритмы - http://algolist.manual.ru/maths/count_fast/intsqrt.php
Выдаёт почти всегда нули ... Загнал его сейчас в VS.NET 2005 - то же самое.

[hd44780]

Да рано мне до дельты ...
Нашёл в интернете алгоритмы - http://algolist.manual.ru/maths/count_fast/intsqrt.php
Выдаёт почти всегда нули ... Загнал его сейчас в VS.NET 2005 - то же самое.

[hd44780]

Дубль поста получился ..

[hd44780]

Это снова я .

Окончательно переехал на целочисленный вариант, как писал выше.
Также гонял, проверял, смотрел цифры через 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. До дельта-импульсов пока не добрался.

Спасибо.

[hd44780]

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

[DRUID3]

Посмотрел на своё изделие под музыку. Реакция хорошая и очень адекватная

Ну пока это "показометр"...

Осталось только с этими всплесками разобраться.

Это не всплески, это спектр такой. Если у Вас АЦП 10 бит то будут палки ~1000 + где-то неверное масштабирование. Что-бы не было длинных палок его нужно еще и логарифмировать.

Сдаётся мне, что её надо раздельно накладывать на каждую половинку исходного спектра.

Оконная функция "развязывает" частотные бины между собой. Вам она не нужна! При такой длине выборки и с int в качестве данных + логарифмирование - нет смысла.

P.S.: 18-19 мс на 512 выборке - пока это очень много что-то.

[hd44780]

Ну пока это "показометр"...

Я его и делаю . На большее я пока не претендую .... Слишком мало знаю.
Но стараюсь, чтобы хоть правдоподобно это выглядело... Сейчас меня только эти 2 палки угнетают - одна вначале (0-й элемент), 2-я ровно посередине - в логе видно. И их длина практически не меняется. Даже если на входе ничего нет ..
Пропускать их при рендеринге что ли ...

Да и сам по себе FFT, на мой взгляд, некий показометр .
В том смысле, что какие частоты он анализирует - одному Богу известно, к тому же это прямо зависит от размера выборки.

Это не всплески, это спектр такой. Если у Вас АЦП 10 бит то будут палки ~1000 + где-то неверное масштабирование. Что-бы не было длинных палок его нужно еще и логарифмировать.

Да, АЦП 10 бит. Но в буфер помещаются 8-битные (2 младших бита я отбрасываю сдвигом вправо). Может в этом дело?

Оконная функция "развязывает" частотные бины между собой. Вам она не нужна! При такой длине выборки и с int в качестве данных + логарифмирование - нет смысла.

Спасибо. Тогда и делать её не буду .
Про логарифмирование вопрос - как вообще его делать?

Я попытался при расчёте спектра писать
Spectrum [i] = 20*(int)log10 ( (float) sqrt_int ( s ) ); // типа dB
вместо
Spectrum [i] = sqrt_int ( s );
Но это явно фигня - что-то еле заметное внизу болтается . Если 20 оттуда убрать, там вообще ничего не видно. Почему - понятно.

P.S.: 18-19 мс на 512 выборке - пока это очень много что-то.

Посмотрю ещё ... Может там вывод в UART подмешался. Там же буфера аппаратного нету, всё в лоб ....
Один раз я такую "гаву" словил уже .

[hd44780]

Почистил мусор всякий - расчёт FFT вместе со спектром - 7-8ms по 512-байтовой выборке.
От палок этих вроде избавился, пока на этом остановлюсь.

Много чего ещё надо сделать - оцифровку обоих каналов, логарифмический индикатор уровня хотел добавить...
DRUID3, ещё раз большое Вам спасибо.

[blackfin]

Почистил мусор всякий - расчёт FFT вместе со спектром - 7-8ms по 512-байтовой выборке.

Всё равно медленно. BF527 делает FFT на 512 точках за 11 мкс.

[hd44780]

Всё равно медленно. BF527 делает FFT на 512 точках за 11 мкс.

А на какой частоте он у Вас работает? Глянул на страничке - он до 600 MHz. Детальнее не смотрел. Может ещё и плавающий сопроцессор есть.
Тогда там вообще эта ерундистика с целочисленным табличным синусом/косинусом плюс необходимые коррекции и на фиг не нужна ...
Я тут собираюсь заказать STM32F4DISCOVERY, там частота вроде 192 MHz (сорри, если ошибся) + FPU - вот и проверю.

ЗЫ.
Core2Quad (даже простой 2-х ядерный пень) ещё быстрее считает

[DRUID3]

Почистил мусор всякий - расчёт FFT вместе со спектром - 7-8ms по 512-байтовой выборке.
От палок этих вроде избавился, пока на этом остановлюсь.

...рано. Что-то уж очень большие цифры. Тем более, что мы(Вы) используем RealFFT.

По-поводу логарифмирования и корня:
sqrt() - это же число в степени 0.5 - а она легким движением руки переносится в показатель степени основания логарифма или просто выносится за логарифм. Потому никакой целочисленный корень не нужен.

Кстати, это оффтоп, но у меня есть веские основания полагать почему так развита была теория "кепструмов" для распознавания речи в эпоху малопроизводительных ЭВМ без float - а просто легкое целочисленное вычисление - смотрим дальше по тексту.

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


Это не я такой умный, а Pilot с "винграда". У него на форуме я такое нашел. Правда тогда я и не допер этот способ, а потом допер и вышла лабуда. А почему лабуда? А потому, что мы ищем не логарифм от 1 до 32768, а логарифм 0.001 до 0.999, а уже эту сетку домножаем на как по x 32768 так и по y.

Эскизы прикрепленных изображений