ДПФ/БПФ/ДКП на Cortex-M3 (1986ВЕ94Т), Вычисление БПФ на ARM Опции

[uwboy]

Господа!
Есть необходимость решить задачу спектрального анализа квазипериодического сигнала. Частота дискретизации не менее 80 кГц (изначально 100 кГц, но возможно выбрать исходя из удобного числа отсчётов). Поступает сигнал непрерывно. Спектральный анализ проводится сначала приблизительный — на первой четверти выборки; затем точный — на всей выборке. Длина полной выборки 1 секунда. Для спектрального анализа следует выполнить одно из дискретных преобразований Фурье (как несведущий затрудняюсь определить).
Для ЦОС собираюсь использовать микросхему (1986ВЕ94Т, 1986ВЕ93Т или 1986ВЕ92Т). Предварительно планирую систему следующим образом. Буфер для исходных отсчётов собираюсь организовать на отдельных двух микросхемах (1645РУ4АУ): первая хранит первую четверть выборки и результаты ДПФ, вторая — остаток выборки. Перенаправлять вывод с внешнего АЦП собираюсь внешней же логикой. Сразу после захвата первой четверти отсчётов МК следует начать ДПФ на первой четверти выборки. Закончить следует быстрее, чем 250 мс (соответствует времени захвата четверти выборки), при этом должно остаться немного времени на анализ спектра (поиск пиков). По завершении захвата остатка выборки, МК выполняет ДПФ на первой четверти и остатке менее, чем за секунду.

Вопросы такие:

1. Успеет МК выполнить ДПФ?
2. Какой случай ДПФ более предпочтителен (вход — действительные целые, выход — амплитуды целые/дробные)?
3. Значительно ли повысит быстродействие количество отсчётов являющееся степенью числа 2?
4. Возможно ли использовать внутренний АЦП микросхемы 1986ВЕ94Т и единственную микросхему СОЗУ и копировать из внутреннего АЦП во внешнее СОЗУ отсчёты по прерываниям АЦП во время вычисления ДПФ? Не сильно ли это притормозит основную программу (ДПФ)?
5. Значительно ли вырастет производительность ДПФ при использовании внутренней СОЗУ микросхемы 1986ВЕ94Т для части операций?

Поправьте, если я неправильно задаю вопрос или не в тот раздел обращаюсь. Поиск использовал, но просветления не достиг. Заранее спасибо, джентльмены!

[Golikov A.]

А вы с преобразование Фурье хорошо знакомы?

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

если вы передадите отсчеты подряд 0 1 2 3 4 5 6 7 , то значения по частотам получите в битово - реверсивном порядке, то есть
0 4 2 6 1 5 3 7, где в каждой ячейке будет стоять амплитуда спектра частоты выборки делить на 2 в степени индекс. как то так.. При этом в 0 ячейке лежит постоянка с кривой амплитудой.

Если хотите частоты нормальные, то надо так же загнуть первоначальные отсчеты. Именно поэтому в ДСП что рассчитаны на преобразование фурье есть спец модуль который железным образом делает такую индексацию.


Для ускорения спектрального анализа надо применять БПФ - быстрое преобразование фурье, самое распространенное это так называемая схема бабочка, она требует количество входных отсчетов равных степени 2, потому ответ да степень 2 серъезно ускорит процесс.

Потом есть еще проблема при вырезании четверти сигнала, на концах вы получите обрезанный спектр, то есть локализцете сигнал во времянной области, а принцип неопределенности Гизенберга гласит что при этом вы получите бесконечный спектр, другими словами вырезав часть сигнала и посчитав спектр вы получите спектр чего угодно, но не исходного сигнала.

Для устранения проблем используют оконное преобразование фурье, оно немного гадит амплитуду постоянки, но в целом это можно учесть. Весь входной сигнал надо умножить на функцию возрастающую от 0 до 1 и обратно спадающую, возрастание и падение плавное, по какой нить экспоненте.

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

[Aleksandr Barano...]

Господа!
Частота дискретизации не менее 80 кГц (изначально 100 кГц, но возможно выбрать исходя из удобного числа отсчётов). Поступает сигнал непрерывно. Спектральный анализ проводится сначала приблизительный — на первой четверти выборки; затем точный — на всей выборке. Длина полной выборки 1 секунда. Для спектрального анализа следует выполнить одно из дискретных преобразований Фурье (как несведущий затрудняюсь определить).

Получается 80000выброк.с*0.75с = 60000точек. Это как-то фантастично выглядит. Какое же Вам нужно частотное разрешение?

[uwboy]

А вы с преобразование Фурье хорошо знакомы?

Из библиотеки FFTW запускал. Её документацию читал, но не вчитывался. В общих чертах представляю, что такое ДПФ и что именно оно вычисляет, но могу вчитаться в книжку Лайонса за подробностями, правда сомневаюсь, что что-то принципиально новое узнаю.

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

Прошу прошения, но что имеется в виду?

Если хотите частоты нормальные, то надо так же загнуть первоначальные отсчеты. Именно поэтому в ДСП что рассчитаны на преобразование фурье есть спец модуль который железным образом делает такую индексацию.

Вот про такое был не в курсе.

ответ да степень 2 серъезно ускорит процесс.

Хорошо, спасибо. Попробую 65 или 131 кГц.

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

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

Вообщем в вашей концепции кроме выбора микросхем есть еще куча математических заморочек

Если это может послужить оправданием, то у меня это в первый раз.

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

FFTW с оптимизациями делает это на современном x86-м довольно долго, так что ARM захлебнётся в аналогичных вычислениях. Вы меня прямо опечалили.

Получается 80000выброк.с*0.75с = 60000точек. Это как-то фантастично выглядит. Какое же Вам нужно частотное разрешение?

Нет, для точный расчётов первая четверть тоже используется, т.е. отсчётов в большом преобразовании будет 80000 или 100000. Разрешение по частоте должно быть 1 Гц. Именно так.

целочисленных вариантов я что-то не припоминаю, ну кроме варианта ооочень больших целых типа бит на 128

А вот за этот вариант что скажете?

[Golikov A.]

скажу что нужен сопроцессор, ПЛИС какую нибудь, и на ней сделать параллельные вычислители, или специализированный считатель типа вот в соседней теме обсуждают
http://electronix.ru/forum/index.php?showtopic=114799
на этой платке стоит какой то могучий спец процессор который прям только считает

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

вот неплохой труд по БПФ
http://www.wl.unn.ru/ftp/public/analog/5.pdf


Оцифровывая с частотой 80 КГц, вы разглядите сигналы до 40 КГц, хорошо бы кстати на входе поставить фильтр на эту частоту, чтобы высшие частоты не портили сигнал. Но почему вы хотите для расчета использовать именно все точки из 1 секундного периода? Какая у вас минимальная частота в сигнале?
Считая дискретное преобразование вы получите частоты в зависимости от числа точек

