Перейти к содержанию
Поиск
    

Полосовой фильтр с помощью CMSIS DSP

Ответить

smk    0

smk
Опубликовано · Жалоба
45 минут назад, Михась сказал:

Начните с Winfilter, действительно.

Начинал. Надо понять как должно быть. как понять работает или нет? Нужен рабочий пример на котором я все увижу. Дальше разберусь.

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

TU-104    0

TU-104
Опубликовано (изменено) · Жалоба
7 hours ago, smk said:

С изменениями разобрался. Всеравно не то. Пример бы рабочий.

Вот рабочий пример, на паскале правда, возможно поможет, синтаксис там понятный.

Делал себе эмулятор когда-то, чтоб разобраться.

В программе можно сгенерировать 2 частоты F0+F1, чтоб посмотреть, как алгоритм отрабатывает на таких сигналах.

Для начала можно задать одну F0, тогда амплитуду F1 задать =0

Кнопка "Генерировать" - создаёт массив точек сгенерированного сигнала (число точек задаётся, Fs у меня всегда 8000 Гц) и выводит график.

Дальше задаётся длина интервала оценки, на которые будет разбит сигнал и на каждом будет запущен алгоритм Герцеля, задаётся искомая частота и нажимается кнопка "Gortc proc". Синим цветом у меня вычисляется амплитуда всего сигнала за вычетом DC на интервале, и при одной частоте видно как к концу интервала амплитуда равна числу, которое выдает алгоритм герцеля (он выдаёт амплитуду на искомой частоте)

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

image.thumb.png.e146fc4fdebe43c641231361f80b0b26.png

 

В архиве ехе-шник и .pas, часть нужного кода вставляю сюда: 

Procedure Gortc_Proc();
  //const a1 = 0.618;      // 2•cos(2PI•k)   k=f/Fd
  var D0, D1, D2 : real;
   iD0, iD1, iD2 : Int64; //integer;
      i : integer;
      P : real;
      Ff: integer;
      a1: real;
      iA1: Int64;
      NumSamples : integer;  // число выборок для интервал оценки. с новым интервалом фильтр обнуляется

begin

  Fmain.Chart1.Series[3].Clear;
  Fmain.Chart1.Series[4].Clear;
  Fmain.Chart1.Series[5].Clear;

  Ff := Fmain.SpinEditFfind.Value;
  a1 :=2*cos(2*Pi*Ff/Fd);

  NumSamples := Fmain.SpinEditInterval.Value;

//  D0 := 0;
  D1 := 0;
  D2 := 0;
  for i:= 0 to NumPoints-1 do begin
    D0 := X_IN[i]+ a1*D1-D2;
    D2 := D1;
    D1 := D0;

    P := sqrt(D1*D1+D2*D2 - a1*D1*D2);
    P := P/(NumSamples/2);

//    pp := round(P) div 10;
    Fmain.Chart1.Series[4].AddXY(i, P);
      // сброс коэффициентов через N мс
    if (i mod NumSamples) =0   // 80*125us = 10ms
    then begin
      mLog(IntToStr(i)+'. AMP(f) на интервале '+IntToStr(round(P)));
      D1:=0;
      D2:=0;
    end;

  end;

    

   //вариант с целыми
  iA1 := round(a1*256);
  iD1 := 0;
  iD2 := 0;
  for i:= 0 to NumPoints-1 do begin
    iD0 := X_IN[i]+ ((iA1*iD1) div 256) -iD2;
    iD2 := iD1;
    iD1 := iD0;

    P := sqrt(abs(iD1*iD1+iD2*iD2- iA1*iD1*iD2/256));
    P := P/(NumSamples/2);

//    pp := round(P) div 10;
    Fmain.Chart1.Series[5].AddXY(i, P);
      // сброс коэффициентов через N мс с новым интервалом
    if (i mod NumSamples) =0   // 80*125us = 10ms
    then begin
      mLog(IntToStr(i)+'. AMP(f) на интервале '+IntToStr(round(P)));
      iD1:=0;
      iD2:=0;
    end;

  end;



