FFT на Си для AVR

[ARV]

можно залить на yputube и разместить здесь...

[YS]

Персональное спасибо YS за спокойный и вразумительный ликбез!

Не за что, не за что. Спокойствие - одна из конфуцианских добродетелей.

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

[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Дб?

[WiseLord]

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

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

Цикл 64 раза по 64 раза. Вот и выходит, если частота одна - на ней и выходит максимум (4096, если без ошибок округления), если комбинация частот - суммарная амплитуда тоже будет такой.

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

[YS]

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

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

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

[INA]

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

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

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

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

P.S. А это: https://youtu.be/WyIumlJ8TME то, что было сделано раньше на ATMega32...

[INA]

И ещё немного...
https://youtu.be/L4_-yqQE4t8

[Rtmip]

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

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

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

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

Игорь Николаевич, я простил, да! Но где же она?
"Софтинка" где? Интересно попробовать.

[INA]

Доброй ночи!

Да вот она...
Правда, думаю что не на всех планшетах будет нормально выглядеть интерфейс...

[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, получается что то непонятное, собственно подскажите пожалуйста как обработать результат что бы засунуть его в шим и скажем помигать диодом? Как я понимаю чем частота ближе к искомой тем чаще он должен моргать, правда это при небольшой отстройку по частоте.

[NebelWefer]

Товарищи, правильно ли я понял: если мне при анализе звукового сигнала не интересны частоты выше 15кГц, то я могу снизить частоту дискретизации до 30кГц?

[WiseLord]

Да. Но перед этим нужно отфильтровать частоты выше 15, чтобы их на входе и не появлялось. Иначе 16кГц будет отображаться в 14кГц, 20кГц будет видно как 10кГц и т.д.

[NebelWefer]

Ааа, ясно, иными словами, частота дискретизации ВСЕГДА должна вычисляться из максимальной частоты сигнала, независимо от потребностей. Спасибо за пояснение!
После прочтения темы возникла мысль, что не обязательно получать весь массив амплитуд частот на выходе.
То есть:
Получили массивы RE и IM, затем только берем только те k , которые соответствую нужным частотам, исходя из примера товарища 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
Верно ли я вывел?

[КРАМ]

Для получения нескольких гармоник не нужно делать БПФ. ДПФ в таком случае экономичнее.

[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м случае влияние соседних частот меньше?

[КРАМ]

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

[NebelWefer]

Энергия сигнала попадающего ровно между фильтрами поровну распределится между этими фильтрами.

Получается при большем массиве результат точнее?
625 1250 1875
312 625 937 1250,
если в сигнале есть составляющая 900 Гц, то в первом случае ее присутствие отразится и на 625 и на 1250, а во 2м случае на 625 и 937 не попав в 1250.
Верно?

[КРАМ]

Естественно, чем уже полоса анализа, тем точнее можно определить спектральный максимум сигнала. Но при этом надо понимать, что реальные сигналы имеют достаточно размытый спектр, хотя бы в силу их ограниченности во времени. То есть у короткого во времени сигнала не может быть узкого спектра. А строго монохроматический спектр имеет бесконечная во времени синусоида. Сиречь, НИКАКАЯ реальная синусоида не может быть строго монохроматичной.

[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

работает действительно быстрее, но на такой массив уходит куча памяти, а еще и для других К нужны.
Как вы обходите этот момент?

[КРАМ]

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