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

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

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

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

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

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

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

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

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

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

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

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

Добрый день!
Ещё вопрос. Как я говорил - при попытке выполнить преобразование Фурье на микроконтроллере с фиксированной точкой я закладывал поворачивающие множители в виде:
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) - округление

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

Спасибо.

Добрый день!
Вот очередной работающий вариант программы исследования преобразования Фурье.
Может, кому пригодится. Если понадобится исходник - напишите. Не кладу, просто, чтобы не загромождать сайт (CBuilder6).
Пока остался непонятым вопрос о соотношении сигнал/шум. При увеличении числа точек
соотношение сигнал/шум растёт не как корень из полосы сигнала к пjлосе разрешения преобразования, а гораздо медленнее.
Ошибки предыдущих вариантов были в неправильном использовании целочисленной арифметики. (Ошибки округления и, даже, переполнения.)
FFT написано так, как его видит MSP430 - 16-битная целочисленная арифметика, входной сигнал - 12 бит.
Кто заинтересуется - первые две закладки - имитация сигнала.
Вторые две закладки - работа с реальным сигналом со входа Sound Blaster-a.
Пятая закладка - не по делу - дальнейшие исследования.
Амплитуда должна быть не менее 200 милливольт. (Посмотрите на закладке, как рисуется сигнал).
Диапазон 3500-4500 Гц (По условиям моей задачи). При работе с реальным сигналом при увеличении
частоты максимальная составляющая спектра ползёт в сторону уменьшения - эффект от децимации.
Программа выдаёт файл отчёта и файл с перестановочным массивом и поворотными
множителями для выбранного количества точек.

Нужно будет свое тоже опубликовать... я 3 года назад тоже много чего понаписывал...

Добрый день! Выкладываю программку, немного облегчающую написание FFT. Программка считает таблицу перестановок и поворотные множители при различныхь параметрах. Для провекрки просчитавется FFT для полученных значений. Результат запасается в файл в директории Мои документы/W_Test. Результирующие массивы можно прямо из файла в исходный текст программы на C вставлять. Написана на Del[hi. Кто захочет улучшить - пишите, вышлю исходный код.