end;

 

PS. Наверно можно и в экселе сделать всё то же самое, скопировав формулы в сотню строк

 

dsp_detect_d10.zip

Изменено пользователем TU-104

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

TU-104    0

TU-104
Опубликовано (изменено) · Жалоба

PPS: Вот Си вариант, работает под cortex-M7 с float, под целые так и не переделал, не понадобилось 

// буфер 240(30мс) s16 элементов накапливаю для обработки Герцелем//  3 цикла ДМА
#define RX_DSP_LEN        240
// #define a1 :=2*cos(2*Pi*Ff/Fd);
#define a1_1200         (float)1.1755705046   //1.1755705045849462583374119092781
#define a1_1600         (float)0.6180339888   //0.61803398874989484820458683436564


// вызывается из codec_work() при накоплении буфера 30мс
void DSP_gertzel(){
  
  s32 Summ = 0;
  s16 DCoff;
  u16 Amp;
  u32 i;
  
  float D0, D1, D2;
  float Amp12, Amp16;
  
  for (i=0;i<RX_DSP_LEN; i++){
    Summ += Buf_DSP[i];
  }
  DCoff = Summ/RX_DSP_LEN; // среднее арифм. DC offset
  
  Summ = 0;
  for (i=0;i<RX_DSP_LEN; i++){
    Summ += abs(Buf_DSP[i]-DCoff);
  }  
  Amp = (u16) round(((float)(Summ/RX_DSP_LEN)*PI/2));

  
   // если общий уровень не в воротах, то просто считаю детекцию нулевой и выполняю что нужно
  if ((Amp<MIN_LVL_ADASE)||(Amp>MAX_LVL_ADASE)){
    Amp12 = 0;
    Amp16 = 0;
    //d0printD(DEBUG_DSP, "\rDC=%d Amp=%d       no proc xxx             xxx ", DCoff, Amp);
  }
  else {
    // *******   сам алгоритм   *********
    // **********   1200   **************
    D1 = 0;
    D2 = 0;
    for (i=0; i<RX_DSP_LEN; i++){
      D0 = Buf_DSP[i]+ a1_1200*D1-D2;
      D2 = D1;
      D1 = D0;
    }
    Amp12 = sqrt(D1*D1+D2*D2 - a1_1200*D1*D2);
    Amp12 = Amp12/(RX_DSP_LEN/2);
    // **********   1600   ************
    D1 = 0;
    D2 = 0;
    for (i=0; i<RX_DSP_LEN; i++){
      D0 = Buf_DSP[i]+ a1_1600*D1-D2;
      D2 = D1;
      D1 = D0;
    }
    Amp16 = sqrt(D1*D1+D2*D2 - a1_1600*D1*D2);
    Amp16 = Amp16/(RX_DSP_LEN/2);
    // *******************************     
    
  } // Amp в норме
  

  // дальше определяем, какая частота присутствует. общий уровень отфильтрован выше
  if((Amp12 > MIN_LVL_ADASE)&&(Amp16> MIN_LVL_ADASE)){
   (...)
    считаем на интервале 30мс обнаруженным
  }
    
    
}

 

Изменено пользователем TU-104

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

spirit_1    5

spirit_1
Опубликовано (изменено) · Жалоба
9 hours ago, smk said:

Начинал. Надо понять как должно быть. как понять работает или нет? Нужен рабочий пример на котором я все увижу. Дальше разберусь.

рабочий пример ниже запускаете в любом компиляторе Си

Пример свипирует диапазон 900-1100 гц, диапазон АЦП 0-2000 (введено смещение 1000)

результат работы ниже , четко виден максимум на 1000 гц

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

найдите оптимальную для вас .  Затем запустите как есть пример на контроллере, он должен отработать в режиме симуляции,   после этого уже можете подсунуть данные из ДМА

image.thumb.png.202f1ec4c2b37289ad34c0233bad1ff1.png


 

