БПФ, может есть у кого исходник? Опции

[Pianist]

Всем добрый вечер. Хочу реализовать простенький анализатор спектра на AVR. Микроконтроллер пока не выбрал. Может есть у кого программа быстрого преобразования фурье на С? Поделитесь пожалуйста...

Сообщение отредактировал Pianist - Mar 14 2009, 20:11

[Alex11]

Посмотрите. Это, конечно, под PC, но там полторы страницы текста, можно оптимизировать и под контроллер.  FFT32k.zip ( 1010 байт ) Кол-во скачиваний: 938

[klen]

тут все уже сделано, используются инструкции арифметики дробная фиксированой запятой

http://elm-chan.org/works/akilcd/report_e.html

[Pianist]

Спасибо. Но все же хочется самому попробовать. Насколько я понял без операции умножения мне далеко не уехать, а из микроконтроллеров AVR аппаратное умножение появилось только у АтМега128?

[Rst7]

умножение появилось только у АтМега128?

Все современные AVR серии Mega имеют аппаратный умножитель. Вы давно не посещали сайт Атмела

[blackfin]

Может есть у кого программа быстрого преобразования фурье на С? Поделитесь пожалуйста...

БПФ

[ARV]

тут все уже сделано, используются инструкции арифметики дробная фиксированой запятой
http://elm-chan.org/works/akilcd/report_e.html

+5
на чистом Си для AVR можно и не начинать - скорости не хватит для реального времени

не могу удержаться... вот что с помощью этой либы сделал

[Rst7]

на чистом Си для AVR можно и не начинать - скорости не хватит для реального времени

Да ладно. Для такого спектроанализатора достаточно было сделать 10 гетеродинных фильтров. Без всякого Фурье. По скорости это, ну, на вскидку, не более 15-20% CPU Load.

[GDI]

Вот, на асме, правда, 64 точки 8 бит DFT. Не моё, сам где то когда то скачал.

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

 DFT864V1.zip ( 3.6 килобайт ) Кол-во скачиваний: 133
 GLOBDEFS.zip ( 6.32 килобайт ) Кол-во скачиваний: 118
 MathCons.zip ( 2.45 килобайт ) Кол-во скачиваний: 121

 

[ARV]

Да ладно. Для такого спектроанализатора достаточно было сделать 10 гетеродинных фильтров. Без всякого Фурье. По скорости это, ну, на вскидку, не более 15-20% CPU Load.

речь идет об анализаторе или о реализации БПФ на Си? мое мнение о нехватке производительности относилось к реализации на Си БПФ для AVR.

на счет вашего варианта не могу сказать ничего, хотя сомнения все-таки есть... особенно если эти фильтры так же на чистом Си делать...

[Rst7]

мое мнение о нехватке производительности относилось к реализации на Си БПФ для AVR.

Ну вот в теме про распознование голоса мы делали DHT, правда, 32хточечный. На чистом Си. Не померли 2000 тактов.

на счет вашего варианта не могу сказать ничего, хотя сомнения все-таки есть... особенно если эти фильтры так же на чистом Си делать...

Да нефиг там делать. Все просто, как стол.

[ARV]

Ну вот в теме про распознование голоса мы делали DHT, правда, 32хточечный. На чистом Си. Не померли 2000 тактов.

гм... нашел ваши исходники... все-таки ваш код можно назвать "чистым Си" с весьма большой натяжкой - там все сделано не только в виде оптимизации алгоритма, но и специально заточено на IAR... в сущности, доказывать, что при глубоком знании особенностей компилятора и умелом их использовании можно получить очень близкий к чисто ассемблерному варианту код - бессмысленно. я же имел ввиду под чистым Си, когда возведение в степень так и делается long result = var * var;
т.е. ваше решение "не в лоб", а я говорил о "лобовой атаке" все расставлено по своим местам.

а вот по поводу полосовых фильтров - увы, в ЦОС не шибко силен... считаете, что 10 цифровых полосовых фильтров будут быстрее БПФ?

[Rst7]

но и специально заточено на IAR...