40КГц 20КГц 10 5 2.5 1.25 ....

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


за вариант скажу вот что
DEST = ((long)(A) * (long)()>>15
то есть таблица синусов в 16 битных интах, а считают то в лонгах, и думаю бит на 64... АРМ без математического модуля не сильно быстрее лонги перемножает по сравнению с флотами. А если брать с модулем математическим, то может и флоты будет за 1 такт считать...

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

[Lmx2315]

Закончить следует быстрее, чем 250 мс (соответствует времени захвата четверти выборки), при этом должно остаться немного времени на анализ спектра (поиск пиков).

..если пиков немного то лучше отказаться от Фурье и считать Алгоритм Гёрцеля.

[Aleksandr Barano...]

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

Нет, для точный расчётов первая четверть тоже используется, т.е. отсчётов в большом преобразовании будет 80000 или 100000. Разрешение по частоте должно быть 1 Гц. Именно так.

Прямоугольное окно сильно испортит картину. Какой тогда смысл в таком высоком разрешении? А сигнал, случайно не узкополосный? Если да, то частоту дискретизации можно сильно уменьшить.
Чтобы оценить время выполнения, можно зарядить 65536 - точечное БПФ на какой-нибудь мощной PC и засечь время.

[Alex11]

Про степень 2, что тут говорили - это относится не к частоте дискретизации, а к количеству отсчетов в окне. Здесь от Вас уже все добиваются - что Вы хотите анализировать и найти: какая нижняя частота сигнала, которая Вас интересует? Если Вам не нужно анализировать сигналы порядка 1 Гц, то нет смысла считать Фурье такой длины. Если нужно проанализировать наличие сигнала в данной выборке, но более высокочастотного, то лучше сделать несколько преобразований по более коротким буферам с перекрытием и "сложить" результаты. Времени потратите существенно меньше. В любом случае, Вам нужно использовать Фурье с окном. Выбор окна зависит от задачи - какие у Вас пики и что требуется получить в результате анализа - разрешение по частоте или анализ амплитуд гармоник.

И еще. Вам нужно разрешить пики, отстоящие по частоте на 1 Гц, возможность определить частоту пика с точностью 1 Гц или у Вас нижняя частота интересующего Вас сигнала 1 Гц? Это существенно разные вещи. Длинная выборка строго необходима только в первом и третьем случаях. Да и то, 1 сек не хватит для 1 Гц. Если пики отстоят друг от друга достаточно далеко, то можно обойтись гораздо меньшей кровью и подсчитать положение пика дополнительной обработкой спектра после Фурье.

[DASM]

А вы с преобразование Фурье хорошо знакомы?

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

если вы передадите отсчеты подряд 0 1 2 3 4 5 6 7 , то значения по частотам получите в битово - реверсивном порядке, то есть
0 4 2 6 1 5 3 7, где в каждой ячейке будет стоять амплитуда спектра частоты выборки делить на 2 в степени индекс. как то так.. При этом в 0 ячейке лежит постоянка с кривой амплитудой.

Если хотите частоты нормальные, то надо так же загнуть первоначальные отсчеты. Именно поэтому в ДСП что рассчитаны на преобразование фурье есть спец модуль который железным образом делает такую индексацию.


Для ускорения спектрального анализа надо применять БПФ - быстрое преобразование фурье, самое распространенное это так называемая схема бабочка, она требует количество входных отсчетов равных степени 2, потому ответ да степень 2 серъезно ускорит процесс.

Потом есть еще проблема при вырезании четверти сигнала, на концах вы получите обрезанный спектр, то есть локализцете сигнал во времянной области, а принцип неопределенности Гизенберга гласит что при этом вы получите бесконечный спектр, другими словами вырезав часть сигнала и посчитав спектр вы получите спектр чего угодно, но не исходного сигнала.

Для устранения проблем используют оконное преобразование фурье, оно немного гадит амплитуду постоянки, но в целом это можно учесть. Весь входной сигнал надо умножить на функцию возрастающую от 0 до 1 и обратно спадающую, возрастание и падение плавное, по какой нить экспоненте.

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

Странно, зачем же тогда Техас библиотеки под 64 дсп писал почему только флоат то?

[Golikov A.]

я не говорил только, но исходил из того что дабл на арме считать сильно дольше чем флоат. Причем без мат модуля получалось что 2 лонга перемножить чуть ли не дольше чем 2 флота, а считать в 16-32 битных интах мне кажется не хватит разрядности. Могу ошибаться, это на уровне интуиции... в добавок автор хочет еще и массивы запихивать в 80К точек, думаю там любая фиксированная точка к концу расчета переполнится....

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

Кстати АРМ с мат модулем даблы за сколько перемножает?

[uwboy]

..если пиков немного то лучше отказаться от Фурье и считать Алгоритм Гёрцеля.

Надо наперёд знать частоту. Мне наперёд частоты пиков неизвестны, даже приблизительно. А вот комбинировать с неточным БПФ вполне здраво: сначала вычислить приблизительное положение пиков через БПФ на четверти выборки, а затем вычислить найти точные частоты пиков через Алгоритм Гёрцеля.

Прямоугольное окно сильно испортит картину.

Лишние регулярные составляющие будут?

Какой тогда смысл в таком высоком разрешении?

Частоту измерять.

А сигнал, случайно не узкополосный?

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

Если да, то частоту дискретизации можно сильно уменьшить.

Возможно, верхняя частота исследуемого сигнала может лежать в диапазоне от пары десятков герц, до 15 кГц. А вот оценить частоту нужно с точностью до 1 Гц.

Чтобы оценить время выполнения, можно зарядить 65536 - точечное БПФ на какой-нибудь мощной PC и засечь время.

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

Про степень 2, что тут говорили - это относится не к частоте дискретизации, а к количеству отсчетов в окне.

Ну, да. А у меня длина маленького окна 250 мс, а большого 1 секунда.

Здесь от Вас уже все добиваются - что Вы хотите анализировать и найти: какая нижняя частота сигнала, которая Вас интересует?

Допустим, 21 герц, но нужно знать, что это не 20 Гц и не 22 Гц, с вероятностью 50% и более.

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

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

возможность определить частоту пика с точностью 1 Гц

Именно это. Причём нужно это как можно быстрее. Темп выдачи составляет одну секунду. Начальное запаздывание на секунду мне простительно.

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

Допустим, что соседние пики отстоять могут на 3 Гц. Кстати, буду очень признателен, если расскажете мне о «дополнительной обработке спектра после Фурье».

Сообщение отредактировал uwboy - Aug 20 2013, 07:17

[prig]

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

А Вы сами хорошо знакомы?

И откуда требования к точности?
Целочисленные(с фиксированной точкой) 32-битные варианты, включая варианты с групповой плавающей точкой, более чем достаточны для подавляющего числа приложений. Во всяком случае, для 64К БПФ с 32-битной групповой плавающей точкой надо ещё поискать тот сигнал, что бы точность вычисления начала сказываться.

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

[Golikov A.]

А Вы сами хорошо знакомы?

И откуда требования к точности?
Целочисленные(с фиксированной точкой) 32-битные варианты, включая варианты с групповой плавающей точкой, более чем достаточны для подавляющего числа приложений. Во всяком случае, для 64К БПФ с 32-битной групповой плавающей точкой надо ещё поискать тот сигнал, что бы точность вычисления начала сказываться.

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

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

[uwboy]

Спасибо, господа, за ваши советы.
Уже прочёл про полифазное БПФ.
Покрутил частотный план всего измерителя. Предел, который достижим при допустимой погрешности следующий:

• частота дискретизации 32768 Гц (соответственно четверть выборки 8192 отсчёта, целое, 12 или 16 разрядов, плюс полная выборка в 32768 отсчётов);
• до шести узкополосных составляющих (каждая состоит из несущей и паразитной 10-20% АМ в пределах от 2 до 200 Гц); несущие располагаются случайно в диапазоне от 100 Гц до 16 кГц; разброс мощностей составляющих до 10 дБ (или ±5 дБ).

Требуется определить частоты несущих с точностью до 1 Гц.
Какие есть у вас рекомендации?

Сообщение отредактировал uwboy - Aug 20 2013, 12:05

[DASM]

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

Вы цену и производительность фиксед пойнт и дабла сравнили?

[_pv]

http://www.embeddedsignals.com/ARM.htm
там бенчмарки есть, для длины 4096 - сотня тактов на отсчёт.
32 bit FFT increases dynamic range by 90 dB , needs extra 20% to 50% cycles

[Golikov A.]

Вы цену и производительность фиксед пойнт и дабла сравнили?

нет.

[MSP430F]

Всем доброго времени суток!
Не балует гугл ответами на запрос о 1986ВЕ94Т. Эта тема - одна из трех (!) на весь инет ссылок! Миландр рулит!
У меня аналогичная задача как у топикстартера, и то же на этом МК. Только требования поскромнее. Надо посчитать БПФ для 4096 точек за 1 секунду. Чувствую, что МК справится. Хочу считать на float (32 бита), double не нужен, АЦП у меня 16 бит. В этой работе www.rf.unn.ru/NATO/1ws/SfP_Perov.pdf‎ люди посчитали, что для АЦП с разрядность от 11 до 23 float достаточно.

1986ВЕ94Т - это Cortex-M3, есть встроенный умножитель за 1 такт 32x32. Ответьте, пожалуйста, кто знает, сколько тактов занимает в таком МК умножение двух float ? Судя по http://habrahabr.ru/sandbox/58129/ даже целочисленные вычисления зависят от компилятора, поэтому если есть несколько результатов, приведите все. Спасибо.

Извиняюсь, правильная ссылка вот

[Golikov A.]

как мне тут указали несколькими постами выше, что и 32 бит с фиксированной точкой должно хватить, и умножение у вас будет делаться за 1 такт, может попробовать так?

[MSP430F]

как мне тут указали несколькими постами выше, что и 32 бит с фиксированной точкой должно хватить, и умножение у вас будет делаться за 1 такт, может попробовать так?

Я не понимаю, как можно реализовать целочисленное БПФ с 32 разрядными числами и получить хоть какую-нить приемлемую точность?

[Golikov A.]

32 инт это
от +- 2 147 483 648
соответственно можно считать что это

+- 21474.83648

или
+- 214748.3648

в целом правда разрядов дофига... в той статье что вы приводили они использовали 16 битный инт, и посчитали его пригодным для рассчета 11 бит АЦП, а он +-32768. С 32 битным получается на порядок больше целых, и еще 4 десятичных знака.. Главное не потерять точность в вычислениях.

[prig]

Я не понимаю, как можно реализовать целочисленное БПФ с 32 разрядными числами и получить хоть какую-нить приемлемую точность?

Вообще-то, 32 разряда дают динамический диапазон, который в реальной практике не встречается.
Ну, потеряются при БПФ несколько битов от 32-х, как это скажется на результате? Да никак. С 16-ю битами тоже не всё плохо.

В своё время приходилось делать оценку 16-бит версий БПФ с групповой плавающей точкой для 8К. Потеря точности (оценка по шуму) - примерно 1...2 бита.
При нормировании данных к 15-ти битам, шум от 8К 16-бит радикс-2 БПФ сопоставим с шумом 14-бит оцифровки.
Т.е. SNR преобразования что-то около 86дБ для шумоподобного сигнала. А вот 4-й радикс шумней на 6дБ из-за дополнительной потери 1-го бита при постадийном нормировании.

В том, что касается сабжа, на перемножителе 32*32 можно сделать вполне приемлимый вариант улучшенного 16-битного(16 данные+16 защита) БПФ на 64К с фикс. точкой. Лучше 90 дБ по шуму вполне достижимо. Потерь точности из-за нормирования данных ни на 2-м ни на 4-м радиксе не будет.

А вот со временем выполнения всё несколько сложнее. Эффективность кода по сравнению с сигнальниками будет существенно ниже.
8К 16-бит радикс-2 БПФ с групповой плавающей точкой на 533МГц блэкфине, со всем что можно в L1-памяти, выполняется за 1мс.
Грубая прикидка для 64К на 80МГц М3 даёт что-то около 350мс (1*10*7*5, где 5 - опять же, грубая оценка эффективности к блэкфину).

Так что, выбор 1986ВЕ94Е для данной задачи ТС несколько сомнителен.
Но что бы сказать наверняка, можно попробовать найти готовый радикс-4 и проверить на реальной железке или симуляторе.

Если искать что-нибудь другое, следует учитывать большой объём данных. Если не размещать данные и коэффициенты в памяти/кэше первого уровня, производительность может очень сильно упасть. На том же 533МГц блэкфине для 64к будет использоваться внешняя SDRAM. Результат будет на уровне 100мс вместо ожидаемых 10мс.

...Надо посчитать БПФ для 4096 точек за 1 секунду. Чувствую, что МК справится. Хочу считать на float (32 бита), double не нужен, АЦП у меня 16 бит. В этой работе www.rf.unn.ru/NATO/1ws/SfP_Perov.pdf‎ люди посчитали, что для АЦП с разрядность от 11 до 23 float достаточно.

Фикс точка на 32 бита (20 данные + 12 защита) Вам более чем хватит. 1986ВЕ94Е сгодится.
float (32 бита) - 24 значащих бита и экспонента. Это явно избыточно и слишком затратно. Проще разобраться с фикс. точкой.

Сообщение отредактировал prig - Aug 21 2013, 17:23

[_pv]

А вот со временем выполнения всё несколько сложнее. Эффективность кода по сравнению с сигнальниками будет существенно ниже.8К 16-бит радикс-2 БПФ с групповой плавающей точкой на 533МГц блэкфине, со всем что можно в L1-памяти, выполняется за 1мс. Грубая прикидка для 64К на 80МГц М3 даёт что-то около 350мс (1*10*7*5, где 5 - опять же, грубая оценка эффективности к блэкфину).

533000 тактов на 8к отсчётов это 65тактов на отсчёт,
ftp://ftp.analog.com/pub/www/technology/d...kfin_FftDct.zip :

Код

Performance     :
                Code Size   : 504 Bytes.
                Cycle Count : 3 * N * M   +   20 * M   -   N   +   17
                          where N = FFT length and M = log(N) to the base 4.
                      137 cycles for FFT size of   16.
                      589 cycles for FFT size of   64.
                     2913 cycles for FFT size of  256.
                    14453 cycles for FFT size of 1024.

соответственно 1507505 для длины 65536 или 23 такта на отсчёт, только вот памяти под коэффициенты не напасёшся.
в ссылке выше fft для M3 получалось как раз в 5 раз больше около сотни тактов на отсчёт при 4096.

[MSP430F]

Вообще-то, 32 разряда дают динамический диапазон, который в реальной практике не встречается.
Ну, потеряются при БПФ несколько битов от 32-х, как это скажется на результате? Да никак. С 16-ю битами тоже не всё плохо.

Фикс точка на 32 бита (20 данные + 12 защита) Вам более чем хватит. 1986ВЕ94Е сгодится.
float (32 бита) - 24 значащих бита и экспонента. Это явно избыточно и слишком затратно. Проще разобраться с фикс. точкой.

Вот здесь человек пишет буквально следующее
" Перед автором стояла задача написать преобразование Фурье на 2048 точек при разрядности исходных данных 16 бит. Из-за отсутствия арифметического сопроцессора пришлось делать целочисленное преобразование, что создало некоторые трудности. При разрядности исходных данных 16 бит разрядность коэффициентов должна быть не менее 16, чтобы не происходило потери точности. Их произведение содержит 32 разряда. 2048 точек дают еще 11 дополнительных разрядов, а это значит, что в 32-разрядное процессорное слово промежуточные данные не помещаются. Вычисление каждой “бабочки” ведется с точностью 64 разряда, а результат округляется до 32 разрядов. "

Как Вы тогда можете это прокомментировать ?

[_pv]

Вот здесь
Как Вы тогда можете это прокомментировать ?

были данные 16 разрядов сделали преобразование Фурье, стало вдруг 32 разряда, если теперь обратное сделать чтобы исходные данные получить для результата уже 64 бита понадобится?
кто мешал по ходу округлять чтобы из 32 бит не вылезать.

[Golikov A.]

Вот здесь человек пишет буквально следующее
" Перед автором стояла задача написать преобразование Фурье на 2048 точек при разрядности исходных данных 16 бит. Из-за отсутствия арифметического сопроцессора пришлось делать целочисленное преобразование, что создало некоторые трудности. При разрядности исходных данных 16 бит разрядность коэффициентов должна быть не менее 16, чтобы не происходило потери точности. Их произведение содержит 32 разряда. 2048 точек дают еще 11 дополнительных разрядов, а это значит, что в 32-разрядное процессорное слово промежуточные данные не помещаются. Вычисление каждой “бабочки” ведется с точностью 64 разряда, а результат округляется до 32 разрядов. "

Как Вы тогда можете это прокомментировать ?

это математическая оценка. Которая во всех случаях даст 100% точность без потери
в общем случае перемножение двух 16 битных чисел дает 32 битное число, если вы перемножите два максимальных числа, но если вы перемножите два числа в серединке, то получится бит 20-24...
этим и пользуются практики.

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


Кстати у меня вопрос. Когда делаете 16 битное фурье на тучу точек, вычисления ведут в 20-32 битах, но время от времени округляют до 16? Потому что просто прям суммировать 16 бит точно переполнится, даже при 11 битах входных...
А процедура округления и нормировки не создает лишних тактов больше чем если бы сразу во флотах считать?

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

[MSP430F]

были данные 16 разрядов сделали преобразование Фурье, стало вдруг 32 разряда, если теперь обратное сделать чтобы исходные данные получить для результата уже 64 бита понадобится?
кто мешал по ходу округлять чтобы из 32 бит не вылезать.

Я не понял, здесь люди что-нибудь в теме понимают, или, как часто на форумах бывает, дилетанты от нечего делать язык чешут. Чуток хоть надо в топике понимать.
Итак, имеем БПФ на 4096 точек. Используем целочисленную арифметику. Это 12 этапов (log2 4096) умножений в общем случае комплексных (это сейчас не принципиально) чисел. Если АЦП на 16 разрядов, то, чтобы не потерять точность при умножениях, коэффициенты должны быть не хуже по разрядности, пусть тоже будет 16 бит. Первый этап. Умножаем два таких числа. В худшем случае получается 32 битное число. А впереди еще 11 этапов. Чтобы не потерять точность при округлениях, результат надо разделить на 2 в 15, и мы получили промежуточный результат уже 17-ти разрядный. Умножаем его на коэффициент, который опять-таки 16-разрядный и получаем 33-разярядное число. К концу этапов, коих 12, имеем 44 разряда.

Если округлять на каждом этапе до 32 бит, это будет эквивалентно в итоге использованию 4-разрядных коэффициентов ! Круто, да ? И кому будет нужна такая точность ?

[prig]

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

...
В худшем случае получается 32 битное число. А впереди еще 11 этапов. Чтобы не потерять точность при округлениях, результат надо разделить на 2 в 15, и мы получили промежуточный результат уже 17-ти разрядный. Умножаем его на коэффициент, который опять-таки 16-разрядный и получаем 33-разярядное число. К концу этапов, коих 12, имеем 44 разряда.

Если округлять на каждом этапе до 32 бит, это будет эквивалентно в итоге использованию 4-разрядных коэффициентов ! Круто, да ? И кому будет нужна такая точность ?

Ну и каша у Вас в голове. Результат рассуждений соответствующий.
Попробуйте разобратся с тем, что такое фикс. точка, и как её используют.

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

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

...
fft для M3 получалось как раз в 5 раз больше около сотни тактов на отсчёт при 4096.

Я примерно так и оценивал, но это действительно грубая оценка.
Кто знает, сколько придётся обрабатывать эту заготовку напильником, и что из этого получится?

[Golikov A.]

Я не понял, здесь люди что-нибудь в теме понимают, или, как часто на форумах бывает, дилетанты от нечего делать язык чешут. Чуток хоть надо в топике понимать.
Итак, имеем БПФ на 4096 точек. Используем целочисленную арифметику. Это 12 этапов (log2 4096) умножений в общем случае комплексных (это сейчас не принципиально) чисел. Если АЦП на 16 разрядов, то, чтобы не потерять точность при умножениях, коэффициенты должны быть не хуже по разрядности, пусть тоже будет 16 бит. Первый этап. Умножаем два таких числа. В худшем случае получается 32 битное число. А впереди еще 11 этапов. Чтобы не потерять точность при округлениях, результат надо разделить на 2 в 15, и мы получили промежуточный результат уже 17-ти разрядный. Умножаем его на коэффициент, который опять-таки 16-разрядный и получаем 33-разярядное число. К концу этапов, коих 12, имеем 44 разряда.

Если округлять на каждом этапе до 32 бит, это будет эквивалентно в итоге использованию 4-разрядных коэффициентов ! Круто, да ? И кому будет нужна такая точность ?

да это верно, и не верно одновременно...

4 разрядные коэффициенты это
0 1 2 3

хотим посчитать
1/2 + 1/2 + 1/2 + 1/2 = 0.5 + 0.5 + 0.5 + 0.5= 2;
в этих терминах 4 разрядных чисел
0.5 == 1
1/2 + 1/2 + 1/2 + 1/2 = 0.5 + 0.5 +0.5 + 0.5 = 1 + 1 + 1 + 1 = 4== 0;



если мы возьмем 3 разрядный
0 1 2 3 4 5 6 7
и будем считать что у нас младший разряд после точки, то получим
0 0.5 1 1.5 2 2.5 3 3.5

берем наш пример
1/2 + 1/2 + 1/2 + 1/2
переводим его в 3 разрядные числа
1/2 + 1/2 + 1/2 + 1/2 = 0.5 + 0.5 +0.5 + 0.5 = 1 + 1 + 1 + 1 = 4
переводим результат обратно в 2 разрядные
4 == 2 .

посчитано точно, причем в начале и конце 2 разрядные числа в середине 3 разрядные.

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


Теперь что происходит с флоат.

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

то есть это как
2 + 4 + 8 + 4 == 2*(1+2+4+2)

именно этим и пользуются при счете в целых, как я понимаю....

[Ruslan1]

У меня аналогичная задача как у топикстартера, и то же на этом МК. Только требования поскромнее. Надо посчитать БПФ для 4096 точек за 1 секунду. Чувствую, что МК справится.

У меня PIC32 (ядро MIPS32, 40MHz) справлялся с FFT8192 за примерно 530 мс (оптимизация кода на минимуме, это отладка была). Все расчеты делаются в floaf, причем плавающей точки в процессоре нет (все программно, стандартными функциями компилятора). Это совершенно без ухищрений и допиливаний, сишные исходники FFT тоже стандартные из интернета. Скорость устраивает- значит не лезу улучшать.

Уверен, что ARM должен быть шустрее.

[MSP430F]

Кто-нибудь может привести проверенный рабочий код на C (можно ссылку) целочисленного БПФ, использующего 32-разрядную (не выше) арифметику для 16-разрядных чисел ?

И чтобы без потери точности при вычислениях.

Сообщение отредактировал MSP430F - Aug 22 2013, 11:16

[uwboy]

сишные исходники FFT тоже стандартные из интернета

Прямо с Википедии, или адресок подскажете?

[Aleksandr Barano...]

Кто-нибудь может привести проверенный рабочий код на C (можно ссылку) целочисленного БПФ, использующего 32-разрядную (не выше) арифметику для 16-разрядных чисел ?

И чтобы без потери точности при вычислениях.

Попробуйте поискать гуглом с ключевым словом fix-fft.

[prig]

И чтобы без потери точности при вычислениях.

Вообще без потери точности? Такого не найдёте, т.к. с точки зрения реальных задач это полный абсурд.

Сообщение отредактировал prig - Aug 22 2013, 13:19

[MSP430F]

Вообще без потери точности? Такого не найдёте, т.к. с точки зрения реальных задач это полный абсурд.

Я имею ввиду, что, если на входе даные с разрядностью 16 бит, то хотелось бы и на выходе иметь данные с той же и разрядностью и точностью.

Все что я смог найти, лишь на 1024 точки. Думаю, что все-таки сделать БПФ на 4096 точки с честными 16 разрядами на целочисленной 32-битной орифметике невозможно в принципе, потому и прошу, покажите мне такой код, а то больше смахивает на троллинг.

Сообщение отредактировал MSP430F - Aug 22 2013, 14:01

[Ruslan1]

Прямо с Википедии, или адресок подскажете?

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

Кстати, что любопытно: я даже от таблицы отказался, синусы-косинусы на ходу считаю. Сначала была таблица, но снес (большая и, как я говорил уже, рекорды скорости не ставил).
Сравнивал алгоритм с результатами Матлабовской функции: разница не превышает рассчетную из-за разной битности данных.

CODE

#ifndef FFT_CALC_H
#define FFT_CALC_H
typedef f32 TREAL;
//#define M_PI 3.1415926535
extern void FFT(u16 n,TREAL *x,TREAL *y);
#endif //FFT_CALC_H

void FFT(u16 N,TREAL *x,TREAL *y)
{
f32 c,s,t1,t2,t3,t4,u1,u2,u3;
int i,j,p,l,L,M,M1,K;
L=N;
M=N/2;
M1=N-1;
while(L>=2)
{
l=L/2;
u1=1.;
u2=0.;
t1=M_PI/(f32)l;
c=cos(t1);
s=(-1)*1*sinf(t1);
// j = FFT_N_SAMPLES/L;
// c=Sinewave[(j + FFT_N_SAMPLES/4)&(FFT_N_SAMPLES-1)]; //c=cos(t1);
// s=(-1) * Sinewave[j]; //s=(-1)*1*sin(t1);
for(j=0; j<l;j++)
{
for(i=j;i<N;i+=L)
{
p=i+l;
t1=*(x+i)+*(x+p);
t2=*(y+i)+*(y+p);
t3=*(x+i)-*(x+p);
t4=*(y+i)-*(y+p);
*(x+p)=t3*u1-t4*u2;
*(y+p)=t4*u1+t3*u2;
*(x+i)=t1; *(y+i)=t2;
}
u3=u1*c-u2*s;
u2=u2*c+u1*s; u1=u3;
}
L/=2;
}
j=0;
for(i=0;i<M1;i++)
{
if(i>j)
{
t1=*(x+j); t2=*(y+j);
*(x+j)=*(x+i); *(y+j)=*(y+i);
*(x+i)=t1; *(y+i)=t2;
}
K=M;
while(j >=K)
{
j-=K;K/=2;
}
j+=K;
}
}

[prig]

Я имею ввиду, что, если на входе даные с разрядностью 16 бит, то хотелось бы и на выходе иметь данные с той же и разрядностью и точностью.

Все что я смог найти, лишь на 1024 точки. Думаю, что все-таки сделать БПФ на 4096 точки с честными 16 разрядами на целочисленной 32-битной орифметике невозможно в принципе, потому и прошу, покажите мне такой код, а то больше смахивает на троллинг.

Любой код с фикс. точкой 1.31 (нотация по ADI; с форматом и тем, как с ним работать, предлагаю разобраться самому).
Размещаете данные в младших 19-ти битах и знаковом разряде (не забыв правильно заполнить биты, куда не попали исходные данные). Коэффициенты 32 бита.
После преобразования данные нормализуются в нужный формат.

20 бит - это 122дБ. ДД преобразования будет примерно на уровне 110дБ.
Это эквиваленто снижению SNR исходного идеального 16-бит сигнала(98дБ) на величину менее 1 дБ.
Для сигнала на уровне 90дБ о привнесённом шуме можно будет вообще забыть.
Вы собираетесь работать со значениями SNR 98дБ? Сильно сомневаюсь.

Это моё понимание вопроса.
Могу несколько ошибаться в значениях (таки, лучше промоделировать для конкретной системы). Опустил некоторые детали.
Подход достаточно стандартный. Работает на 100%. Проверено на реальном железе.

А вот что Вы понимаете под "той же точностью" и чего хотите добиться, вопрос тёмный.

ПС Естественно, для 64К надо будет использовать 16 бит из 32 или прикручивать дополнительную нормализацию между стадиями, что и есть групповая плавающая точка.

Сообщение отредактировал prig - Aug 22 2013, 16:47

[Alex11]

По поводу дополнительной обработки спектра после Фурье. 1. Перед расчетом нужно использовать окно, отличное от прямоугольного. В зависимости от выбранной функции можно оптимизировать различные параметры. Мы использовали окно Гаусса, т.к. требовалось получить максимально точные значения амплитуд. В любом случае форма линии в спектре повторяет форму оконной функции. Если пик одиночный, то можно расчитать его положение по 3 - 5 точкам у вершины гораздо точнее одного бина, используя метод наименьших квадратов. Вы получаете параметры функции - положение и амплитуду. В
Вашем случае, когда может быть непредсказуемая модуляция, она может вносить ошибку, но, полагаю, незначительную по частоте. В случае наложения двух пиков в пределах ширины ответ будет произвольным и разрешить эти линии не удасться.

[Golikov A.]

Любой код с фикс. точкой 1.31 (нотация по ADI; с форматом и тем, как с ним работать, предлагаю разобраться самому).
Размещаете данные в младших 19-ти битах и знаковом разряде (не забыв правильно заполнить биты, куда не попали исходные данные). Коэффициенты 32 бита.
После преобразования данные нормализуются в нужный формат.

ПС Естественно, для 64К надо будет использовать 16 бит из 32

Вы наверное очень в теме, а мы совсем не в теме. Мне на данный момент БПФ с фикс точкой не нужно, но для общего развития интересно. Не могли бы вы немного спуститься до нашего уровня и разъяснить что есть:

фикс.точка 1.31, да еще и в нотации ADI.
какие данные надо разместить в младших 19 битах, и чем правильно заполнить биты? Какой формат нужный?
И сколько бит вы предлагали использовать для 4К точек, 19? тогда для 64К почему 16?

Вы написали много текста, но пропустив некоторые для вас "очевидные" вещи, сделали для нас простых смертных его полностью не читаемым...

хотя может я один такой, не понимающий.

[MSP430F]

Любой код с фикс. точкой 1.31 (нотация по ADI; с форматом и тем, как с ним работать, предлагаю разобраться самому).
Размещаете данные в младших 19-ти битах и знаковом разряде (не забыв правильно заполнить биты, куда не попали исходные данные). Коэффициенты 32 бита.

А я о чем говорю? Если коэффиценты у Вас 32 бита (очень хорошие, кстати, коэффициенты), то как Вы собираетесь реализовывать арифметику в 32 разрядах ? Я же не спорю о том, что целочисленные вычисления - это круто. Я только говорю, что если у меня данные честные 16-битные и я хочу получить на выходе спектр с диапазоном тоже в 90 дБ, то целочисленных вычислений в 32 разрядов для этого не достаточно. Для этого-то человек ЗДЕСЬ и использует вычисления с точностью 64 разряда.

[prig]

...
Не могли бы вы немного спуститься до нашего уровня и разъяснить что есть...
...
фикс.точка 1.31, да еще и в нотации ADI.
какие данные надо разместить в младших 19 битах, и чем правильно заполнить биты? Какой формат нужный?
И сколько бит вы предлагали использовать для 4К точек, 19? тогда для 64К почему 16?

Вы написали много текста, но пропустив некоторые для вас "очевидные" вещи, сделали для нас простых смертных его полностью не читаемым...

хотя может я один такой, не понимающий.

Знаете, сильно лень. И так мого букв, и так много времени потрачено. Разве что, из чисто альтруистических соображений.
Так что, считайте это только подсказками, когда будете разбираться.

- ADI я упомянул сознательно, и именно потому, что у них всё это хозяйство с фикс. точкой описано достаточно хорошо и подробно.
Полистайте их доки на сигнальники, статьи и книжки по DSP. У TI всё то же самое, но описано несколько хуже.
А в длинных и умных книжках могут так запутать простые вещи, что без пол-банки...

Если вкратце и грубо:
1.31 - запись вещественных чисел со знаком (фактически - целое число). Диапазон типа +-0.9999... Умножение в этом формате (фактически, это целочисленное со знаком 32*32 с дополнительным сдвигом влево на один бит и округлением результата до 32-битного числа) даёт результат в этом же формате. Т.е., если умножается 0.0999.. на 0.9999.., должно получаться 0.0999... С точностью до младшего бита, есстественно.

- Коммент с заполнением полей при переводе формата и нормировке для БПФ - не самый удачный.
Корректней описать так:
16 бит целое со знаком (или же 1.15, что в общем-то одно и то же) копируется в старшие 16 бит формата 1.31, а младшие 16 бит заполняются нулями. После чего делается арифметический сдвиг вправо всего 32-бит числа для выделения "запасных" битов).

