FFT на Си для AVR

[INA]

Спасибо Вам за быстрые ответы и вообще за оперативное общение!

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

1. Каким образом влияет частота дискретизации на результаты обработки данных БПФ?
2. Как мне рассчитать таблицы синуса и косинуса для своего (какого-то...) приложения, если мне, например вздумалось выделить из спектра три любые частоты (скажем, 1000Гц, 1200Гц, 1500Гц)?
С помощью какого-нибудь стороннего софта я могу это сделать? Слышал, что это можно сделать в Экселе...
Есть МатЛаб - быть может он может помочь в генерации таблиц?

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

P.S. Да, проект Чана я видел... Люди говорят, что очень оптимальный вариант БПФ, но мне пока нужно разобраться в сути...
Не хочется просто в тупую использовать чужие библиотеки, не понимая как всё это работает...

[pcb]

2. Как мне рассчитать таблицы синуса и косинуса для своего (какого-то...) приложения, если мне, например вздумалось выделить из спектра три любые частоты (скажем, 1000Гц, 1200Гц, 1500Гц)?

Целиком задача как звучит?

1. Каким образом влияет частота дискретизации на результаты обработки данных БПФ?

Чем выше частота дискретизации, тем меньше шаг по частотной линии.

[YS]

Спасибо Вам за быстрые ответы и вообще за оперативное общение!

Не за что, просто я сегодня весь день за компом.

1. Каким образом влияет частота дискретизации на результаты обработки данных БПФ?

Сдвигает максимальную частоту анализа и меняет разрешение по частоте.

2. Как мне рассчитать таблицы синуса и косинуса для своего (какого-то...) приложения, если мне, например вздумалось выделить из спектра три любые частоты (скажем, 1000Гц, 1200Гц, 1500Гц)?

ДПФ не выделяет отдельные частоты, оно выделяет все разом. Нельзя так посчитать ДПФ, чтобы получить только, например, три результата - это будет уже не ДПФ, а несколько вычислений скалярного произведения/коэффициента корреляции. Как я уже говорил, отдельная выборка ДПФ это результат рассчета корреляции сигнала с отдельной частотой. Вещественная часть - с косинусом, мнимая - с синусом.

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

мне пока нужно разобраться в сути...

Суть следующая:

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

Если есть сигналы [1,2,3] и [4,5,6] их скалярное произведение будет 1*4 + 2*5 + 3*6 = 32.

2. Когда нам надо построить спектр (то есть проверить кучу частот), можно действовать как в п.1, но это будет долго. Такое вычисление "в лоб" (или почти "в лоб") называется дискретным преобразованием Фурье. Оно дает некоторое приближение спектра сигнала.

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

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

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

[INA]

Целиком задача как звучит?

Приветствую, Вас!

Да целиком задача звучит достаточно банально - понять как работает БПФ, применительно к AVR и пр. железкам...
Просто хочу своими руками всё это дело ощутить, чтоб не пользоваться сторонними библиотеками, не понимая как всё работает.

Да, про частоту дискретизации я понял, что чем она выше, тем меньше шаг.
Я, наверное не правильно задал вопрос...
Меня вот какой момент интересует: допустим частота дискретизации сигнала с АЦП 10000Гц, тогда для реализации алгоритма определения амплитуды нескольких частот из этого спектра (пусть те же - 1000, 1200, 1500 Гц) мне нужны таблицы синусов и косинусов (для ускорения работы алгоритма)...
Так вот, на этом моменте вычисления этих таблиц (или массивов) я должен учитывать то, что мои данные с АЦП идут с частотой дискретизации 10000Гц?
Т.е., другими словами, при вычислении самих таблиц я должен оперировать этой (10000Гц) частотой?
Ведь иначе, как я понимаю данные не будут совпадать?

[YS]

понять как работает БПФ, применительно к AVR и пр. железкам...

Поймите главное - преобразование Фурье это преобразование Фурье, все остальное (расчет таблиц, операция-"бабочка") - всего лишь методы оптимизации, чтобы вычислять исходное преобразование быстрее. Конкретно про алгоритм Кули-Тьюки неплохо написано тут, но я бы не забивал этим голову; оптимизация - задача специалистов в области вычислительной математики, а не инженеров. Инженеру важно понимать, что преобразование Фурье дает на выходе и как это применять. Например, с помощью БПФ можно лихо формировать сигналы с OFDM-модуляцией.

Так вот, на этом моменте вычисления этих таблиц (или массивов) я должен учитывать то, что мои данные с АЦП идут с частотой дискретизации 10000Гц?

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

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

В общем, это уже зависит от выбранной стратегии оптимизации.

А собственно сгенерировать таблицу можно в любой среде. Хоть на том же C, хоть на Python, Lua и пр., да и в Excel тоже.

[INA]

Сдвигает максимальную частоту анализа и меняет разрешение по частоте.

Т.е., для работы БПФ совершенно без разницы, с какой частотой исходный сигнал был оцифрован?

Суть следующая:

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

