реклама на сайте
подробности

 
 
2 страниц V   1 2 >  
Reply to this topicStart new topic
> Не получается рассчитать/подобрать фильтр., LPF
misyachniy
сообщение Nov 30 2014, 14:52
Сообщение #1


Знающий
****

Группа: Свой
Сообщений: 716
Регистрация: 27-05-05
Из: Kyiv
Пользователь №: 5 454



Я использовал winfilter для генерации исходника
http://www.winfilter.20m.com/

Код вставлял в Builder без изменений (FIR 4).

Исходно идет оцифровка сигнала 17300 Гц 4-х кратной чатотой.
Затем вычисляются sin/cos и суммирование за 64 периода

Код
balance_value =  ADC_GetInjectedConversionValue(ADC3, ADC_InjectedChannel_1);
  switch (i & 3)
  {
   case 0: Re += balance_value; break;
   case 1: Im += balance_value; break;
   case 2: Re -= balance_value; break;
   case 3: Im -= balance_value; break;
  }  
  
  i++;
  
  if (i == MAX_AVERAGE_BALANCE_SIGNAL )
  {
   i  = 0;
    
   Re /= (MAX_AVERAGE_BALANCE_SIGNAL/4);
   Im /= (MAX_AVERAGE_BALANCE_SIGNAL/4);

    re_value = Re;
    im_value = Im;

    amp_series_samples = (Re * Re + Im * Im);
  }


Чтобы результат (Re * Re + Im * Im) поместился в 32 бита накопленные значение делятся на 4.
Но на расчет фильтра, по моему, не должны влиять

Для расчета фильтра нужна частота семплирования и верхняя граница пропускания.

Частоту семплирования я рассчитал как 17300 * 4 / 64 = 1082Гц.
Частоту верхней границы ставил от 1 до 10Гц но "красивого" гладкого сигнала не получил.
меньше 1 Гц программа не позволяет использовать
Пробовал 1082Гц увеличивать/уменьшать 2 и 4 раза - не помогло.

Пример плохого сигнала "small_fe.PNG", хорошего "big_cu.PNG", сдвоеного "1_cop_ua.PNG" в архиве.
В нем же csv данные.
Одну картинку прикладываю отделльно для оперативного просмотра



Эскизы прикрепленных изображений
Прикрепленное изображение
 

Прикрепленные файлы
Прикрепленный файл  _________________________________.rar ( 133.88 килобайт ) Кол-во скачиваний: 10
 
Go to the top of the page
 
+Quote Post
Maverick
сообщение Dec 1 2014, 11:40
Сообщение #2


я только учусь...
******

Группа: Модераторы
Сообщений: 3 447
Регистрация: 29-01-07
Из: Украина
Пользователь №: 24 839



Цитата(misyachniy @ Nov 30 2014, 16:52) *
Я использовал winfilter для генерации исходника
http://www.winfilter.20m.com/

Код вставлял в Builder без изменений (FIR 4).

Исходно идет оцифровка сигнала 17300 Гц 4-х кратной чатотой.
Затем вычисляются sin/cos и суммирование за 64 периода

Код
balance_value =  ADC_GetInjectedConversionValue(ADC3, ADC_InjectedChannel_1);
  switch (i & 3)
  {
   case 0: Re += balance_value; break;
   case 1: Im += balance_value; break;
   case 2: Re -= balance_value; break;
   case 3: Im -= balance_value; break;
  }  
  
  i++;
  
  if (i == MAX_AVERAGE_BALANCE_SIGNAL )
  {
   i  = 0;
    
   Re /= (MAX_AVERAGE_BALANCE_SIGNAL/4);
   Im /= (MAX_AVERAGE_BALANCE_SIGNAL/4);

    re_value = Re;
    im_value = Im;

    amp_series_samples = (Re * Re + Im * Im);
  }


Чтобы результат (Re * Re + Im * Im) поместился в 32 бита накопленные значение делятся на 4.
Но на расчет фильтра, по моему, не должны влиять

Для расчета фильтра нужна частота семплирования и верхняя граница пропускания.

Частоту семплирования я рассчитал как 17300 * 4 / 64 = 1082Гц.
Частоту верхней границы ставил от 1 до 10Гц но "красивого" гладкого сигнала не получил.
меньше 1 Гц программа не позволяет использовать
Пробовал 1082Гц увеличивать/уменьшать 2 и 4 раза - не помогло.