Spoiler 
/****************************************************************************
 *                                                                          *
 * Filename: hello.c                                                        *
 *                                                                          *
 * Purpose : Standard C sample for Pelles C for Windows.                    *
 *                                                                          *
 * History : Date      Reason                                               *
 *           02-09-10  Created                                              *
 *                                                                          *
 ****************************************************************************/


#include <stdio.h>
#include <stdlib.h>
#include <time.h>
 #include <math.h>

#define PI 3.14159265358979323846 
#define buffer_size 1000

#define DC_offset 200

int main()
{
int samples[buffer_size]; 
float FS = 4000.0; //fréquence d'échantillonnage 
float FDETECT = 1000.0; // частота фильтрации

int K; 
float coefficient; 
float W; 
float sine; 
float cosine; 
float Q0, Q1, Q2; 
float real; 
float imag; 
float magnitude; 
float scalingFactor; 
int i; 

K = (int) (0.5 + ((buffer_size * FDETECT) / FS));
W = (2.0 * PI * K) / buffer_size; 
cosine = cos(W); 
sine = sin(W); 
coefficient = 2 * cos(W); 
scalingFactor = buffer_size / 2.0; 


for (float Fin = 900.0; Fin < 1100.0; Fin= Fin+5.0){// перестройка частоты с шагом 10 гц в диапазоне 900- 1100гц

//---- заполняем массив тестовыми значениями синусоиды в диапазоне +_1000
//float Fin = 1000.0; // входная частота  
float Kor = (FS/Fin)/(2*PI);
for(int i = 0; i < buffer_size; i++){
samples  = (sin(i/Kor))* 1000 + DC_offset;
//printf("sample = %d  Amplituda= %d \n", i, samples );    //
    }


Q0 = 0; 
Q1 = 0; 
Q2 = 0; 
    
for (i=0 ; i<buffer_size ; i++) 
 {
  Q0 = samples + (coefficient * Q1) - Q2;
  Q2 = Q1;
  Q1 = Q0;
        //printf("test = %f \n",  (float)Q1 );
 }
real = (Q0 - (Q1 * cosine)) / scalingFactor; 
imag = (- Q1 * sine) / scalingFactor;    
magnitude = sqrt(real * real + imag * imag);

printf(" Freq_Hz= %f     Amplituda = %f \n", Fin , (float)magnitude );
}
}

 

Изменено пользователем haker_fox
Для оформления кода есть кнопка <>. Длинный код прячется под спойлер.

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

smk    0

smk
Опубликовано · Жалоба

Пример работает. Вопрос как заплоненный ДМА буфер (16 бит) правильно  подставить. 

