FFT на Си для AVR

[YS]

Вроде как всё правильно сделал, но в Протеусе не хочет работать...

Я правильно понимаю, что в прошивке вы, как и следует, не использовали вещественных чисел и тригонометрических функций, а сделали на фиксированной точке и с предрассчитанными массивами? Вещественные вычисления на AVR работают ОЧЕНЬ медленно по причине отсутствия математического сопроцессора. Тот пример на Lua, что я приводил выше, хорош в основном для демонстрационных целей. Для реальной работы его надо адаптировать.

[INA]

Добрый день!
Рад, что тема движется...
Вот как раз сейчас хотел спросить у Вас каким образом вот здесь
были посчитаны коэффициенты для частоты 150Гц?
Опять заблудился...


А сам код такой:

Код: Выделить всё

//Программа инициализации АЦП
  void adc_init()
{
	ADMUX = (1<<REFS0); //выбираем источник питания АЦП 5v
	ADMUX |= (0<<MUX0);
	ADCSRA=(1<<ADEN)|(1<<ADPS2)|(1<<ADATE)|(1<<ADIE); //включаем АЦП с делителем = Fcpu/16
	ADCSRA |= (1<<ADSC);
 
}
// 

Ну и обработка прерывания:

Код: Выделить всё

 ISR(ADC_vect) //подпрограмма обработки прерывания от АЦП
{
 //cli();
 PORTD = 1;
 buffer[i] = (ADCL + (ADCH<<8))-512;
 i++;
 count++;
 if(count == BUF_SIZE) {cli(); i=0; count=0;}

   PORTD = 0;

}

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

[YS]

Не совсем понял, про какие коэффициенты вы говорите. Ну а код просто настраивает АЦП и заполняет массив. Обработки сигнала в приведенных фрагментах нет.

Кстати, запись

Код: Выделить всё

buffer[i] = (ADCL + (ADCH<<8))-512;

не особо хороша. Дело в том, что в Си не определен порядок вычисления операндов выражения, в то время как ADCL и ADCH требуют совершенно определенного порядка доступа (даташит, стр. 258):

When an ADC conversion is complete, the result is found in these two registers.
When ADCL is read, the ADC Data Register is not updated until ADCH is read. Consequently, if
the result is left adjusted and no more than 8-bit precision is required, it is sufficient to read
ADCH. Otherwise, ADCL must be read first, then ADCH.
The ADLAR bit in ADMUX, and the MUXn bits in ADMUX affect the way the result is read from
the registers. If ADLAR is set, the result is left adjusted. If ADLAR is cleared (default), the result
is right adjusted.

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

То есть, я бы написал как-то так:

Код: Выделить всё

int16_t tmp;

...

tmp=ADCL;
tmp|=((uint16_t)ADCH)<<8;
tmp-=512;

buffer[i]=tmp;

...

[ARV]

То есть, я бы написал как-то так:

Код: Выделить всё

int16_t tmp;

...

tmp=ADCL;
tmp|=((uint16_t)ADCH)<<8;
tmp-=512;

buffer[i]=tmp;

а почему бы не делать ПРАВИЛЬНО? вы пользуетесь WinAVR (или AVR-GCC), так почему не пользоваться для получения значения АЦП предопределенной переменной ADC? в CodeVision для аналогичных целей используется ADCW.
к чему морочить себе голову с доступом к половинкам регистра АЦП, если компилятор готов все сделать за вас?

Код: Выделить всё

int tmp = ADC - 512;

[INA]

Не совсем понял, про какие коэффициенты вы говорите.

Да вот тут он пишет:

Например имеем 512 отсчетов АЦП нужно посчитать мнимую и действительную части для 150Гц при частоте дискретизации 19200 Гц:

во вложении картинка с коэффициентами...

Не пойму по какой формуле он всё это считал?...
То, что потом всё умножается на 127 - это понятно..
А как изначально?...
Пробовал повторить этот расчёт - ничего не получается...

[YS]

во вложении картинка с коэффициентами...

Дык это просто синус и косинус.

Помните я писал

Код: Выделить всё

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;

Вот он, если переносить на мой пример, заранее посчитал N значений

Код: Выделить всё

math.cos(2*math.pi*REL_FREQ*(k/N))

и

Код: Выделить всё

math.sin(2*math.pi*REL_FREQ*(k/N))

Ну и округлил/нормировал, разумеется.

[INA]

Код: Выделить всё

math.cos(2*math.pi*REL_FREQ*(k/N))

Пробовал таким образом - ничего не получается...
Не сходится с его таблицей...

К примеру 8-е значение в таблице косинусов - 117, а у меня получается с округлением - 65...

[YS]

К примеру 8-е значение в таблице косинусов - 117, а у меня получается с округлением - 65...

Так абсолютные значения зависят от нормировки. Главное - вид функции. Постройте его массив, ваш массив, и сравните.

[INA]

К примеру 8-е значение в таблице косинусов - 117, а у меня получается с округлением - 65...

Так абсолютные значения зависят от нормировки. Главное - вид функции. Постройте его массив, ваш массив, и сравните.

Всё равно не пойму почему не сходятся результаты... )
Видимо, не мой день... )

[INA]

Ткните просто меня носом, как он вычислил эти коэффициенты?...
Хожу вокруг да около и ничего не понимаю....
Автор пишет, что вычислял коэфф. для частоты 150Гц при частоте дискретизации 19200...
Как всё это связать, чтоб получилось вычислить то же самое?

[Morroc]