Пример плохого сигнала "small_fe.PNG", хорошего "big_cu.PNG", сдвоеного "1_cop_ua.PNG" в архиве.
В нем же csv данные.
Одну картинку прикладываю отделльно для оперативного просмотра

в матлабе не проще подобрать фильтр? например с помощью fdatool


--------------------
If it doesn't work in simulation, it won't work on the board.

"Ты живешь в своих поступках, а не в теле. Ты — это твои действия, и нет другого тебя" Антуан де Сент-Экзюпери повесть "Маленький принц"
Go to the top of the page
 
+Quote Post
Lmx2315
сообщение Dec 1 2014, 11:48
Сообщение #3


отэц
*****

Группа: Свой
Сообщений: 1 729
Регистрация: 18-09-05
Из: Москва
Пользователь №: 8 684



QUOTE (Maverick @ Dec 1 2014, 14:40) *
в матлабе не проще подобрать фильтр? например с помощью fdatool

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


--------------------
b4edbc0f854dda469460aa1aa a5ba2bd36cbe9d4bc8f92179f 8f3fec5d9da7f0
SHA-256
Go to the top of the page
 
+Quote Post
misyachniy
сообщение Dec 4 2014, 17:13
Сообщение #4


Знающий
****

Группа: Свой
Сообщений: 716
Регистрация: 27-05-05
Из: Kyiv
Пользователь №: 5 454



Попробовал "чистый" matlab.
Фильтр Баттерворта первого порядка дает отличный результат файл "matlab butter 1 order.PNG"
Судя по запаздыванию Matlab соединил около сотни одиночных. Это для меня много.

При использовании заготовки программы winfilt08 с 7-tap файл "winfilt08.PNG"
Результат не отличается от простого усредняющего

Код
int filter(short NewSample)
{

// Фильтр с сайта EasyElectronics.ru
//  Y(n) = (15*Y(n-1) + X(n)) >> 4

y = y*15 + NewSample >> 4;
return y;

return fir(NewSample);
}

Рисунок отфильтрованного сигнала "simple_average.PNG"

Эскизы прикрепленных изображений
Прикрепленное изображение
Прикрепленное изображение
Прикрепленное изображение
 
Go to the top of the page
 
+Quote Post
rx9cim
сообщение Dec 25 2014, 16:56
Сообщение #5


Участник
*

Группа: Участник
Сообщений: 33
Регистрация: 2-07-12
Пользователь №: 72 593



Посмотрите прогу rxdisp для SDR, там есть исходники, в том числе нужные вам, правда для float.
Go to the top of the page
 
+Quote Post
misyachniy
сообщение Jan 2 2015, 11:56
Сообщение #6


Знающий
****

Группа: Свой
Сообщений: 716
Регистрация: 27-05-05
Из: Kyiv
Пользователь №: 5 454



Наиболее хорошо получился результат с такой прогаммой расчета коэффициентов.
http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html

Код
//-----------------------------------------------------
#define NZEROS 2
#define NPOLES 2

// bessel 2 order  Fs =423, Fc = 5
float gain = 4.854200800e+02;
float coef_0 =  -0.8489920322;
float coef_1 =  1.8407517468;

float filter_bessel_5Hz_AMP(float input_value)  // 11.6 uS
  {
      static float xv[NZEROS+1], yv[NPOLES+1];  
       xv[0] = xv[1]; xv[1] = xv[2];
        xv[2] = input_value / gain;
                                            
        yv[0] = yv[1]; yv[1] = yv[2];
        
        yv[2] = xv[0] + xv[2] + 2 * xv[1]
                      + (coef_0 * yv[0]) + (coef_1 * yv[1]);

        return  yv[2];
  }

Исходник сгенерированный программой пришлось перекроить и вынести коэффициенты в переменную.
cortex m3 имеет систему команд в 16 бит шириной и с константами работает медленнее.( ИМХО :-))

Результат на рисунке:
Прикрепленное изображение


На верхнем графике амплитуда и фаза до фильтрации
На нижнем амплитуда после фильтрации.

Возник вопрос по определению наличия сигнала.
Вот увеличенный отфильтрованный сигнал:
Прикрепленное изображение


Явно виден тренд к повышению уровня постоянной составляющей.
Если ее отсечь, то определить пики сигнал легко.

Порылся в интернете по слову "тренд" - предлагают фильтры для торговли на бирже.
Или низкочастотные или текущее среднее предлагают.
Есть ли другие фильтры для поиска сигнала? Как они называются?

Go to the top of the page
 
+Quote Post
ViKo
сообщение Jan 2 2015, 13:34
Сообщение #7


