Посоветуйте алгоритм (точнее его реализацию на С или С++) для ARM9 (NXP LPC3250) с ОЗУ 64МБ.
Один рабочий нашел вот тут. "Прицепил" к своему проекту. Все работает, но до величины окна 65536. Больше - и программа виснет. Под "кучу" выделил 15МБ ОЗУ. В принципе если постараться, то можно еще 40МБ добавить.

Интересует вопрос - каковы требования к объему памяти разных алгоритмов?




Тему можно удалить. Проблема не в БПФ, а во мне массив для получения результата не того размера подготовил

Хотя нет тему не стоит удалять. В общем у меня получилось примерно так.

тип double (8 байт) на реальную часть и тип double (8 байт) на мнимую часть. Если величина окна 2^20=1048576, то получается на исходный массив выборок - 16МБ, плюс 8МБ на массив WStore (пока не вник, не знаю точно можно от него избавится или нет). В итоге примерно 24МБ - вот такие требования. Мне пока не хватает примерно 4МБ, поэтому окно 2^19.

Считает это дело мое устройство примерно секунд 50 (100МГц тактовая, 64МБ 100МГц SDRAM). Есть у кого алгоритмы получше?



Где то вроде видел что-то такое: сделать много FFT последовательно накладывая окна, а потом эти результаты подвергнуть еще раз какому-то преобразованию. В итоге будет большое окно с маленькими требованиями к памяти. Такое возможно?

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

нет. есть оцифрованный сигнал. надо увидеть спектр (не в динамике, а именно график) на большом временном промежутке. примерно 1 сек. При этом частоты могут быть от 100Гц до 1МГц, примерно. Вот получается, что окно должно быть 2*10^6 выборки.

Вот сейчас в голову пришла такая "идея". Если разрешение по частоте допустимо порядка 100Гц, то брать последовательно по шкале времени участки (N штук) по 16384 элементов, вычислить для каждого FFT, просуммировать коэффициенты одинаковых частот и поделить на количество участков N. Получим спектральную характеристику продолжительного участка, но с небольшими требованиями к памяти и небольшим размером файла результата.

Как идея?

Одно "НО". Как суммировать коэффициенты? векторно (как комплексные числа) или по модулю?

По сути это вариант "скользящего" окна, но с суммированием.

Идея называется метод периодограмм.
Складывать надо квадраты модулей, т.е. мощности одноименных бинов N участков. Результат - квадратный корень из суммы, деленной на N.

Проверил метод на практике. Не понравилось. На плате есть источник шума 125кГц (DC-DC). При амлитуде входного сигнала 3 мВ и фиксированной частоте он раз в 100 меньше на спектрограмме, чем этот входной сигнал. Но если входной сигнал имеет широкополосную ЧМ модуляцию, то в итоге, из-за своей стабильной частоты, на периодограмме сигнал наводки увеличивается пропорционально количеству просуммированных окон. Получается что хоть полезный сигнал по амплитуде и гораздо выше, но преобладающим в периодограмме будет стационарный сигнал наводки с DC-DC. Интересная особенность, где-то может и полезная, а в моем случае это не то, что нужно.

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

может попробовать считать арифметикой с фиксированной точкой(исключить double и float, следя за переполнениями и вовремя сдвигая результат)?

А 6-байтовых real или float в этом компиляторе не существует. Порой удавалось на Single (4-х байтовых real) БПФ на 65000 точек заделать. Считал мгоновенно (нетбук). Алгоритм древний из Голда-Рабинера.Синусы и косинусы там интерактивно считались. Когда таблично пробовал задать, считало медленнее, как ни странно. Еще на электронике60 работало.

Хочется готовое 100% работающее решение. Имел ввиду что-то, что кардинально изменит требования к памяти и скорость, а не использование других типов. Ну и чтобы кто-нибудь поделился своим мнением и опытом по поводу различных реализаций БПФ.

А float наверно можно использовать.

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

Если же складывать когерентно, но с пропусками, то шаг между бинами Фурье конечно уменьшится и будет определяться максимальным расстоянием между отсчетами. Но частотная функция окна (преобразование Фурье от временной) будет очень широко размазана и убьет всё разрешение. Т.е. значения соседних бинов будут очень сильно размазаны.



Если использование длинного окна обусловлено требованиями к разрешению, кардинально ничего Вы не измените.
Скорость немного можно подправить использованием целых типов, а требование по памяти немного можно уменьшить как сказано выше, если считать синусы/косинусы итеративно, а не держать таблицу. Короткий целый тип тоже даст экономию в двое.
Кроме того можно найти реализации без всяких больших рабочих массивов WSTORE, в том числе и "на месте", т.е. такие чтобы входной массив совпадал с выходным. Теретически можно избавиться от всех таблиц (синусных и инверсной адресации) и провести преобразование "на месте". Тогда получите N/2 от ОЗУ


PS. Есть один такой чувак который лет 20 уже пишет книгу и библиотек про быстрые преобразования. Там посмотрите, может найдете "на месте". Там их столько, что огого. http://www.jjj.de/fxt/

А сколько нужно? 49 секунд устроит?

И какая нужна точность по амплитуде на выходе?

Для целочисленных алгоритмов число "шумящих" разрядов растет примерно как 0.5*[log2(N)-2], так что для FFT на 65536 точек с фиксированной запятой Вы потеряете семь младших бит.

Точнее можно промоделировать в MATLAB'е.

PS. Кстати, для вещественных данных: Efficient_FFT_Computation_of_Real_Input

нет. 49 тоже многовато до 5 секунд было бы хорошо.
АЦП 16 бит - поэтому очень жалко 7 бит терять. Или я неправильно понял?




Нет. Требование к разрешению примерно 100 Гц. 16384 элемента вполне хватает. Требование по времени, которое охватывает БПФ. Т.е. спектр за 1 секунду, например. Итого - 2млн элементов. Разрешение 1Гц мне не надо.

За ссылку спасибо

Всё правильно. Если процессор целочисленный, можно посчитать FFT в целых числах - в формате "fract 24.8".

Нашел в интернете интересную статью как раз о том, какие переменные при какой разрядности входных данных надо использовать. Мне при 16 битах надо float

Еще можно за 1 проход делать 2 действительных БПФ, но пока я не совсем понял, как это реализовать. Кому интересно - вот статейка




кстати ссылка о том же но на англицком

апну тему старую

я в своё время использовал fixed point для оценки спектра НЧ сигнала.Код красиво ложился на AT91 и DSP TI. Всё красиво получалось. На фоне этого статья которая говорит, что использовать int недопустимо, выглядит некомпетентно....
или я чего-то не понял

Интовое БПФ - это до сих пор дисертабельная тема.