Чисто интуитивно как то так - 19200 / 150 - число отсчетов, теперь укладываем в них полный период синуса/косинуса - оно ?

в excel =ОКРУГЛ(COS(1/( 19200 / 150 )*8*360*3,14/180)*127;4)

8й отсчет в последовательности из 128

[YS]

а почему бы не делать ПРАВИЛЬНО? вы пользуетесь WinAVR (или AVR-GCC), так почему не пользоваться для получения значения АЦП предопределенной переменной ADC?

В дидактических целях. Кроме AVR и AVR-GCC существуют и другие контроллеры и компиляторы, и там тоже бывают подобные требования, но не всегда определены соответствующие "прямые" возможности. Ну а в данном конкретном случае, конечно, можно и так.

Ткните просто меня носом, как он вычислил эти коэффициенты?...

Ох. Ну окей, я написал сходу и "влоб":

Код: Выделить всё

SAMPLING_FREQ_HZ	= 19200;
SINE_FREQ_HZ		 = 150;
SAMPLES_NO			= 128;

timestep = 1/SAMPLING_FREQ_HZ;
sine_table = {};
k=0;

while k<SAMPLES_NO do

	ph=2*math.pi*SINE_FREQ_HZ*k*timestep;
	sine_table[k+1]=math.floor(127*math.sin(ph));
	k=k+1;

end;

k=0;
while k<SAMPLES_NO do
	io.write(sine_table[k+1] .. ",");

	if ((k+1)%16)==0 then
		print("");
	end;

	k=k+1;
end;

Результат:

Код: Выделить всё

0,6,12,18,24,30,36,42,48,54,59,65,70,75,80,85,
89,94,98,102,105,108,112,114,117,119,121,123,124,125,126,126,
127,126,126,125,124,123,121,119,117,114,112,108,105,102,98,94,
89,85,80,75,70,65,59,54,48,42,36,30,24,18,12,6,
0,-7,-13,-19,-25,-31,-37,-43,-49,-55,-60,-66,-71,-76,-81,-86,
-90,-95,-99,-103,-106,-109,-113,-115,-118,-120,-122,-124,-125,-126,-127,-127,
-127,-127,-127,-126,-125,-124,-122,-120,-118,-115,-113,-109,-106,-103,-99,-95,
-90,-86,-81,-76,-71,-66,-60,-55,-49,-43,-37,-31,-25,-19,-13,-7

Сравним:



Как я уже говорил, это просто массив синуса, без всяких изысков. Надеюсь, разница в единицы (ошибка округления) вас не смущает?

[INA]

Доброе утро всем!

Вчера не смог уже ответить - уснул за столом...

YS, спасибо за очередную подсказку (точнее ликбез)!

Но у меня опять ничего не выходит...

Считал просто на калькуляторе:

ph = 6.28 * 150 * 8 * 0.000052 = 0.391872;
sine_table[k+1] = 127 * cos(0.391872) = 126.99; - даже близко ничего похожего...

Что я не так делаю?
Не могу понять...

Дело в вычислении косинуса на калькуляторе...
В этом моменте я делаю что-то не так!...

Пересчитал на этом калькуляторе - всё правильно!...

[INA]

нашёл причину...
Угол на калькуляторе не поставил в радианы...

[YS]

Ага, бывает.

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

[INA]

Ага, бывает.

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

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

Поэтому взял табличное значение, просто на выбор 8-е и попытался его получить так как описано у автора...
Ну и забуксовал из-за этих градусов/радианов...

На данный момент вычислил таблицы для частоты 1000Гц, подключил отдельным файлом в программу, скорректировал последнюю и теперь вижу как это работает а Протеусе...
Протеус для таких вещей конечно же не есть хорошо, но то, что он показывает на данном этапе мне достаточно...
Далее нужно спаять платку (прошлый раз сделать это на работе не получилось ) и доводить всё это до ума.

Вот... как-то так...

[INA]

Добрый день!

В общем, в Протеусе нарисовал вот такую схемку:



С генератора подаю сигнал, а на выходе ЦАП получаю аналоговый уровень...
В общем, для теста вполне достаточно...
Сейчас паяю платку...

[YS]

INA, поскольку мне больше нечего сказать по существу, кроме того, что я рад, что у вас, как я понял, все работает, можно я немного поработаю за модераторов, пока они вас сами не нашли?

Дело в том, что правила форума запрещают цитировать сообщение полностью, равно как и создавать подряд несколько сообщений одному автору. За это модераторы могут обратить свое око на вас (забанить, скорее всего, не забанят, но предупреждение вынесут). Картинок более 800 пикселей в ширину вы пока не выкладывали, но знайте, что это тоже частая ошибка - пару раз и я попадался на этом.

[INA]

Добрый вечер!

Да нет, конечно!
Как я могу быть против, если судя по надписи над аватаром,
я "Первый раз сказал Мяу!"...

[INA]

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

В общем, разобрался я со своей задачей.
Спасибо большое всем за помощь в понимании не понимаемого...
Персональное спасибо YS за спокойный и вразумительный ликбез!...

Чтобы как-то закрепить информацию и увидеть как всё это работает в реальности,
собрал на Меге32 небольшой 5-ти полосный анализатор спектра...
Собрал его на китайской дырявой платке, но тем не менее всё работает как положено...
Чему, естественно безмерно рад!...
Записал видео его работы, но файл большой и здесь выложить не получится, однако если кому-нибудь будет интересно, могу скинуть на какой-нибудь файлообменник...
Ещё раз спасибо всем!