- Результат вычисления бабочки радикс-2 может вырасти на один бит (в два раза) по сравнению с исходными значениями. Если использовать 1.31 без постадийного нормирования, то для 4К приходится оставлять 12 "запасных" бит, а для 64К - 16. Т.о. получается, что в первом случае округление будет вестись по фактическим 20-ти разрядам, а во втором - по 16-ти.

Если использовать постадийное нормирование (групповую плавающую точку), требуется 1 "запасной" бит для радикс-2 и 2 бита для радикс-4. Производительность упадёт, но всё-равно будет лучше, чем для 32-битного числа с плавающей точкой. А точность выше оного.

П.С. Поймите меня правильно. Я не пытаюсь устроить здесь ликбез или что-то кому-то показать. Время это занимает много, толку для меня самого почти никакого.
Так что, ещё раз повторюсь, считайте это только подсказками для желающих изучить вопрос. Не более.

А я о чем говорю?
...

О, теперь понял, о чем это Вы.
Оказывается, Вы не в курсе, что при целочисленном умножении 32*32 получается 64-битный результат, который может благополучно округляться.

[Golikov A.]

Знаете, сильно лень. И так мого букв, и так много времени потрачено.

Пишите много, но не емко...
и того.

- ADI я упомянул сознательно, и именно потому, что у них всё это хозяйство с фикс. точкой описано достаточно хорошо и подробно.
Полистайте их доки на сигнальники, статьи и книжки по DSP. У TI всё то же самое, но описано несколько хуже.
А в длинных и умных книжках могут так запутать простые вещи, что без пол-банки...