Универсальный солдатик
******

Группа: Модераторы
Сообщений: 8 634
Регистрация: 1-11-05
Из: Минск
Пользователь №: 10 362



Цитата(misyachniy @ Jan 2 2015, 14:56) *
Явно виден тренд к повышению уровня постоянной составляющей.
Если ее отсечь, то определить пики сигнал легко.

Порылся в интернете по слову "тренд" - предлагают фильтры для торговли на бирже.
Или низкочастотные или текущее среднее предлагают.
Есть ли другие фильтры для поиска сигнала? Как они называются?

ФВЧ? biggrin.gif
Go to the top of the page
 
+Quote Post
misyachniy
сообщение Jan 2 2015, 18:11
Сообщение #8


Знающий
****

Группа: Свой
Сообщений: 716
Регистрация: 27-05-05
Из: Kyiv
Пользователь №: 5 454



Цитата(ViKo @ Jan 2 2015, 15:34) *
ФВЧ? biggrin.gif

Бесселя ФВЧ не бывает.
Другие фильтры дадут выбросы.
Go to the top of the page
 
+Quote Post
ViKo
сообщение Jan 3 2015, 07:50
Сообщение #9


Универсальный солдатик
******

Группа: Модераторы
Сообщений: 8 634
Регистрация: 1-11-05
Из: Минск
Пользователь №: 10 362



Цитата(misyachniy @ Jan 2 2015, 21:11) *
Бесселя ФВЧ не бывает.
Другие фильтры дадут выбросы.

CR - цепочка.
Go to the top of the page
 
+Quote Post
des00
сообщение Jan 3 2015, 17:02
Сообщение #10


Вечный ламер
******

Группа: Модераторы
Сообщений: 7 248
Регистрация: 18-03-05
Из: Томск
Пользователь №: 3 453



Извините за глупый вопрос, а откуда вообще берется постоянка, да еще и фиксированным трендом? Во входном сигнале ее не видно.


--------------------
Go to the top of the page
 
+Quote Post
misyachniy
сообщение Jan 4 2015, 11:06
Сообщение #11


Знающий
****

Группа: Свой
Сообщений: 716
Регистрация: 27-05-05
Из: Kyiv
Пользователь №: 5 454



Цитата(des00 @ Jan 3 2015, 19:02) *
Извините за глупый вопрос, а откуда вообще берется постоянка, да еще и фиксированным трендом? Во входном сигнале ее не видно.


Перед АЦП спользуется усилитель с однополярным питанием и виртуальной землей.
Прикрепленное изображение


Тренд не фиксированый.
Это нагрев схемы и/или датчика.
Go to the top of the page
 
+Quote Post
ViKo
сообщение Jan 4 2015, 11:41
Сообщение #12


Универсальный солдатик
******

Группа: Модераторы
Сообщений: 8 634
Регистрация: 1-11-05
Из: Минск
Пользователь №: 10 362



Что вам мешает вместо (вместе с) R8 поставить конденсатор?
Go to the top of the page
 
+Quote Post
misyachniy
сообщение Jan 4 2015, 12:54
Сообщение #13


Знающий
****

Группа: Свой
Сообщений: 716
Регистрация: 27-05-05
Из: Kyiv
Пользователь №: 5 454



Цитата(ViKo @ Jan 4 2015, 13:41) *
Что вам мешает вместо (вместе с) R8 поставить конденсатор?


Поставить то я могу.
Но пливет опорное напряжение АЦП - в процессоре нет выделеного вывода опорного.
К тому же отрицательные выбросы будут на входе АЦП.
Go to the top of the page
 
+Quote Post
ViKo
сообщение Jan 4 2015, 12:58
Сообщение #14


Универсальный солдатик
******

Группа: Модераторы
Сообщений: 8 634
Регистрация: 1-11-05
Из: Минск
Пользователь №: 10 362



Цитата(misyachniy @ Jan 4 2015, 15:54) *
Поставить то я могу.
Но пливет опорное напряжение АЦП - в процессоре нет выделеного вывода опорного.
К тому же отрицательные выбросы будут на входе АЦП.

Сместите в середину диапазона АЦП, с помощью резисторного делителя. Как сделали на входе ОУ.
Опорное напряжение у АЦП вряд ли плывет так сильно. Это выход ОУ плывет.

Кстати, что за ОУ, какой у него входной ток, какое падение создадут входные токи на тех резисторах, что стоят?
Go to the top of the page
 
