Преобразование Фурье для AVR

Больше
12 года 5 мес. назад - 12 года 5 мес. назад #1 от Gudd-Head
Gudd-Head создал эту тему: Преобразование Фурье для AVR
Хочется чего-нибудь такого замуть своего, умного...

Итак, если звёзды сойдутся, я когда-нибудь доберусь до момента нахождения спектра сигнала, преобразованного 8-ми битным АЦП, используя силы RISC МК. И вот что я надумал, основные идеи:
1. использовать синусно-косинусное преобразование
2. отбрасывать всё лишнее, стараться чтобы переменные по размеру не вылезали за 8 бит.

Теория.
Пусть у нас есть последовательность s[n], состоящая из Nfft=2^m 8-ми битных отсчётов входного сигнала.
Тогда мы можем представить каждый отсчёт последовательности s[n] в виде суммы Nfft/2 слагаемых вида
a(k)*sin(2*pi*k*n/Nfft)+b(k)*cos(2*pi*k*n/Nfft), (1)
где n=1...Nfft - номер отсчёта, k=0...(Nfft/2-1) - номер гармоники.
Таким образом, зная a(k) и b(k) мы можем найти амплитуду k-ой гармоники по формуле
sqrt(a(k)^2+b(k)^2), (2)
т.е. вычислить составляющие спектра.
Сами коэффициенты вычисляются так:
a(k)=sum(n=1...n=Nfft)[a(k)+s(n)*sin(2*pi*k*(n-1))/Nfft)] (3)
b(k)=sum(n=1...n=Nfft)[a(k)+s(n)*cos(2*pi*k*(n-1))/Nfft)] (4)
, причём k=0 - тривиальный случай: постоянная составляющая.
Реализация.
Естественно, каждый раз вычислять значение синуса или косинуса мы не будем, вобьём всё в таблицу. И с плавающей точкой работать мы тоже не будем, перейдём к целочисленной 8-ми битной арифметике чтобы было побыстрей.
Алгоритм, приготовления:
Считаем Nfft значений синуса через 2*pi/Nfft радиан, начиная с нуля, умножаем их на 127, округляем (чтобы в последствии использовать только аппаратный умножитель) и вбиваем в ОЗУ. Для Nfft=64 будет что-то вроде
[0, 12, 25, 37, ..., -49, -37, -25, -12] (5)
и лежать в диапазоне от -127 до +127, т.е. занимать 8 бит.
Вбиваем в ОЗУ Nfft значений нашего сигнала, предполагая его беззнаковым в диапазоне 0...255 (1 байт).
Очищаем (обнуляем) ещё Nfft байт ОЗУ.
Алгоритм, суть:
Заводим два счётчика, k от 1 до Nfft/2-1 и n от 1 до Nfft (пример на m-языке):
for k=1:(Nfft/2)
for n=1:Nfft
a(k)=a(k) + s(n)*sin(2*pi*k*n/Nfft);
b(k)=b(k) + s(n)*cos(2*pi*i*n/Nfft);
end;
end;
И получаем коэффициенты a(k) и b(k).

А вот самый фарш алгоритма с реализацией на АСМе на напишу позже (сейчас посто времени нет).
Последнее редактирование: 12 года 5 мес. назад от Gudd-Head.

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

Больше
12 года 5 мес. назад #2 от Gudd-Head
Gudd-Head ответил в теме Re: Преобразование Фурье для AVR
Фарш.