ADI означает фирму, и надо обращаться к их программным кодам

Если вкратце и грубо:
1.31 - запись вещественных чисел со знаком (фактически - целое число). Диапазон типа +-0.9999... Умножение в этом формате (фактически, это целочисленное со знаком 32*32 с дополнительным сдвигом влево на один бит и округлением результата до 32-битного числа) даёт результат в этом же формате. Т.е., если умножается 0.0999.. на 0.9999.., должно получаться 0.0999... С точностью до младшего бита, есстественно.

1.31 означает 1 бит знака, 31 бит - значащие

- Коммент с заполнением полей при переводе формата и нормировке для БПФ - не самый удачный.
Корректней описать так:
16 бит целое со знаком (или же 1.15, что в общем-то одно и то же) копируется в старшие 16 бит формата 1.31, а младшие 16 бит заполняются нулями. После чего делается арифметический сдвиг вправо всего 32-бит числа для выделения "запасных" битов).

наверное все таки влево, иначе вы просто обнулите число или я чего то не понимаю? Да и вообще сдвигая 32 бита на 32 бита вы получите либо 0 либо 64 бита, тут хотелось бы определенности. Потому что если 0, то это глупо, а если 64, то откуда у нас 32 битная арифметика?

- Результат вычисления бабочки радикс-2 может вырасти на один бит (в два раза) по сравнению с исходными значениями. Если использовать 1.31 без постадийного нормирования, то для 4К приходится оставлять 12 "запасных" бит, а для 64К - 16. Т.о. получается, что в первом случае округление будет вестись по фактическим 20-ти разрядам, а во втором - по 16-ти.