IAR пока является наиболее оптимальным компилятором для AVR (несмотря на некоторые свои забобоны). Так что имею право. А если бы тут был ARM, то функции, помеченные как hard-coded функции заменяеются на прямой код без извращений и все нормально работает.

Т.е. например

Код

//Hardcoded for AVR, return (var*c)>>8;
#pragma inline=forced
static int FDMUL(int var, UREG c)
{
  int v=__multiply_signed_with_unsigned(var>>8,c);
  unsigned int v2=__multiply_unsigned(var,c);
  v+=v2>>8;
  return v;
}

заменяется на

Код

#pragma inline=forced
static int FDMUL(int var, UREG c) {return (var*c)>>8;}

И всем славно.

считаете, что 10 цифровых полосовых фильтров будут быстрее БПФ?

Конечно. Потому что отсчетов для БПФ надо брать намного больше - необходимые полосы на каждой интересующей частоте разной ширины - это легко решается в случае фильтров.

[ARV]

Потому что отсчетов для БПФ надо брать намного больше - необходимые полосы на каждой интересующей частоте разной ширины - это легко решается в случае фильтров.

больше 32? я обошелся 32 отсчетами... но с полосовыми хотелось бы разобраться... что-нибудь посоветуете? если можно, "для средних умов"

[Rst7]

больше 32?

Ну и как тогда получить с частотами

32, 64, 125, 250, 500 Гц и 1, 2, 4, 8, 15,5

правильные полосы? 32 точки с нижней частотой в 32Гц имеют верхнуюю полосу в 1024Гц. А никак не 15.5.

[ARV]

правильные полосы? 32 точки с нижней частотой в 32Гц имеют верхнуюю полосу в 1024Гц. А никак не 15.5.

частоты правильные. я же не говорил, что я из одного комплекта в 32 отсчета все частоты получаю
а что на счет моей просьбы?

[Rst7]

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

А это что:

я обошелся 32 отсчетами...

Так сколько у Вас там БПФ делается по 32 точки каждое?

но с полосовыми хотелось бы разобраться... что-нибудь посоветуете? если можно, "для средних умов"

А че там разбираться. Инициализируем sum_sin=0 и sum_cos=0. На каждом отсчете выполняем

Код

sum_sin+=signal*sin(phase*2*pi);
sum_cos+=signal*cos(phase*2*pi);
phase+=freq;

После n отсчетов делаем power=sqrt(sum_sin^2+sum_cos^2);

Итого, имеем полосовой фильтр плюс выпрямитель RMS с полосой Fдискретизации/n и средней частотой Fдискретизации*freq (вроде нигде не ошибся).

При реализации нескольких октавных фильтров можно оптимайзить в районе нехранения phase и freq для каждой частоты (для каждого более высокочастотного фильтра phase[n+1]=phase[n]*2).

[ARV]

Так сколько у Вас там БПФ делается по 32 точки каждое?

2

А че там разбираться. Инициализируем sum_sin=0 и sum_cos=0. На каждом отсчете выполняем

Код

sum_sin+=signal*sin(phase*2*pi);
sum_cos+=signal*cos(phase*2*pi);
phase+=freq;

После n отсчетов делаем power=sqrt(sum_sin^2+sum_cos^2);

Итого, имеем полосовой фильтр плюс выпрямитель RMS с полосой Fдискретизации/n и средней частотой Fдискретизации*freq (вроде нигде не ошибся).

При реализации нескольких октавных фильтров можно оптимайзить в районе нехранения phase и freq для каждой частоты (для каждого более высокочастотного фильтра phase[n+1]=phase[n]*2).

спасибо, покумекаю на досуге...

[Rst7]

2

И дальше что? Полосы суммируете? Т.е. из 32х мощностей берете 1,2+3,4+5+6+7,... и т.д.?

[ARV]

И дальше что? Полосы суммируете? Т.е. из 32х мощностей берете 1,2+3,4+5+6+7,... и т.д.?

нет, не суммирую, просто беру нужные гармоники.

[Rst7]

нет, не суммирую, просто беру нужные гармоники.