Поскольку готовой программы на АСМе у меня ещё нет, буду спускаться всё ниже и ниже по уровню языка программирования от человеческого до машинного, одновременно продумывая текст программы.
Итак, у нас есть сигнал [0...255]×Nfft, отсчёты синуса [-127...127]×Nfft и Nfft свободных байт ОЗУ.
Помимо этого понадобится несколько РОН и регистровых пар — косвенных указателей адреса.
Алгоритм основного цикла:
перемножаем отсчёт сигнала с соотв. значением синуса — получаем в итоге двухбайтное число со знаком [-32385...32385(7E81h)]. Отбрасываем младший байт, что эквивалентно целочисленному делению (с отбрасыванием дробной части) на 256, старший байт [-126...126] прибавляем к двухбайтному аккумулятору со знаком - 2 РОН. То же самое проделываем со значением косинуса - ещё 2 РОН.
После Nfft умножений и Nfft-1 сложений для синуса и такого же количества операций для косинуса, получаем в каждом аккумуляторе число [-Nfft*126...Nfft*126]. Т.е. необходимая разрядность аккумулятора m+8, 2^m=Nfft. Сохраняем в ОЗУ старшие 8 значимых бит результата. Проделываем это Nfft/2(-1) раз (что-то я сейчас не могу никак сообразить, одинаковые или разные получатся отсчёты Nfft/2 и Nfft/2-1).
Мдя, тяжеловато получилось: порядка Nfft^3. :( Но, будем надеяться, хотя бы компактно выйдет :)
Спустимся на один уровень.
Запускаем основной цикл: прогоняем k от 1 до Nfft/2(-1)
{
Запускаем вспомогательный цикл: n от 1 до Nfft
{{
Для умножения будем использовать команду MULSU для перемножения числа со знаком и числа без знака. Потому s(n) - n-ый отсчёт сигнала и sin(n·k) — n·k-ый отсчёт синуса загружаем в регистры номер 16...23, перемножаем. Старший байт результата будет в регистре R1. Складываем/вычитаем R1 с двухбайтным аккумулятором.
Пока у нас где-то в R16...R23 висит s(n), перемножим его с косинусом n·k. Для этого достаточно просто добавить 90° к аргументу синуса, т.е. Nfft/4! Т.о. подгружаем sin(n·k+Nfft/4) вместо sin(n·k), перемножаем, старший байт складываем/вычитаем с другим 16-битным аккумулятором.
Увеличиваем n, если n≠Nfft запускаем цикл заново.
}}
Сохраняем в ОЗУ a(k) и b(k), обнуляем аккумуляторы.
Увеличиваем k, если k≠Nfft/2(-1) запускаем цикл заново.
}

Очень удобно, что у AVR три регистровые пары X, Y и Z для косвенной адресации — одна будет адресовать отсчёты синуса, другая — входного сигнала, третья — коэффициенты a(k) и b(k). Если вторая и третья пары будут просто инкрементироваться, то с первой будет чуть хитрее.
Для получения a(1) и b(1), s[n] перемножается на каждый отсчёт синуса, тут всё просто.
Для получения a(2) и b(2), s[n] перемножается на каждый второй отсчёт синуса в диапазоне от 0 до 4*pi*(Nfft-1)/Nfft, в то время как таблица у нас только до 2*pi. Здесь мы воспользуемся периодичностью синуса и кратностью степени двойки количества отсчётов: 4*pi*(Nfft-1)/Nfft = 2*pi*(Nfft-2)/Nfft. Говоря языком ассемблера, будем увеличивать адрес ОЗУ отсчётов синуса не на единицу, а на k. А чтобы адрес таблицы не вылез на пределы самой таблицы, будем накладывать маску вида 0m (напр., 00111111 для Nfft=64) после каждого увеличения.

Пока что всё. На следующий уровень спустимся попозже.

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

Больше
12 года 5 мес. назад #3 от ARV
ARV ответил в теме Re: Преобразование Фурье для AVR
да вам надо статью писать :)

я не ленивый, я энергосберегающий...

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

Больше
12 года 5 мес. назад #4 от Gudd-Head
Gudd-Head ответил в теме Re: Преобразование Фурье для AVR

ARV пишет: да вам надо статью писать :)

Да я думал... времени нет.
А так очень удобно: что-то дома, что-то на работе дописать.
Если в итоге получится что-то дельное, напишу (большую часть можно будет скопировать прямо отюда).
Кстати, я прогонял этот алгоритм (с округлениями) в МатЛабе и сравнивал с обычным БПФ — отличия мизерные :)

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

Больше
12 года 5 мес. назад #5 от ARV
ARV ответил в теме Re: Преобразование Фурье для AVR
так алгоритм - классика в общем-то...
а применение какое видится?

я не ленивый, я энергосберегающий...

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

Больше
12 года 5 мес. назад #6 от Gudd-Head
Gudd-Head ответил в теме Re: Преобразование Фурье для AVR

ARV пишет: а применение какое видится?

Эээ... Возможно, нарисовать спектр сигнала на экране осциллографа.

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

Больше
12 года 5 мес. назад #7 от Gudd-Head
Gudd-Head ответил в теме Re: Преобразование Фурье для AVR
Набросал программку на АСМе. Получилось около 100 байт + коэффициенты синуса (ещё Nfft байт). Будет время - поотлаживаю и протестирую.

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

Больше
12 года 5 мес. назад #8 от ARV
ARV ответил в теме Re: Преобразование Фурье для AVR
как-то уж очень подозрительно мало... 50 команд всего лишь, что ли?!