+Quote Post
des00
сообщение Jan 4 2015, 15:18
Сообщение #15


Вечный ламер
******

Группа: Модераторы
Сообщений: 7 248
Регистрация: 18-03-05
Из: Томск
Пользователь №: 3 453



Цитата(misyachniy @ Jan 4 2015, 18:06) *
Тренд не фиксированый.
Это нагрев схемы и/или датчика.

банальный корректор постоянки : измерение узким фильтром уровня постоянки и вычитание его из сигнала. Ну а затем пиковый детектор на скользящем окне наблюдения?


--------------------
Go to the top of the page
 
+Quote Post
Hose
сообщение Jan 8 2015, 12:08
Сообщение #16


Частый гость
**

Группа: Участник
Сообщений: 82
Регистрация: 7-01-15
Пользователь №: 84 450



Интегрирование за 64 определит частотную характеристику, умножение на син и кос это перенос спектра, перед переносом достаточно почистить зеркальные каналы. Если бы знать все характеристики сигнала и задачи.... Можно было бы найти решение проще.
Go to the top of the page
 
+Quote Post
misyachniy
сообщение Jan 12 2015, 16:49
Сообщение #17


Знающий
****

Группа: Свой
Сообщений: 716
Регистрация: 27-05-05
Из: Kyiv
Пользователь №: 5 454



Цитата(des00 @ Jan 4 2015, 17:18) *
банальный корректор постоянки : измерение узким фильтром уровня постоянки и вычитание его из сигнала. Ну а затем пиковый детектор на скользящем окне наблюдения?


Это слишком "жирно" у меня STM32F103 на 72МГц.

Простое решение это нужно усреднять значение сигнала за какой-то период и смотреть на отклонение сигнала от усредненного.
Вопрос в том как есть ли методы нахождения периода усреднения и порога для принятия решения.


Для начала я написал функцию «скользящего среднего»
Код
// адаптивная подсторойка "нулевого" уровня
void adaptive_calc_average(int new_sample, statistic *stat, MD_PARAMETER *md_parameter)
{
int difference = stat->average - new_sample;

if (abs(difference) < md_parameter->amp_difference )
{
   stat->average *=31;
   stat->average +=new_sample;
   stat->average /=32;
}
}


Предполагалось, что период усреднения можно будет адаптировать в зависимости от сигнала. Пока реализован только механизм отброса заведомо больших отклонений.
if (abs(difference) < md_parameter->amp_difference ) , такие отклонения в расчет не попадают.

Затем была написана функция определения объекта по отклонению сигнала от текущего среднего.

Код
//---------------------------------------------------
// определить наличие объекта по статистике сигнала в покое
int pulse_width_detecting(int new_sample, statistic *stat, MD_PARAMETER *md_parameter)
{

int difference, abs_difference;

adaptive_calc_average(new_sample, stat, md_parameter);

difference = stat->average - new_sample;

if (abs(difference) < md_parameter->amp_difference) // сигнал не превысил порог
  {
    pulse_width =0;
    pulse_sign =0;
    return 0;
  }

  if (pulse_sign ==0) // первое превышение порога
  {
   if (difference > 0) pulse_sign = 1;
   if (difference < 0) pulse_sign = -1;
   return 0;
  }

  if ((difference <0) && (pulse_sign < 0)) // отрицательный импульс
  {
    pulse_width++;

    if (pulse_width > md_parameter->amp_len) return -1;
    else return 0;
  }

  if ((difference > 0) && (pulse_sign >0)) // положительный импульс
  {
    pulse_width++;
    if (pulse_width > md_parameter->amp_len) return 1;
    else return 0;
  }

  pulse_width=0;
  pulse_sign =0;
  return 0;
}


Определяется длина непрерывного отклонения сигнала от текущего среднего, превышающее порог.

Решил привязать значение порога к статистике сигнала – среднее/среднеквадратичное отклонения. Для простоты взял 128 (int len) значений сигнала

//---------------------------------------------------
// подсчитать статитику фонового сигнала
void calc_stat(int *data, int len, statistic *stat)
{
int i, max, min;
float sum, sum2;

max = min = data[0];
sum =0;
for (i=0; i<len; i++)
{
if (max < data[i]) max = data[i];
if (min > data[i]) min = data[i];
sum += data[i];
}

stat->average = sum/len; // среднее значение

max = abs(max - stat->average);
min = abs(stat->average - min);

if (max > min) stat->max_diff = max; // максимальный выброс был в плюс
else stat->max_diff = min; // максимальный выброс был в минус

sum = sum2 =0;
for (i=0; i<len; i++)
{
sum += abs(stat->average - data[i]);
sum2 += abs(stat->average - data[i]) * abs(stat->average - data[i]);
}

stat->mean_diff = sum/len;
stat->sco = sqrt(sum2)/len;
}