Если есть сигналы [1,2,3] и [4,5,6] их скалярное произведение будет 1*4 + 2*5 + 3*6 = 32.

2. Когда нам надо построить спектр (то есть проверить кучу частот), можно действовать как в п.1, но это будет долго. Такое вычисление "в лоб" (или почти "в лоб") называется дискретным преобразованием Фурье. Оно дает некоторое приближение спектра сигнала.

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

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

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

Вот тут я что-то совсем заблудился...
Если при выделении просто нескольких частот достаточно корреляции сигнала с таблицами синуса и косинуса (ну и дальнейшего нахождения амплитуды),
то при БПФ что за таблицами пользуются?
Такой вопрос возник потому, что в исходниках многих проектов используются таблицы.
Одни явно описаны как таблицы синусов и косинусов для определённых частот, а другие содержат (как у Чана) просто данные разной длины...
И у него нет отрицательных значений в таблицах. (?)
Возможно, создатели проектов как раз и пользовались разными подходами к решению подобных задач?

[INA]

Вот тут
есть статья, в которой человек берёт именно таблицы синусов и косинусов.

Умножая их потом на 127, но это уже частности...
В тексте он пишет, что "...имеем 512 отсчетов АЦП нужно посчитать мнимую и действительную части для 150Гц при частоте дискретизации 19200 Гц:..."

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

[YS]

Т.е., для работы БПФ совершенно без разницы, с какой частотой исходный сигнал был оцифрован?

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

То есть, если есть буфер длиной 128 выборок, это будет 65-я выборка (нулевая не учитывается, в ней постоянная составляющая).

Если эти 128 выборок получены при частоте 44100 Гц, то 1-я будет соответствовать частоте примерно 344.5 Гц, вторая - 681.1 Гц, и так далее, до 65-й, которая будет соответствовать частоте 22050 Гц.

Если они получены при, например, 8 кГц, 1-я будет соответствовать 62.5 Гц, вторая - 125 Гц, и так далее, до 4 кГц (65-я выборка).

Если при выделении просто нескольких частот достаточно корреляции сигнала с таблицами синуса и косинуса

Не с таблицами, а вообще с синусом и косинусом. Синус и косинус можно вычислять рядами Тейлора, можно в уме, можно на сопроцессоре или еще как-то, а можно таблицами. Таблицы - просто один из способов ускорения. Это не есть обязательный элемент. Не зацикливайтесь на нем.

то при БПФ что за таблицами пользуются?

Зависит от конкретной оптимизации. В таблицы складывают то, что долго считать. А алгоритмов БПФ есть много, алгоритмом Кули-Тьюки список способов не кончается.

И у него нет отрицательных значений в таблицах. (?)

Зависит от представления данных.

Возможно, создатели проектов как раз и пользовались разными подходами к решению подобных задач?

Они по-разному оптимизируют одну и ту же операцию.

Умножая их потом на 127, но это уже частности...

Таблицы - это тоже частность. Поставить вычисление рядом Тейлора мешает только то, что это долго.

Но он не пишет каким образом он их сгенерировал... ?

Сгенерировать их можно на любом языке по очевидной формуле - s[k]=sin(2*п*p*(k/N)), где N - количество отсчетов, максимум для k, p - количество периодов на N отсчетов.

[INA]

Вот спасибо большое!
Теперь всё начинает вставать на свои места...
Мне нужен тайм-аут, чтобы всё вышеизложенное переварить...
Ну и поэкспериментировать с софтом...
Вопросы ещё появятся!...
Ещё раз большое спасибо!

[YS]

Не за что, разбирайтесь.

[INA]

По ходу раздумий ещё один вопрос появился...
Уже в практической плоскости...
И вопрос, похоже не один...

1. Насчёт дискретизации: получается, что если я точно не знаю частоту дискретизации (значения в массиве, считанные с АЦП уже есть), то я просто-напросто не смогу определить какой выборке соответствует какая частота?

2. Есть тот же самый массив данных с АЦП... Теперь, чтобы получить амплитуды гармоник мне нужно найти синус и косинус какой-то выборки и по формуле Пифагора A = sqrt(real[10]*real[10] + imag[10]*imag[10]); (10 - это для примера 10-й отсчёт какой-то частоты) найти её амплитуду?

И вот тут по ходу написания вопросов возник ещё один:
У меня есть выборка... скажем та же 10-я, те adc_data[10]... это, естественно число какой-то разрядности...
Мне нужно найти его синус и косинус, чтоб воспользоваться формулой Пифагора?
Никаких накоплений в виде суммы произведений как тут x_real[k] += x_n[n] * Math.Cos( (-2) * Math.PI * k * n / 1024 ); не нужно?

[YS]

1. Насчёт дискретизации: получается, что если я точно не знаю частоту дискретизации (значения в массиве, считанные с АЦП уже есть), то я просто-напросто не смогу определить какой выборке соответствует какая частота?

Ага. Массив он и есть массив.

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

Никаких накоплений в виде суммы произведений как тут x_real[k] += x_n[n] * Math.Cos( (-2) * Math.PI * k * n / 1024 ); не нужно?