Почему опять у вас 4К БПФ стало 20 разрядным, а 64К - 16 разрядным? Или вы считаете число отбрасываемых разрядов?

П.С. Поймите меня правильно. Я не пытаюсь устроить здесь ликбез или что-то кому-то показать. Время это занимает много, толку для меня самого почти никакого.
Так что, ещё раз повторюсь, считайте это только подсказками для желающих изучить вопрос. Не более.

Ну вас трудно понять, потому и времени занимает много. Были бы ближе к народу, было бы быстрее, может в этом для вас был бы толк? Много раз убеждался объясняя самые простые вещи всегда и для себя что-то новое понимаешь...


П.С. кстати все время упускаются процы с арифм модулем которые флоты за 1 такт перемножают, а таких не мало...

А я о чем говорю? Если коэффиценты у Вас 32 бита (очень хорошие, кстати, коэффициенты), то как Вы собираетесь реализовывать арифметику в 32 разрядах ? Я же не спорю о том, что целочисленные вычисления - это круто. Я только говорю, что если у меня данные честные 16-битные и я хочу получить на выходе спектр с диапазоном тоже в 90 дБ, то целочисленных вычислений в 32 разрядов для этого не достаточно. Для этого-то человек ЗДЕСЬ и использует вычисления с точностью 64 разряда.

