Простые вопросы по FFT, Реализация на микроконтроллерах без плавающей точки. Опции

[JohnKorsh]

Добрый день | вечер | утро! Не поможет ли кто понять простые вопросы по FFT. Пытаюсь сделать FFT для MSP430. Для начала написал на C простое FT, затем, FFT. Сравнил результаты. Совпадают. Пытаюсь перевести в MSP430 без переменных типа float. (Побыстрее, думаю, будет). Не могу понять простую вещь: чтобы использовать поворотные множители создал массив. Пока он на C в виде float (cos и sin (2*Pi*k/N)) - всё нормально, но если попытаться умножить на константу (для того, чтобы перевести cos и sin в значения int - больше единицы) FT продолжает работать исправно (с ростом амплитуды на выходе), а FFT - распадается - появляются новые гармоники. Наверное так не должно быть, ведь FFT это лишь способ считать, или я не прав?

[fontp]

Добрый день | вечер | утро! Не поможет ли кто понять простые вопросы по FFT. Пытаюсь сделать FFT для MSP430. Для начала написал на C простое FT, затем, FFT. Сравнил результаты. Совпадают. Пытаюсь перевести в MSP430 без переменных типа float. (Побыстрее, думаю, будет). Не могу понять простую вещь: чтобы использовать поворотные множители создал массив. Пока он на C в виде float (cos и sin (2*Pi*k/N)) - всё нормально, но если попытаться умножить на константу (для того, чтобы перевести cos и sin в значения int - больше единицы) FT продолжает работать исправно (с ростом амплитуды на выходе), а FFT - распадается - появляются новые гармоники. Наверное так не должно быть, ведь FFT это лишь способ считать, или я не прав?

Алгоритм один, но ошибки накопляются по разному. Особенно если коэффициенты плохо округляются и имеют постоянное смещение

[Alex11]

Опять же переполнение надо смотреть.

[DRUID3]

...результат каждого умножения сдвигать на "число разрядов" - 1. Можно определить тип данных(fix point т.н.1.15

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

) где сдвиг будет автоматический. Но сэкономить "мипсу" оно может только тогда когда есть аппаратная поддержка этого формата - например АЦП выдает данные уже так...

[alex_os]

Добрый день | вечер | утро! Не поможет ли кто понять простые вопросы по FFT. Пытаюсь сделать FFT для MSP430. Для начала написал на C простое FT, затем, FFT. Сравнил результаты. Совпадают. Пытаюсь перевести в MSP430 без переменных типа float. (Побыстрее, думаю, будет). Не могу понять простую вещь: чтобы использовать поворотные множители создал массив. Пока он на C в виде float (cos и sin (2*Pi*k/N)) - всё нормально, но если попытаться умножить на константу (для того, чтобы перевести cos и sin в значения int - больше единицы) FT продолжает работать исправно (с ростом амплитуды на выходе), а FFT - распадается - появляются новые гармоники. Наверное так не должно быть, ведь FFT это лишь способ считать, или я не прав?

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

[JohnKorsh]

Пока непонятно. Для уверенности в правильности реализации FFT вывел в файл результаты FT и FFT и программно сравнил (для разных частот 1024 точки) - совпадают. Но наткнулся на ещё один непонятый эффект - при использовании окна (любого, отличного от прямоугольного) для FT всё работает, ка и должно - "размывается" центральный пик (немного), но хорошо подавляются все боковые частоты. А для той же ситуации в случае FFT появляются дополнительные гармоники. Видимо, что-то недопонимаю. Прилагаю программку с помощью которой патаюсь понять. Первые две закладки - симуляция, вторые две - сигнал с SoundBlaster-а, последняя закладка - не относится к FFT.

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

 S_Ft.rar ( 382.25 килобайт ) Кол-во скачиваний: 112

 

[DRUID3]

...это что за "троян"?

[JohnKorsh]

Извините, не проверял, может, где и "подцепил". Хотя, Касперский молчит.

[DRUID3]