Ну и что это получается? Что на частоте, скажем 8кГц вы обрабатываете только полосу в 1кГц. Т.е. между полосами у Вас дырки, причем эти дырки шире самих полос. Фигня-с

[GDI]

А че там разбираться. Инициализируем sum_sin=0 и sum_cos=0. На каждом отсчете выполняем

Код

sum_sin+=signal*sin(phase*2*pi);
sum_cos+=signal*cos(phase*2*pi);
phase+=freq;

После n отсчетов делаем power=sqrt(sum_sin^2+sum_cos^2);

Уважаемый rst7 это не классический ли ДПФ получается, только там где то нужно умножить на 1/N. И с аргументом синусов и косинусов что то не то... Надо сделать как то так:

Код

sum_sin=0;
sum_cos=0;
for (i=0; i <N; i++){
  sum_sin+=signal*sin(freq*2*pi/N);
  sum_cos+=signal*cos(freq*2*pi/N);
}
power=(1/N)*sqrt(sum_sin^2+sum_cos^2);

Где N это длина ДПФ

[Rst7]

Уважаемый rst7 это не классический ли ДПФ получается

Конечно. Это одна гармоника ДПФ.

только там где то нужно умножить на 1/N.

Ну это один раз после всех операций. Сложность аж О(1), так что можно не учитывать. Но помнить.

[ARV]

Ну и что это получается? Что на частоте, скажем 8кГц вы обрабатываете только полосу в 1кГц. Т.е. между полосами у Вас дырки, причем эти дырки шире самих полос. Фигня-с

во-первых, для поставленной цели - самое нормальное. во-вторых, ширина индицируемой полосы (т.е. добротность фильтра обычными словами) вполне удовлетворительная, чтобы захватывать не 8К, а диапазон скажем от 7,5К до 8,2К. наконец, целочисленная арифметика дает, мягко говоря, небольшую погрешность, еще более размазывающую полосу... ну и, в-четвертых, наложение оконной функции никто не отменял (хотя мне и без нее хорошо)...
Кстати, как вы себе представляете спектроанализатор без дырок между полосами?! Это уже будет показометр - как в китайских магнитолах, просто фиксированные "картинки"

P.S. выходит, ваш "полосовой фильтр" ничем от БПФ не отличается... а я уж подумал было...

[Rst7]

наконец, целочисленная арифметика дает, мягко говоря, небольшую погрешность, еще более размазывающую полосу...

Для 32х точек эта погрешность исчезающе мала.

А правильный октавный спектроанализатор должен показывать мощность сигнала в каждой полосе от 0.707*F до 1.41*F, где F - центральная частота. Так что "показометр" именно у Вас.

От противного - октавный эквалайзер имеет именно такие полосы регулировки. А не полосу в 1кГц на центральной частоте 8кГц.

выходит, ваш "полосовой фильтр" ничем от БПФ не отличается... а я уж подумал было...

Отличается. Еще и как. Курите глубже

[ARV]

Для 32х точек эта погрешность исчезающе мала.

не знаю, не знаю... она заметна

А правильный октавный спектроанализатор должен показывать мощность сигнала в каждой полосе от 0.707*F до 1.41*F, где F - центральная частота.

нигде не встречал требований к "правльному" спектроанализатору

Так что "показометр" именно у Вас.

вы первый, кому не понравилось

От противного - октавный эквалайзер имеет именно такие полосы регулировки. А не полосу в 1кГц на центральной частоте 8кГц.

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

Отличается. Еще и как. Курите глубже

я подумаю (хотя это сложно).

[Rst7]

не знаю, не знаю... она заметна

Плохо накодили.

нигде не встречал требований к "правльному" спектроанализатору

Ваш имеет дырки между полосами. Т.е. зоны нечувствительности. Точнее, видимо, из-за ошибок в математике, какой-то мусор пролезает, но не более.

В идеальном случае, если Вы подадите на вход синусоидальный (для простоты) сигнал с любой частотой, проссумировав мощности во всех полосах, Вы должны получить число, зависящее только от мощности входного сигнала, но не от частоты (в пределах рабочей полосы, у Вас от 32Гц до 16кГц).

