FFT на Си для AVR

Обсуждаем контроллеры компании Atmel.
Ответить
ARV
Ум, честь и совесть. И скромность.
Аватара пользователя
Сообщения: 18682
Зарегистрирован: Чт дек 28, 2006 08:19:56
Откуда: Новочеркасск

Сообщение ARV »

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

Мой уютный бложик... заходите!
Контактная информация:
Реклама
YS
Друг Кота
Аватара пользователя
Сообщения: 7518
Зарегистрирован: Вс мар 29, 2009 22:09:05

Сообщение YS »

Персональное спасибо YS за спокойный и вразумительный ликбез!
Не за что, не за что. :) Спокойствие - одна из конфуцианских добродетелей. :)

По пути я и сам во многом разобрался, ибо тоже пришлось подумать. :)))
Разница между теорией и практикой на практике гораздо больше, чем в теории.
Контактная информация:
Реклама
Встал на лапы
Сообщения: 116
Зарегистрирован: Ср апр 28, 2010 17:49:14

Сообщение zenek »

Здравствуйте!
Есть сигнал который подвергается дискретному преобразованию фурье.
На выходе ДПФ есть мнимая и действительная части, корень из суммы квадратов которых есть амплитуда определенной частоты в сигнале.

Вот написал небольшую программку, которая генерирует сигнал и проводит над ним прямое преобразование фурье:

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

#include "stdafx.h"

#include <stdio.h>
#include <math.h>
#include <string.h>
#define PI 3.141592653589793238462643383279


#define SAMPLERATE 8000 
#define signal1_HZ 125  
#define signal2_HZ 375
#define DFT_SIZE 64     // êîëè÷åñòâî âûáîðîê äëÿ ÄÏÔ


signed char buff[SAMPLERATE];
float fbuf[SAMPLERATE];


void Generate (void)
{
    FILE* myfile;
    FILE* file; 
    int i;
    int ch;
    myfile = fopen("C:\\file.wav", "w"); 

    for(  i = 0; i < SAMPLERATE ;i++)
	{
		
		fbuf[i]=64*((sin(2*PI*i*signal1_HZ/SAMPLERATE))+(sin(0.8+2*PI*i*signal2_HZ/SAMPLERATE))); // äâà ñèíóñà
		//fbuf[i]=42*((sin(2*PI*i*signal1_HZ/SAMPLERATE))+(sin(0.8+2*PI*i*signal2_HZ/SAMPLERATE))+(sin(2*PI*i*1500/SAMPLERATE))); // òðè ñèíóñà
		fputc((signed char)fbuf[i],myfile);
		//printf("%f\n",fbuf[i]); 
    } 

    fclose(myfile);
 
}


float Re[DFT_SIZE],Im[DFT_SIZE],Ampl[DFT_SIZE],summ=0;


int main(void)
{
	int i=0,j=0;

	Generate();

	for(i=0;i<DFT_SIZE;i++){

		for(j=0;j<DFT_SIZE;j++){
			Re[i]+=(float)fbuf[j]*cos(2*PI*i*j/DFT_SIZE);
			Im[i]-=(float)fbuf[j]*sin(2*PI*i*j/DFT_SIZE);
		}

		Ampl[i] = (float)sqrt(Re[i]*Re[i] + Im[i]*Im[i]);
	}


	for(i=0;i<DFT_SIZE;i++){
		printf("Amp[%d]\t\t(%d Hz)\t = %f\n",i,(SAMPLERATE/DFT_SIZE)*i,Ampl[i]);
	}

	while(1);
	return 0;

}
Вот результат работы программы для сигнала - синус 125Гц + синус 375Гц:
Изображение

А вот для просто синуса 375Гц:
Изображение

Для понимания, вот сами сигналы:
Изображение

Частота дискретизации сигнала - 8000Гц.
Теперь вопрос!
Как определить реальную амплитуду? В результатах работы программы амплитуда пляшет для двух синусов - по 2048 на частоту, для одного синуса - 4064 единиц. Что это за единицы? Как из привести к нормальному виду? И как из них получить в дб, какое значение принимать за 0Дб?
Друг Кота
Аватара пользователя
Сообщения: 4905
Зарегистрирован: Чт апр 11, 2013 11:19:59
Откуда: Минск

Сообщение WiseLord »

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