я не ленивый, я энергосберегающий...

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

Больше
12 года 5 мес. назад #9 от Gudd-Head
Gudd-Head ответил в теме Re: Преобразование Фурье для AVR

ARV пишет: как-то уж очень подозрительно мало... 50 команд всего лишь, что ли?!

Ну да :blush: Два вложенных цикла.
Подразумевается, что всё уже в ОЗУ.
Цикл занесения массива коэфициентов синуса съест ещё пару байт + сами коэффициенты, как я уже говорил.
На выходе получаем амплитуды синусов и косинусов. Чтобы найти сумму их квадратов, надо будет замутить ещё один цикл, это съест ещё памяти :whistle:

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

Больше
12 года 5 мес. назад - 12 года 5 мес. назад #10 от Gudd-Head
Gudd-Head ответил в теме Re: Преобразование Фурье для AVR
Работает! :woohoo:
Пока что проверил только на первой гармонике синуса с Nfft = 16. :blush:
Правда, выполняется... более 4000 тактов :(
Зато занимает немногим более 100 байт B)
А всё-таки ещё одного счётчика косвенной адресации для ускорения не хватает: чтобы один был на отсчёты синуса, другой - на отсчёты косинуса.
Но, чувствую, в проге есть пара мест, где можно немного ускорить.
Последнее редактирование: 12 года 5 мес. назад от Gudd-Head.

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

Больше
12 года 5 мес. назад #11 от ARV
ARV ответил в теме Re: Преобразование Фурье для AVR
а сколько счетчиков косвенной адресации вам надо? их же три: X, Y и Z... не хватает?

я не ленивый, я энергосберегающий...

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

Больше
12 года 5 мес. назад - 12 года 5 мес. назад #12 от Gudd-Head
Gudd-Head ответил в теме Re: Преобразование Фурье для AVR

ARV пишет: а сколько счетчиков косвенной адресации вам надо? их же три: X, Y и Z... не хватает?

Я знаю что их три — я сам об этом писал =)
А ещё я написал что трёх счётчиков хватает - просто с четырьмя было бы чуток быстрее.
Кстати, рожицы blush oops одинаковые (по кр. мере, у меня): :blush: :oops: UPD: вместо УПС вставляется смайлик БЛУШ.
Теоретически оперативки хватит на 256-точечное преобразование... на которое уйдёт примерно 16 млн. тактов :S
Последнее редактирование: 12 года 5 мес. назад от Gudd-Head.

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

Больше
12 года 5 мес. назад #13 от ARV
ARV ответил в теме Re: Преобразование Фурье для AVR
как-то вы неправильно считаете... БПФ тем и славится, что при удвоении числа семплов не происходит удвоения вычислительных затрат... точную формулу не помню, но вы ее наверняка знаете. и из практики: 128-точечное 16-битное (!!!) БПФ в библиотеке ЧЕНА выполняется примерно 17 миллисекунд при тактовой 8 мГц, так что ваш вариант на 256 точек ну никак не может выполняться 1 секунду при тактовой 16 мГц!

я не ленивый, я энергосберегающий...

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

Больше
12 года 5 мес. назад - 12 года 5 мес. назад #14 от Gudd-Head
Gudd-Head ответил в теме Re: Преобразование Фурье для AVR

ARV пишет: как-то вы неправильно считаете... БПФ тем и славится, что при удвоении числа семплов не происходит удвоения вычислительных затрат... точную формулу не помню, но вы ее наверняка знаете. и из практики: 128-точечное 16-битное (!!!) БПФ в библиотеке ЧЕНА выполняется примерно 17 миллисекунд при тактовой 8 мГц, так что ваш вариант на 256 точек ну никак не может выполняться 1 секунду при тактовой 16 мГц!

А я и не говорил что у меня БПФ :silly: У меня простое ПФ... зато компактное. Да, насчёт 16 млн операций для 256-точечного я загнул — их 1,1 млн... и 280 тыщ для 128-точечного (вычислительная сложность Nfft²), т.е. 35 мс на 8 МГц при 100 байт кода. ЙА гений??? :woohoo:
Последнее редактирование: 12 года 5 мес. назад от Gudd-Head.

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

Больше
12 года 5 мес. назад #15 от ARV
ARV ответил в теме Re: Преобразование Фурье для AVR
на счет гениальности: поживем - увидим :)

а почему не БПФ? будет ведь намного быстрее!

я не ленивый, я энергосберегающий...

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

Работает на Kunena форум