опять же, где эти нормы прописаны? на счет добротности фильтров эквалайзера? и насколько эти нормы таковыми являются?

Нигде. Однако есть здравый смысл, опирающийся на опыт. Вы когда-нибудь слушали сигнал, в котором поднята одна очень узкая полоса? А такое будет при регулировке тембра очень узкими полосами.

И с аргументом синусов и косинусов что то не то... Надо сделать как то так:

Да все то. Просто по привычке допустил пару тривиальных упрощений. Зачем делить на N каждый раз, если можно разделить 1 раз, и freq и phase. Да и цикл не интересует. Поэтому сделано накопление фазы.

Умножение фазы на 2pi подразумевает, что аргумент phase имеет значение 0...1. Отнюдь не реальное умножение. В реальной жизни phase суть байт, 0...255, следовательно функция задается таблицей со значениями table_sin[index]=sin(index*pi/128).

[ARV]

Rst7, не в порядке спора, а чисто для пояснения.
боюсь, анализатор спектра в том виде, как его задумал я, вряд ли станет удовлетворять эстетические потребности, если его сделать по вашему алгоритму, хотя умозрительно он верный... ведь по-вашему, если на входе присутствует сигнал хоть какой частоты, то анализатор должен что-то показать... так вот, будет очень неубедительно, если флейта станет "индицироваться" столбиками от 1 до 8 кимлогерц - хоть какого уровня... пропадет динамика картины, а именно это главное в подобных девайсах не так ли? мы ведь не о каком-то научном спектральном анализе ведем речь...
что касается математических вычислений, то я их не делал сам, сэкономил силы/время, а воспользовался ранее указанной ссылкой, т.е. ЭлмЧеновской библиотекой - всего лишь подкорректировал ее под свои задачи. А его библиотека, хотя и ведет расчеты в 16-битных целых, страдает погрешностями, например, корень из 9 может оказаться равным 4... (я точно не помню. но когда тестировал - сильно удивлялся). на больших числах погрешность меньше, но все равно где-то до 10-15% запросто достигает... если младшие биты отбрасывать, получается незаметно внешне - я и смирился. кстати, быстродействие ченовской библиотеки ничуть не ниже вашего варианта, и, наверное, даже выше - для 128-точечного БПФ я получал без малого 100fps для своего индикатора 8х8 - причем с учетом семплов АЦП (около 32 ksps)... если память не изменяет, весь цикл был то ли 11, то ли 16 (могу ошибаться) миллисекунд... так что в угоду скорости принесена точность...

[Rst7]

так вот, будет очень неубедительно, если флейта станет "индицироваться" столбиками от 1 до 8 кимлогерц - хоть какого уровня... пропадет динамика картины, а именно это главное в подобных девайсах не так ли? мы ведь не о каком-то научном спектральном анализе ведем речь...

Не далее, как несколько постов назад Вы обвиняли меня в том, что, якобы, у меня "показометр" -

Кстати, как вы себе представляете спектроанализатор без дырок между полосами?! Это уже будет показометр - как в китайских магнитолах, просто фиксированные "картинки"

Я требую сатисфакции

Кстати, о "качестве" картинки. Если Вам потреблястов развлекать мигающими светодиодиками - не возражаю. А вот если Вы свой девайс воткнете в какой-нибудь микшерский пульт - потом может быть стыдно за свою поделку. Объясняю почему - допустим, мы увеличили усиление в какой-то полосе (ну так нам захотелось), при этом пытаемся сориентироваться по Вашему "октавному спектроанализатору", не перегрузим ли мы дальнейший тракт. Однако, оказалось, что максимальная спектральная плотность сигнала совсем не на той частоте, которую показывается Ваш "показометр", а всего-лишь рядом. В результате: Ваш "показометр" показывает, что все окей, перегруза нет, на местном контроле в ушах вроде все хорошо (потому что запасы обычно немалые), а дальше вполне хватило для перегрузки, например, АЦП/усилителя в тракте записи.

[rubic]