for(i=0; i<DFT_SIZE; i++) {
  for(j=0; j<DFT_SIZE; j++){
Цикл 64 раза по 64 раза. Вот и выходит, если частота одна - на ней и выходит максимум (4096, если без ошибок округления), если комбинация частот - суммарная амплитуда тоже будет такой.

Так что надо результат (и действительную часть, и мнимую) на DFT_SIZE разделить, чтобы получить амплитуды в тех же "попугаях", в которых были исходные сигналы.
Контактная информация:
Реклама
Эиком - электронные компоненты и радиодетали
YS
Друг Кота
Аватара пользователя
Сообщения: 7518
Зарегистрирован: Вс мар 29, 2009 22:09:05

Сообщение YS »

zenek, ваш листинг - наглядная демонстрация того, почему я всегда всем советую писать комментарии только на английском.

Вывод консольной программы нет смысла прикладывать картинкой. Его лучше просто скопировать из консоли.

По теме - в реальности для установления соответствия результата FFT и реального входного сигнала лучше провести калибровку.
Разница между теорией и практикой на практике гораздо больше, чем в теории.
Контактная информация:
Реклама
INA
Первый раз сказал Мяу!
Аватара пользователя
Сообщения: 27
Зарегистрирован: Пн мар 23, 2015 17:18:47

Сообщение INA »

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

Продолжив тему про FFT, написал небольшую софтинку под Android... :)

Вот то, что получилось: https://youtu.be/eDiQ7_FXkuk

Это немножко не в эту тему форума... надеюсь, меня простят... :)

P.S. А это: https://youtu.be/WyIumlJ8TME то, что было сделано раньше на ATMega32... :)
Реклама
INA
Первый раз сказал Мяу!
Аватара пользователя
Сообщения: 27
Зарегистрирован: Пн мар 23, 2015 17:18:47

Сообщение INA »

И ещё немного... :)
https://youtu.be/L4_-yqQE4t8
Держит паяльник хвостом
Сообщения: 963
Зарегистрирован: Чт апр 03, 2014 23:16:55
Откуда: Россия

Сообщение Rtmip »

INA писал(а):Добрый вечер всем!

Продолжив тему про FFT, написал небольшую софтинку под Android... :)

Вот то, что получилось: https://youtu.be/eDiQ7_FXkuk

Это немножко не в эту тему форума... надеюсь, меня простят... :)
Игорь Николаевич, я простил, да! :music: Но где же она? :dont_know:
"Софтинка" где? Интересно попробовать.
INA
Первый раз сказал Мяу!
Аватара пользователя
Сообщения: 27
Зарегистрирован: Пн мар 23, 2015 17:18:47

Сообщение INA »

Доброй ночи!

Да вот она...
Правда, думаю что не на всех планшетах будет нормально выглядеть интерфейс...
Вложения
Spectre.rar
(483.74 КБ) 233 скачивания
Родился
Сообщения: 1
Зарегистрирован: Пн июл 21, 2014 21:08:02

Сообщение Ypir »

Добрый вечер! Попробовал реализовать FFT, но возникло несколько вопросов.
Вопрос 1: При получении спектра результат напоминает пилу, т.е. переменная am нарастает до 0xFF затем обнуляеться и снова растет, должно ли так быть?

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

{
      for(i1=0; i1<32; i1++)
      {
      A[i1]=read_adc(0x20);
      
      }
      
      re = im = 0x0000;
	for(i2=0;i2<32;i2++)
	{
	   for(i4=0;i4<16;i4++)
	   {
	     smim = A[i2]*sin16[i4];
	     smre = A[i2]*cos16[i4];
	     im = im + smim;	   
	     re = re + smre;
	   }
	}
	
         tempsqr=(re*re)+(im*im);
	am=iqrt(tempsqr); 
	
Вопрос 2: Если результат преобразования (переменная am) засунуть в шим, банально OCR0A=am, получается что то непонятное, собственно подскажите пожалуйста как обработать результат что бы засунуть его в шим и скажем помигать диодом? Как я понимаю чем частота ближе к искомой тем чаще он должен моргать, правда это при небольшой отстройку по частоте.
Грызет канифоль
Сообщения: 274
Зарегистрирован: Чт апр 16, 2009 14:23:59

Сообщение NebelWefer »

Товарищи, правильно ли я понял: если мне при анализе звукового сигнала не интересны частоты выше 15кГц, то я могу снизить частоту дискретизации до 30кГц?
Друг Кота
Аватара пользователя
Сообщения: 4905
Зарегистрирован: Чт апр 11, 2013 11:19:59
Откуда: Минск

Сообщение WiseLord »

Да. Но перед этим нужно отфильтровать частоты выше 15, чтобы их на входе и не появлялось. Иначе 16кГц будет отображаться в 14кГц, 20кГц будет видно как 10кГц и т.д.
Контактная информация:
Грызет канифоль
Сообщения: 274
Зарегистрирован: Чт апр 16, 2009 14:23:59

Сообщение NebelWefer »

Ааа, ясно, иными словами, частота дискретизации ВСЕГДА должна вычисляться из максимальной частоты сигнала, независимо от потребностей. Спасибо за пояснение!
После прочтения темы возникла мысль, что не обязательно получать весь массив амплитуд частот на выходе.
То есть:
Получили массивы RE и IM, затем только берем только те k , которые соответствую нужным частотам, исходя из примера товарища Jordan
jordan писал(а):

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

for (int k = 0; k < 1024; k++)
            {
                for (int n = 0; n < 1024; n++)
                {
                    x_real[k] += x_n[n] * Math.Cos( (-2) * Math.PI * k * n / 1024 );
                    x_imag[k] += x_n[n] * Math.Sin( (-2) * Math.PI * k * n / 1024);
                }

                x_Ampl[k] = Math.Sqrt( x_real[k] * x_real[k] + x_imag[k] * x_imag[k] );

                
            }
Если мне нужны амплитуды, например 1кГц и 5кГц, то я выбрасываю верхний цикл и отдельно считаю для k=(fгарм * SIZE) / fдискр.
к1=1000*1024/40000= 25 или 26,
к2=5000*1024/40000= 128
Верно ли я вывел?
Друг Кота
Аватара пользователя
Сообщения: 25437
Зарегистрирован: Чт янв 10, 2008 22:01:02
Откуда: Московская область, Фрязино

Сообщение КРАМ »

Для получения нескольких гармоник не нужно делать БПФ. ДПФ в таком случае экономичнее.
Грызет канифоль
Сообщения: 274
Зарегистрирован: Чт апр 16, 2009 14:23:59

Сообщение NebelWefer »

NebelWefer писал(а): Если мне нужны амплитуды, например 1кГц и 5кГц, то я выбрасываю верхний цикл и отдельно считаю для k=(fгарм * SIZE) / fдискр.
к1=1000*1024/40000= 25 или 26,
к2=5000*1024/40000= 128
Верно ли я вывел?
Получается верно)
Реализовал этот код "в лоб", работает, частоты выделяются. Теперь дело только за оптимизацией.
Но хотелось бы еще уточнить по выбору размера массива: понятно, что чем больше массив, тем выше разрешающая способность по частотам, а вот на точность результата влияет ли?
Например, есть массивы: 625 Гц, 1250, 1875... и 312 625 937...
Одинаковый результат будет для 625 Гц? Будет ли во 2м случае влияние соседних частот меньше?
Друг Кота
Аватара пользователя
Сообщения: 25437
Зарегистрирован: Чт янв 10, 2008 22:01:02
Откуда: Московская область, Фрязино