float Kor = (FS/Fin)/(2*Pi);
for(int i = 0; i < buffer_size; i++){
samples[i]  = (sin(i/Kor))* 1000 + DC_offset;

Вот сюда.

DC_offset у меня пол питания.

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

spirit_1    5

spirit_1
Опубликовано · Жалоба
On 1/21/2023 at 11:54 AM, smk said:

Пример работает. Вопрос как заплоненный ДМА буфер (16 бит) правильно  подставить. 

float Kor = (FS/Fin)/(2*Pi);
for(int i = 0; i < buffer_size; i++){
samples[i]  = (sin(i/Kor))* 1000 + DC_offset;

Вот сюда.

DC_offset у меня пол питания.

 

2.То что вы показали , в цикле это протосто заполнение буфера тестовым сигналом вместо ДМА,поскольку у вас теперь это сделает ДМА, может удалить или закоментировать этот цил for (float Fin.....){

}

3. Весь инит  поместите перед main 

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
 #include <math.h>

#define PI 3.14159265358979323846 
#define buffer_size 1000

uint16_t samples[buffer_size]; 
float FS = 4000.0; //fréquence d'échantillonnage 
float FDETECT = 1000.0; // частота фильтрации

int K; 
float coefficient; 
float W; 
float sine; 
float cosine; 
float Q0, Q1, Q2; 
float real; 
float imag; 
float magnitude; 
float scalingFactor; 

 

 

 

в прерывании по заполнеии буфкпа ДМА 

void DMA1_Channel1_IRQHandler (void){
	
	DMA1_Channel1->CCR &= ~DMA_CCR1_EN;
	DMA1->IFCR |= DMA_IFCR_CGIF1;


Q0 = 0; 
Q1 = 0; 
Q2 = 0; 

	for ( uint32_t i=0; i<buffer_size; i++) 
	 {
		Q0 = (float)samples1[i] + (coefficient * Q1) - Q2;
		Q2 = Q1;
		Q1 = Q0;
	 }


real = (Q0 - (Q1 * cosine)) / scalingFactor; 
imag = (- Q1 * sine) / scalingFactor;    
magnitude = sqrt(real * real + imag * imag);

printf(" Freq_Hz= %f     Amplituda = %f \n", Fin , (float)magnitude );


	f = 1;
	DMA1_Channel1->CCR |= DMA_CCR1_EN;	 
}

У вас ДМА  в normal или циркуляр?

Если в циркуляр то можете не успеть обработать буфер до следующего вызова. Лучше сделать нормал , и  по необходимости запусать ДМА

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

smk    0

smk
Опубликовано · Жалоба
11 часов назад, spirit_1 сказал:

У вас ДМА  в normal или циркуляр?

Если в циркуляр то можете не успеть обработать буфер до следующего вызова. Лучше сделать нормал , и  по необходимости запусать ДМА

Да, я так и сделал. Беру выборку, обрабатываю. Заткм следующую беру.

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

jcxz    184

jcxz
Опубликовано · Жалоба
14 часов назад, spirit_1 сказал:

в прерывании по заполнеии буфкпа ДМА

printf() в ISR??? Серьёзно???  :shok:

14 часов назад, spirit_1 сказал:

magnitude = sqrt(real * real + imag * imag);

sqrt() - какой аргумент принимает? float, double?

14 часов назад, spirit_1 сказал:

DMA1_Channel1->CCR &= ~DMA_CCR1_EN;

Принудительное выключение работающего DMA-канала - крайне сомнительное "решение".

Не говоря уже о том, что даже оно у вас сделано криво.

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

jcxz    184

jcxz
Опубликовано · Жалоба
14 часов назад, spirit_1 сказал:

DMA1->IFCR |= DMA_IFCR_CGIF1;

Чтение из write-only регистра - плохая идея. И признак быдлокода. Что вы надеетесь оттуда прочитать?

14 часов назад, spirit_1 сказал:

real = (Q0 - (Q1 * cosine)) / scalingFactor; 

imag = (- Q1 * sine) / scalingFactor;

Операция деления на ARM - вещь ресурсозатратная (тактов CPU). Глупо её применять там, где она совершенно не нужна. Тем более - в ISR.

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

spirit_1    5

spirit_1
Опубликовано · Жалоба
On 1/22/2023 at 8:00 AM, smk said:

Да, я так и сделал. Беру выборку, обрабатываю. Заткм следующую беру.

Ну так должно работать.

То что написал  jcxz выше верно но вас в данный момент волновать особо не должно

Если ДМА у вас работает и вы уверенны что в массиве по прерыванию принятый верный сигнал, частота выборки ацп правильная то все должно работать 

а приведением в порядок кода зайиетесь потом

Наверное проще было бы гомнокод индийский HAL использовать и в коллбэке это делать 

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

Присоединяйтесь к обсуждению

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

Гость
Ответить в этой теме...

×   Вставлено с форматированием.   Вставить как обычный текст

  Разрешено использовать не более 75 эмодзи.

×   Ваша ссылка была автоматически встроена.   Отображать как обычную ссылку

×   Ваш предыдущий контент был восстановлен.   Очистить редактор

×   Вы не можете вставлять изображения напрямую. Загружайте или вставляйте изображения по ссылке.

×
×
×
  • Создать...