Нужно. Сама суть - в накоплении, расчете корреляции.

[INA]

Нужно. Сама суть - в накоплении, расчете корреляции.

Тогда, если Вы не против, я Вас ещё немного задержу...
Вы уж простите мой тупизм...

В общем, я сделал так:
Есть схемка... простая на Меге32...
вход АЦП, которой смещён на +2,5В для того, чтобы получить среднюю точку и входные данные шли без искажений.
Теперь при считывании данных в массив из них сначала вычитается число 512, те. data[i] = adc_out - 512;
Таким образом в массиве появляются положительные и отрицательные числа, в соответствии со значением амплитуды в тот момент времени...
Теперь из них нужно получить действительную и мнимую части через формулу суммирования произведений...(???)
Вот тут я встрял...
real_data[i] += data[i] * sin(2 * pi * i * (вот тут застрял)); ???

[YS]

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

Для начала сгенерируем сигнал. Пускай это будет сумма двух синусоид, одна с относительной частотой 5 и немного сдвинутая по фазе, другая - 10 и без сдвига. Ну и шума добавим, куда без него...

Код:

signal = {};

F1=5;
P1=0.8;
F2=10;
k=0;
while k<1000 do
   signal[k]=math.cos(2*math.pi*F1*(k/1000)+P1) + math.cos(2*math.pi*F2*(k/1000)) + 1*(math.random() - 0.5);
   k=k+1;
end;


Выглядит это так:



А теперь попробуем обнаружить сигнал, считая корреляцию:

Код:

REL_FREQ=5;
k=0;
s_re=0;
s_im=0;
while k<1000 do
   s_re=s_re+math.cos(2*math.pi*REL_FREQ*(k/1000))*signal[k];
   s_im=s_im+math.sin(2*math.pi*REL_FREQ*(k/1000))*signal[k];
   k=k+1;
end;

print(s_re,";",s_im,";",math.sqrt(s_re^2+s_im^2));


REL_FREQ - та частота (относительная, разумеется), которую обнаруживаем. На выходе печатаем вещественную часть, мнимую часть, и амплитуду.

Вывод для REL_FREQ=5:

Код:

340.9868715185 ; -352.52017847724 ; 490.45134598816

Видим, что относительная амплитуда - более 490. Много. Кстати, арктангенс отношения мнимой и вещественной части даст нам относительную фазу. atan(-352/340)=-0.8, что по модулю равно P1.

Вывод для REL_FREQ=10:

Код:

487.23896012134 ; 0.57229975367557 ; 487.23929622634

Относительная амплитуда тоже около 490, тоже много.

А теперь попробуем найти какую-нибудь частоту, которой нет в нашей смеси. Например, примем REL_FREQ=15:

Код:

-15.783519852489 ; 8.4533221429504 ; 17.904696428211

На выходе получили всего около восемнадцати. Если бы не было шума, получили бы ноль. Результат с убранным шумом (в первом фрагменте написал 0*(math.random() - 0.5)):

Код:

-2.9354296771089e-013 ; 2.8305136012818e-013 ; 4.0778124817472e-013

Везде минус тринадцатая степень. Ну, то есть, в сущности ноль. То есть, этой частоты в сигнале нет.

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

[INA]

Спасибо Вам ещё раз!
Огромное!

Просто через язык программирования до меня доходит быстрее, чем через математику...

В данном случае (в данном примере) частоты взяты относительные...
Но, по сути частоту можно указать явно?
В Герцах?

[YS]

В данном случае (в данном примере) частоты взяты относительные...
Но, по сути частоту можно указать явно?
В Герцах?

Если f - относительная частота, то частоту в Герцах можно выразить как

(f*F)/N, где F - частота дискретизации, N - количество анализируемых отсчетов (размер буфера).

В самом деле, если на N отсчетов приходится f периодов колебания (смысл относительной частоты в данном примере именно таков), то на один период будет приходиться N/f отсчетов.

С другой стороны, отсчеты берутся через промежутки времени 1/F секунд.

Тогда реальный период колебания составляет N/(f*F) секунд, ну и частота, соответственно, (f*F)/N.

[INA]

Ок!
Спасибо за терпение и подробные разъяснения!

Теперь попытаюсь всё это "вложить" в железо и поиграться наяву...


P.S. На редкость хороший форум и хорошие люди здесь!
А котов и кошек у меня в доме пятеро (!!!)
Трое взрослых и два котёнка...

[YS]

Не за что, не за что. Docendo discimus.

Да, о дружелюбности этого форума многие говорят.

А у меня живности в доме нет.

[INA]

Добрый день всем!

YS (прошу прощения, не знаю имени...), спасибо ещё раз за помощь!

Сегодня пробовал (да и вчера тоже) все полученные знания в железо воплотить...
Вроде как всё правильно сделал, но в Протеусе не хочет работать...
Думаю, что из-за Протеуса... Завтра на работе соберу платку и буду тестить в железе - так надёжнее...

[WiseLord]

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