естественно если у вас 32 битные операнды, то все их виды без потерь перемножить можно только в 64 битном числе

Но если у вас коэффициенты 32 битные, но со значениями от -2^15 .. +2^15, то у вас они фактически 16 битные, и их можно перемножить в 32 битах без потерь. Игра идет на том что реальный сигнал и реальные коэффициенты никогда не занимают весь диапазон иначе бы не помогли ни флоты ни даблы.

если у вас 32 бита то вы можете иметь 2^32 разных чисел, и не важно какой формат, во флоте их тоже 2^32, потому что просто комбинаций их 32 бит больше нет, потому считать что во флоте, что в целых 32 битных точность на выходе одинаковая, просто для флота она получится сама, а для 32 целых надо принять некоторый набор телодвижений.

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

[Maverick]

Кто-нибудь может привести проверенный рабочий код на C (можно ссылку) целочисленного БПФ, использующего 32-разрядную (не выше) арифметику для 16-разрядных чисел ?

И чтобы без потери точности при вычислениях.

посмотри это

[prig]

-Обращаться к программным кодам ADI не обязательно. Достаточно разобраться с с форматами и арифметикой.

- 1.31(нотация ADI) означает вещественное, 1 бит знака, 31 бит - дробная часть в дополнительном коде

- Для правильной работы БПФ в формате с фиксированной точкой 1.X и произвольными данными, исходные данные арифметически сдвигаются вправо на необходимое количество "защитных" бит N, где N - двоичный логарифм размерности БПФ.

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