Прогнал 4 образца Cu, Al, Fe, Ni и подобрал значение порога ширины импульса.

Прикрепленное изображение


красный – сигнал;
синий – текущее среднее;
зеленый – функция детектирования сигнала.

Вроде бы все хорошо. Решил сменить условия эксперимента.

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

Прикрепленное изображение


Взаимосвязь среднего отклонения прослеживается, при изменении среднего отклонения с 0,2188 до 0,2422 появились ложные цели при тех же параметрах.

Решил проверить зависимость по "старым" измерениям. Собрал статистику о том во сколько раз порог должен превышать значение среднего отклонения. Получилось число от 5 до 64 – то есть корреляции нет. :-(

В аналогичном приборе автор отказался от адаптивной подстройки, и порог выбирается вручную. Оно как бы и логично при пороге 2 и 3 все работает.
У меня в теории, у автора в практике :-)

Есть фильтры Калмана о них пол интернета исписано, но толкового описания нет, да и STM32F103 его не потянет.
Есть фильтры которыми выявляют отклик от эхолотов и радиолокаторов. Но они тоже тяжеловесные.

Вот и думаю, есть проще типовое решение адаптивной подстройки порога или нет?
Go to the top of the page
 
+Quote Post
des00
сообщение Jan 12 2015, 17:18
Сообщение #18


Вечный ламер
******

Группа: Модераторы
Сообщений: 7 248
Регистрация: 18-03-05
Из: Томск
Пользователь №: 3 453



Цитата(misyachniy @ Jan 13 2015, 00:49) *
Это слишком "жирно" у меня STM32F103 на 72МГц.

dc(k) = a * adcin(k) + (1-a)*dc(k-1)

a = 1/2^16;

adcout(k) = adcin(k) - dc(k)

вот это не потянет ? на фиксированной точке в формате s1.30 ?


--------------------
Go to the top of the page
 
+Quote Post
misyachniy
сообщение Jan 12 2015, 19:30
Сообщение #19


Знающий
****

Группа: Свой
Сообщений: 716
Регистрация: 27-05-05
Из: Kyiv
Пользователь №: 5 454



Цитата(des00 @ Jan 12 2015, 19:18) *
dc(k) = a * adcin(k) + (1-a)*dc(k-1)

a = 1/2^16;

adcout(k) = adcin(k) - dc(k)

вот это не потянет ? на фиксированной точке в формате s1.30 ?


По моему это один из вариантов скользящего среднего.

При установившемся режиме работать будет.
Сейчас по включению датчик до 30 минут имеет довольно быстрый разогрев.
Вот хотелось бы простой понятный алгоритм определения параметров усреднения в приведенном случае коэффициента "a".

В аналоговых схемах применяют диференцирующие цепочки.
В матлабе есть диференцирующие фильтры.
Но мне не понравилась работа FDA.
Точнее не достаточно привычки и/или опыта.
Go to the top of the page
 
+Quote Post
des00
сообщение Jan 12 2015, 19:42
Сообщение #20


Вечный ламер
******

Группа: Модераторы
Сообщений: 7 248
Регистрация: 18-03-05
Из: Томск
Пользователь №: 3 453



Цитата(misyachniy @ Jan 13 2015, 02:30) *
По моему это один из вариантов скользящего среднего.

Шарите sm.gif он самый, но в рекурсивной форме, причем абсолютно устойчивый, с минимальным ресурсом и без децимации по выходу. Хотя децимацию можно добавить. Другое название экспоненциальный фильтр.

Цитата
Вот хотелось бы простой понятный алгоритм определения параметров усреднения в приведенном случае коэффициента "a".

АЧХ построите и сразу увидите. а = 2^16 эквивалентно скользящему среднему на 65536 отсчетов.


--------------------
Go to the top of the page
 
+Quote Post

2 страниц V   1 2 >
Reply to this topicStart new topic
1 чел. читают эту тему (гостей: 1, скрытых пользователей: 0)
Пользователей: 0

 


RSS Текстовая версия Сейчас: 20th July 2025 - 19:01
Рейтинг@Mail.ru


Страница сгенерированна за 0.03873 секунд с 7
ELECTRONIX ©2004-2016