Та не... это я так ...просто у меня и Wine нет, как я его посмотрю? Давайте скриншоты тогда уж...

[JohnKorsh]

Вот скриншоты. Моё понимание пока на том же уровне.

[JohnKorsh]

Вот скриншоты. Моё понимание пока на том же уровне.

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

 Furrie.doc ( 256 килобайт ) Кол-во скачиваний: 185

 

[SPACUM]

Вот скриншоты. Не слишком ли сильно "отщипнул" по краям?

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

[JohnKorsh]

Добрый день!
Ещё вопрос. Как я говорил - при попытке выполнить преобразование Фурье на микроконтроллере с фиксированной точкой я закладывал поворачивающие множители в виде:
W_re (k,N) = 32768 * cos (- 2 * Pi * k/N)
W_im (k,N) = 32768 * sin (- 2 * Pi * k/N)
Множитель 32768 – для того, чтобы sin и cos представить в виде, большем 1.
При таком подходе приходится делить (ну, понятно, просто сдвигать) часть результата в “бабочке” на 32768 (15 разрядов “вправо”), иначе FFT “разваливается” (простое преобразование Фурье не чувствительно к умножению = там просто общий результат возрастает в 32768 раз):
A_re = A_re + (W_re * B_re – W_im * B_im)/32768
A_im = A_im + (W_re * B_im + W_im * B_re)/32768
B_re = A_re - (W_re * B_re – W_im * B_im)/32768
B_im = A_im - (W_re * B_im + W_im * B_re)/32768
(Потом, понятно, идёт масштабирование результата, чтобы не было переполнения по алгоритму “блочной плавающей точки”.)
Наверное, есть какое-то общеизвестное более изящное решение. В Интернете пока не нашёл, может кто посоветует статью или исходник?


"Отщипнул" - я имел в виду, что огносительно большая часть энергии сигнала отбрасывается при домножении на оконную функцию. (Представленное на рисунке окно Хемминга, пожалуй, самый "щадящий" вариант окна). Может, по Вашему опыту, я не так её реализовал - должна бть более "плоская" вершина?

[fontp]

Добрый день!
Ещё вопрос. Как я говорил - при попытке выполнить преобразование Фурье на микроконтроллере с фиксированной точкой я закладывал поворачивающие множители в виде:
W_re (k,N) = 32768 * cos (- 2 * Pi * k/N)
W_im (k,N) = 32768 * sin (- 2 * Pi * k/N)
Множитель 32768 – для того, чтобы sin и cos представить в виде, большем 1.
При таком подходе приходится делить (ну, понятно, просто сдвигать) часть результата в “бабочке” на 32768 (15 разрядов “вправо”), иначе FFT “разваливается” (простое преобразование Фурье не чувствительно к умножению = там просто общий результат возрастает в 32768 раз):
A_re = A_re + (W_re * B_re – W_im * B_im)/32768
A_im = A_im + (W_re * B_im + W_im * B_re)/32768
B_re = A_re - (W_re * B_re – W_im * B_im)/32768
B_im = A_im - (W_re * B_im + W_im * B_re)/32768
(Потом, понятно, идёт масштабирование результата, чтобы не было переполнения по алгоритму “блочной плавающей точки”.)
Наверное, есть какое-то общеизвестное более изящное решение. В Интернете пока не нашёл, может кто посоветует статью или исходник?

А это ничего что
32768 * sin (- 2 * Pi * (N/4)/N) = 32768 * cos (- 2 * Pi * 0/N) = -1.0 вместо 1.0?

Попробуйте что нибудь типа
W_re (k,N) = floor(32767 * cos (- 2 * Pi * k/N) + 0.5)
W_im (k,N) = floor(32767 * sin (- 2 * Pi * k/N) + 0.5)
floor() это стандартная функция отсечения влево, соответственно floor(*+0/5) - округление

Блочная экспонента - это как раз самое изящное решение для вычислений с фиксированой точкой, не заморачивайтесь

[JohnKorsh]

Спасибо.