Сообщение отредактировал prig - Aug 23 2013, 11:37

[MSP430F]

посмотри это

Так там на ASM, к тому же 32-бита не под Cortex M3, а только для ARM 9E.

-Обращаться к программным кодам ADI не обязательно. Достаточно разобраться с с форматами и арифметикой.

- 1.31(нотация ADI) означает вещественное, 1 бит знака, 31 бит - дробная часть в дополнительном коде

- Для правильной работы БПФ в формате с фиксированной точкой 1.X и произвольными данными, исходные данные арифметически сдвигаются вправо на необходимое количество "защитных" бит N, где N - двоичный логарифм размерности БПФ.

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

То есть для 16-разрядных данных в 4096 точках "фактическая" разрядность БПФ = 16 - 12 = 4. Так по Вашему ?

[prig]

То есть для 16-разрядных данных в 4096 точках "фактическая" разрядность БПФ = 16 - 12 = 4. Так по Вашему ?

Не исходных данных.
Для 32-х разрядной арифметики с фикс. точкой (напоминаю про 32*32 и проч.) "фактическая" разрядность БПФ = 32 - 12 = 20.
Шум преобразования, приведённый ко входу, приблизительно 6*(20-2)+2дБ.

Сообщение отредактировал prig - Aug 23 2013, 11:35