кто может занимался с этой библиотекой FFT,у меня много непоняток
вопрос :
Функция fft_input(capture, bfly_buff) должна должна сортировать входные данные(capture) для так называемой "бабочки" в (bfly_buff) с нечётными номерами в дну половину буфера с чётными в другую,у меня же полная били...да получается и вообще bfly_buff с адресом выше 64 она заполняет нулями ,может здесь алгоритм какой другой или я что непонимаю

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

 avrfft.zip ( 22.6 килобайт ) Кол-во скачиваний: 69

 

[ARV]

на сколько я помню, перестановки там делаются уже после FFT, в самом конце расчетов. инпут просто буфер заполняет

[rubic]

После преобразвания тоже но это другое зависит от алгоритма,но все исходники которые мне встречались сначала выстраивают матрицу например для FFT_N=16 {a0,a1,a2....a15} в {a0,a8,a4,a12,a2,a10,a6,a14,a1,a9,a5,a13,a3,a11,a7,a15} для удобства вычисления

[ARV]

чен сделал комплексное БПФ, поэтому мнимые части просто обнуляются. по-моему, так. а перестановки у него в библиотеке после всего сделаны были

[rubic]

по-моему

это по твоему,а в реале спектр только 32 значения склько выброк неделай

[ARV]

ну не знаю я ченовскую библиотеку вдоль и поперек перешерстил, и 16 выборок делал и 128... и получал соответственно 8 или 64 полосы в спектре... и по-своему все переделывал... результат одного из вариантов можно видеть тут вживую практицки

[rubic]

хорошо давай вместо ADC смоделируем синусоиду:

FFT_N=256//fft.h
uint8_t y=0;
for(double i=0;i<2*M_PI;i+=2*M_PI/FFT_N)
{
capture[y]=511+511*sin(F*i);//F от 1 до 128
y++;
}

fft_input(capture, bfly_buff);
fft_execute(bfly_buff);
fft_output(bfly_buff, spektrum);

протестируй на симуляторе с разными F,будут заполнятся в [spektrum] адреса больше 32?

[ARV]

я не буду заниматься ерундой. 256 точек в память atmega8/16 не влезет - я такой цели не ставил. а для 128 и ниже - уже проверял. но я переписал все функции по-своему - оставил только само FFT, а ввод и перестановку свои сделал. все отлично делается.

[Dx!]

протестируй на симуляторе с разными F,будут заполнятся в [spektrum] адреса больше 32?

Я тоже игрался с библиотекой чена - работает адекватно. Давал генератором синусоиду, получал один-единственный столб спектра. Работала до 1024 выборок, давала честные 512 отсчётов спектра, в реальном времини.
Единственно - там передискретизация может быть. Доходим до 8кгц (ацп меги с тактом 250кгц) и при дальнейшем повышении частоты источника спектр показывает её уменьшение... Нужно аналогово фильтровать...

[GG Markes]

Могу я поинтересоваться, чем уважаемые форумчане компилируют прогаммы чена?

[ARV]

Могу я поинтересоваться, чем уважаемые форумчане компилируют прогаммы чена?

avr-gcc, он же WinAVR

[GG Markes]

avr-gcc, он же WinAVR

Второй глупый вопрос. Не могли бы Вы пераметры командой строки показать? А то я запутался.

[Dx!]

Вот левой пяткой писанная тестовая программка, в том числе и сабж тестирует.
http://www.xdevs.com/kb/dx/tmp/!PRG.lcd_fft_test.rar

Если вопрос только в FFT - этого много. Но, думаю, разобраться возможно. Тест выводит спектр на LCD от нокии. На частотах до ~8.5Кгц всё адекватно.

[GG Markes]

Вот левой пяткой писанная тестовая программка, в том числе и сабж тестирует.
http://www.xdevs.com/kb/dx/tmp/!PRG.lcd_fft_test.rar

Если вопрос только в FFT - этого много. Но, думаю, разобраться возможно. Тест выводит спектр на LCD от нокии. На частотах до ~8.5Кгц всё адекватно.

Большое спасибо! Разобрался