Сообщение КРАМ »

"Разрешающая способность по частотам" вообще то называется полосой анализа. И эта полоса обратно пропорциональна ВРЕМЕНИ НАКОПЛЕНИЯ. АЧХ фильтра приведенная к полосе анализа всегда одинакова. Поэтому избирательность по соседнему фильтру не зависит от числа фильтров. Причем если бы было иначе, АЧХ анализатора представляла собой волнистую линию.
Энергия сигнала попадающего ровно между фильтрами поровну распределится между этими фильтрами.
Грызет канифоль
Сообщения: 274
Зарегистрирован: Чт апр 16, 2009 14:23:59

Сообщение NebelWefer »

КРАМ писал(а):Энергия сигнала попадающего ровно между фильтрами поровну распределится между этими фильтрами.
Получается при большем массиве результат точнее?
625 1250 1875
312 625 937 1250,
если в сигнале есть составляющая 900 Гц, то в первом случае ее присутствие отразится и на 625 и на 1250, а во 2м случае на 625 и 937 не попав в 1250.
Верно?
Друг Кота
Аватара пользователя
Сообщения: 25437
Зарегистрирован: Чт янв 10, 2008 22:01:02
Откуда: Московская область, Фрязино

Сообщение КРАМ »

Естественно, чем уже полоса анализа, тем точнее можно определить спектральный максимум сигнала. Но при этом надо понимать, что реальные сигналы имеют достаточно размытый спектр, хотя бы в силу их ограниченности во времени. То есть у короткого во времени сигнала не может быть узкого спектра. А строго монохроматический спектр имеет бесконечная во времени синусоида. Сиречь, НИКАКАЯ реальная синусоида не может быть строго монохроматичной.
Грызет канифоль
Сообщения: 274
Зарегистрирован: Чт апр 16, 2009 14:23:59

Сообщение NebelWefer »

Предлагаю еще немного развить тему для полноты картины. Везде упоминаются таблицы sin/cos для ускорения вычислений. Я попробовал так: рассчитал в экселе для K=1, в проект добавил массив, умноженный на 10000:

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

const short sin_k1[N] = {10000, 9998 .......}
теперь кусок:

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

x_real[k] += x_n[n] * Math.Cos( (-2) * Math.PI * k * n / 1024 );
стал просто:

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

x_real[k] += x_n[n] *sin_k1[n]/10000
работает действительно быстрее, но на такой массив уходит куча памяти, а еще и для других К нужны.
Как вы обходите этот момент?
Друг Кота
Аватара пользователя
Сообщения: 25437
Зарегистрирован: Чт янв 10, 2008 22:01:02
Откуда: Московская область, Фрязино

Сообщение КРАМ »

Любая таблица, а особенно большая, должна размещаться в программном фдеше. А МК заточенные, в том числе, под применение таблиц должны иметь инструмент быстрой выборки данных из флеша, для чего область флеша так или иначе отражена в адресное пространство ОЗУ.
В Си для МК имеются различные виды синтаксиса и квалификаторов данных,чтобы компилятор размещал константы в программном флеше.
Ответить

Вернуться в «AVR»