[MSP430F]

Не исходных данных.
Для 32-х разрядной арифметики с фикс. точкой (напоминаю про 32*32 и проч.) "фактическая" разрядность БПФ = 32 - 12 = 20.
Шум преобразования, приведённый ко входу, приблизительно 6*(20-2)+2дБ.

Вкусно звучит! Может код приведете для примера такого чудо-алгоритма ? (Это без иронии).

[Golikov A.]

Наконец то понял что куда сдвигается...

то есть если у нас есть считалка 32 битная
и мы хотим сделать на ней 4К фурье, то нам надо в нее пихать числа не более 20 бит.
а для 64К - 16 бит.

Да я чет как-то тупанул...

Бабочка то никогда не даст больше чем в 2 раза от входа (я помню вы это говорили, но я думал с учетом вида сигнала), Y = X1 + X2 * w
w, X1 и X2 в пределах [-1 : 1], так что ответ 2 в худшем случае...

А весь путь что пройдет сигнал это как раза Log2(N) бабочей, так что больше чем 2^N запас не нужен...

[Maverick]

Так там на ASM, к тому же 32-бита не под Cortex M3, а только для ARM 9E.

а это?
и это?

как я понимаю это должно отвечать Вашим требованиям

Я имею ввиду, что, если на входе даные с разрядностью 16 бит, то хотелось бы и на выходе иметь данные с той же и разрядностью и точностью.

Все что я смог найти, лишь на 1024 точки. Думаю, что все-таки сделать БПФ на 4096 точки с честными 16 разрядами на целочисленной 32-битной орифметике невозможно в принципе, потому и прошу, покажите мне такой код, а то больше смахивает на троллинг.

но я сильно не вникал в код...

[MSP430F]

а это?
и это?

как я понимаю это должно отвечать Вашим требованиям

но я сильно не вникал в код...

К сожалению, в аннотации к эти кодам сказано
"16 bit FFT precision comparable with other fixed point implementation – precision determined by necessary scaling by 0.5 in every FFT stage",
то есть точность НИКАКАЯ ! А "32 bit FFT increases dynamic range by 90 dB , needs extra 20% to 50% cycles" - только на ARM 9E.
Да мне бы просто код на С без ассемблера...

[MSP430F]

Если поворочивающие коэффициенты 16-разрядные, то на каждом этапе БПФ идет накопление ошибки 1/32768 = 0,00003. За 12 этапов (для 4096 точек) получим уже 0,00037 (= 12 * 0,00003), что эквивалентно использованию 12-битного АЦП вместо 16-битного!
Чтобы обеспечить достойную точность, разрядность коэффициентов должна быть больше 16 разрядов. В этом случае результат умножения превысит 32 разряда, что я и пытаюсь доказать. Вот здесь, для примера, человек пишет "32 bit FFT increases dynamic range by 90 dB". Посмотрите код. Там используются 32-разрядные коэффициенты и специфическое умножение SMULL (Sign Long Multiply) с 64-разрядным результатом.

[Golikov A.]

это если вы 16 бит умножая на 16 бит, пытаетесь на 16 битах остаться.

имея умножитель 32 на 32 бита, вы можете за 1 такт считать в 32 битах, для 4К фурье можно считать в 32-12 = 20 битах величинах, и потеряв 4 бита, останетесь в своих 16 битах...

в 16 битах считать естественно смысла нет...

[Maverick]

еще это посмотрите
и это
Мне кажется Вам нужно будет писать функции для работы с большими числами, иначе высокую точность не получите при использовании целочисленной арифметики (мое мнение)

[MSP430F]

это если вы 16 бит умножая на 16 бит, пытаетесь на 16 битах остаться.

имея умножитель 32 на 32 бита, вы можете за 1 такт считать в 32 битах, для 4К фурье можно считать в 32-12 = 20 битах величинах, и потеряв 4 бита, останетесь в своих 16 битах...

в 16 битах считать естественно смысла нет...

И я о том же. Однако, prig, как я понял, сватал мне именно 64-разрядную арифметику, называя ее почему-то 32-х разрядной. Я было даже повелся , думая человек какую-то военную, не иначе, тайну знает. Да, я и говорю, совет экспертов тут, а не форум электронщиков.

[Golikov A.]

И я о том же. Однако, prig, как я понял, сватал мне именно 64-разрядную арифметику, называя ее почему-то 32-х разрядной. Я было даже повелся , думая человек какую-то военную, не иначе, тайну знает. Да, я и говорю, совет экспертов тут, а не форум электронщиков.

он просто очень плохо объясняет... он предлагал 16 битный вход, и 32 битная арифметика, и
коэффиценты считать в 20 битах для 4К фурье

[prig]

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

- Умножение формата 1.31*1.31 в стандартном "С" реализуется только с использованием целочисленного умножения формата 64*64. Большого смысла в такой реализации нет по причине низкой эффективности. Лучше поискать поддержку фикс. точки (для сигнальников - без проблем) или "прикрутить" её самому.
- Если компилятор "не понимает" формат и "не умеет" выполнять операции с фикс. точкой, простейший и наиболее эффективный путь реализации умножения - ассемблерный код.
Для умножения формата 1.31*1.31 с результатом в формате 1.31 используется целочисленное умножение в формате 32*32 с результатом в формате 64 (на ассемблере - без проблем). Промежуточный 64-битный целочисленный результат нормализуется (сдвиг точки и округление) в 32-битный результат с фикс. точкой.
- Даже если компилятор поддерживает типы с фикс. точкой, рекомендуется реализовать на ассемблере, с использованием промежуточных 64-бит, хотя бы "бабочку" БПФ.

П.С.
Способ реализации умножения 1.31*1.31 с результатом в формате 1.31 не отменяет того, что это 32-битная арифметика.
То же самое касается 32-битного формата с плавающей точкой.

он просто очень плохо объясняет... он предлагал 16 битный вход, и 32 битная арифметика, и
коэффиценты считать в 20 битах для 4К фурье

М.б., плохо объясняю. М.б., кое-кто невнимательно читает или не слишком точен в том, что пишет. Предлагал я несколько иное.

Сообщение отредактировал prig - Aug 30 2013, 11:38

[Golikov A.]

Ну или опять кто-то не понятно объяснил
Я имел в виду фикс точка нормированная в 20 бит, промежуточные операции 32 битные...

[prig]

Ну или опять кто-то не понятно объяснил
Я имел в виду фикс точка нормированная в 20 бит, промежуточные операции 32 битные...

Ну да, исходные данные, нормированные в 20 бит формата с фикс. точкой. Скажем так. Для 4К, естественно.

Сообщение отредактировал prig - Aug 30 2013, 11:36