Фильтр Гаусса от А до Я

Все что касается фильтрации
Pavia
Сообщения: 38
Зарегистрирован: 25 апр 2011, 18:45

Фильтр Гаусса от А до Я

Сообщение Pavia »

Давно хочу реализовать фильтр Гаусса при помощи фильтра с БИХ. Содрать я могу, это не проблема. А вот досконально разобраться для меня важнее. Так как по мимо Гауссовского фильтра надо будет проектировать и другие фильтры. Из всего ЦОС я специализируюсь на обработки изображений.

Так теперь к вопросам.

Вопрос первый.
(1)
Как известно можно приближенно описать дробью смотри ссылку [1]

Так вот делая преобразование Лапласа.

Где
s=j*w
j- мнимая единица
a1=a3=a5=0

Так вот передаточная функцию легко построить перебирая частоты.
Но если промоделировать. Возьмем функцию Дирака.
То есть заполним массив нулями и один отсчёт по середине взять не нулевым, я беру 256 но это не суть важно.

Как известно передаточная функция и БИХ фильтр связан следующим преобразованием.
Если рассматривается передаточная функция вида:
Изображение
то соотношение между входом и выходом такой системы должно удовлетворять разностному уравнению:
Изображение

То можно получить что фильтр не работает. Но если произвести разделение как предложено в работе [1] то все работает. Почему это так? С чем связанна "неустойчивость"? Не пойму, где ломается теория?



[1] http://homepage.tudelft.nl/e3q6n/public ... 95TYLV.pdf
В статье [1] есть неточность, связанная с неправильным расчётом граничных условий. Об этом написано в
[2] http://www.cwp.mines.edu/Meetings/Project06/cwp546.pdf

PS. tex подглючиват степень в дроби плохо рисуется
PPS. Следующий вопрос будет завтра. А то этот вопрос я и так устал набирать.

Аватара пользователя
Бахурин Сергей
Администратор
Сообщения: 1114
Зарегистрирован: 05 окт 2010, 19:55
Контактная информация:

Re: Фильтр Гаусса от А до Я

Сообщение Бахурин Сергей »

Насколько я понял при беглом ознакомлении со статьями, то фильтр Гаусса имеет передаточную характеристику . Тогда все полюса расположены в левой полуплоскости плоскости s и фильтр устойчив. Соответственно для перехода в z плоскость надо произвести билинейное преобразование . В этом смысле никакая теория не ломается все фильтры так и рассчитываются только выделяя полюса из левой полуплоскости (подробно в этом разделе). Если же вы в фильтр подставите и полюса в правой полуплоскости, то фильтр будет неустойчивым. Поэтому факторизация передаточной характеристики на и необходима для выделения устойчивого фильтра.

Pavia
Сообщения: 38
Зарегистрирован: 25 апр 2011, 18:45

Re: Фильтр Гаусса от А до Я

Сообщение Pavia »

Вопрос второй.

Решил не затягивать. Есть реализация Гаусса через БИХ.
К примеру от Intel
http://origin-software.intel.com/en-us/ ... xtensions/

Так вот к статье прилагается код.
http://software.intel.com/file/34862

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

void calGaussianCoeff( float sigma,  float *a0, float *a1, float *a2, float *a3, float *b1, float *b2, float *cprev, float *cnext)
{
  float alpha, lamma,  k; 
  // defensive check
  if (sigma < 0.5f)
	  sigma = 0.5f;

  alpha = (float) exp((0.726)*(0.726)) / sigma;
  lamma = (float)exp(-alpha);
  *b2 = (float)exp(-2*alpha);
  k = (1-lamma)*(1-lamma)/(1+2*alpha*lamma- (*b2));
  *a0 = k;
  *a1 = k*(alpha-1)*lamma;
  *a2 = k*(alpha+1)*lamma;
  *a3 = -k* (*b2);
  *b1 = -2*lamma;
  *cprev = (*a0 + *a1)/(1+ *b1 + *b2); 
  *cnext = (*a2 + *a3)/(1+ *b1 + *b2);
}
Поясните откуда берутся формулы cprev cnext? А да что-то я таких коэффициентов ни в одной книги их не видел. :-(
А также откуда взялась константа 0.726 ?
А да если кто может, то напишите по какому методу тут идет расчёт?

Есть еще от ATI/AMD коэффициенты другие, но принцип похожий.

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

RecursiveGaussian::computeGaussParms(float fSigma, int iOrder, GaussParms* pGP)
{
    // pre-compute filter coefficients
    pGP->nsigma = fSigma; // note: fSigma is range-checked and clamped >= 0.1f upstream
    pGP->alpha = 1.695f / pGP->nsigma;
    pGP->ema = exp(-pGP->alpha);
    pGP->ema2 = exp(-2.0f * pGP->alpha);
    pGP->b1 = -2.0f * pGP->ema;
    pGP->b2 = pGP->ema2;
    pGP->a0 = 0.0f;
    pGP->a1 = 0.0f;
    pGP->a2 = 0.0f;
    pGP->a3 = 0.0f;
    pGP->coefp = 0.0f;
    pGP->coefn = 0.0f;

    switch (iOrder) 
    {
    case 0: 
        {
            const float k = (1.0f - pGP->ema)*(1.0f - pGP->ema)/(1.0f + (2.0f * pGP->alpha * pGP->ema) - pGP->ema2);
            pGP->a0 = k;
            pGP->a1 = k * (pGP->alpha - 1.0f) * pGP->ema;
            pGP->a2 = k * (pGP->alpha + 1.0f) * pGP->ema;
            pGP->a3 = -k * pGP->ema2;
        } 
        break;
    case 1: 
        {
            pGP->a0 = (1.0f - pGP->ema) * (1.0f - pGP->ema);
            pGP->a1 = 0.0f;
            pGP->a2 = -pGP->a0;
            pGP->a3 = 0.0f;
        } 
        break;
    case 2: 
        {
            const float ea = exp(-pGP->alpha);
            const float k = -(pGP->ema2 - 1.0f)/(2.0f * pGP->alpha * pGP->ema);
            float kn = -2.0f * (-1.0f + (3.0f * ea) - (3.0f * ea * ea) + (ea * ea * ea));
            kn /= (((3.0f * ea) + 1.0f + (3.0f * ea * ea) + (ea * ea * ea)));
            pGP->a0 = kn;
            pGP->a1 = -kn * (1.0f + (k * pGP->alpha)) * pGP->ema;
            pGP->a2 = kn * (1.0f - (k * pGP->alpha)) * pGP->ema;
            pGP->a3 = -kn * pGP->ema2;
        } 
        break;
    default:
        // note: iOrder is range-checked and clamped to 0-2 upstream
        return;
    }
    pGP->coefp = (pGP->a0 + pGP->a1)/(1.0f + pGP->b1 + pGP->b2);
    pGP->coefn = (pGP->a2 + pGP->a3)/(1.0f + pGP->b1 + pGP->b2);
}

Аватара пользователя
Бахурин Сергей
Администратор
Сообщения: 1114
Зарегистрирован: 05 окт 2010, 19:55
Контактная информация:

Re: Фильтр Гаусса от А до Я

Сообщение Бахурин Сергей »

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

mr.bit
Сообщения: 50
Зарегистрирован: 19 окт 2010, 13:08

Re: Фильтр Гаусса от А до Я

Сообщение mr.bit »

Pavia писал(а):Давно хочу реализовать фильтр Гаусса при помощи фильтра с БИХ.
А у Вас нет исходников расчета коэффициентов фильтра Гаусса для КИХ структуры? Спасибо.

Аватара пользователя
Бахурин Сергей
Администратор
Сообщения: 1114
Зарегистрирован: 05 окт 2010, 19:55
Контактная информация:

Re: Фильтр Гаусса от А до Я

Сообщение Бахурин Сергей »

У меня есть:

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


int filterGaussian(double* g, int Lg, double BT, int bc){   
   if((!g) || (Lg<=1) || (BT < 0.0))
      return 0;
   double dt = (double)bc * 2.0/(double)(Lg-1);   
   double t;
   double beta  = 28.477658649975010867721351422732 * BT * BT;
   t = -(double)bc-dt;
   double sg = 0.0;
   for(int i = 0; i<Lg; i++){
      t += dt;
      g[i] = exp(-beta*t*t);
      sg+=g[i];
   }
   for(int i = 0; i<Lg; i++){
      g[i] /= sg;
   }   
   return 1;
}

Обсуждался как то в этой ветке

Pavia
Сообщения: 38
Зарегистрирован: 25 апр 2011, 18:45

Re: Фильтр Гаусса от А до Я

Сообщение Pavia »

Я лично таким пользуюсь.

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

// Расчет коэффициентов Гауссовского фильтра для реализации через КИХ фильтр
procedure Get_Gausian_FIR_Coeficient(var a:TArrayReal; Sigma:Real);
var i,r:Integer;
s:Real;
begin
s:=0;
r:=round(Sigma * 3*sqrt(2*pi)/4);
SetLength(a,2*r+1);
for i:=0 to Length(a)-1 do
  begin
  a[i]:=1/(sqrt(2*pi)*Sigma)*Exp(-sqr(i-r)/(2*Sigma*Sigma));
  s:=s+a[i];
  end;
for i:=0 to Length(a)-1 do
  a[i]:=a[i]/s;
end;

mr.bit
Сообщения: 50
Зарегистрирован: 19 окт 2010, 13:08

Re: Фильтр Гаусса от А до Я

Сообщение mr.bit »

Бахурин Сергей
Pavia
Спасибо большое :)

BVU
Сообщения: 4
Зарегистрирован: 08 июн 2011, 12:03

Re: Фильтр Гаусса от А до Я

Сообщение BVU »

Добрый день!
Вопрос следующий: в dsp.dll библиотеке которая свободно распространяется с этого сайта так же имеется реализация расчета коэффициентов фильтра Гауса. Замечу еще раз коэффициентов фильтра! Которые в графическом виде представляют собой импульсную характеристику фильтра Гауса (Гаусова кривая - колокол)... А где же сам процесс фильтрации сигнала? Значит, так такового фильтра Гауса в библиотеке - НЕТ! Как же программно организовать процесс фильтрации?!!!

Аватара пользователя
Бахурин Сергей
Администратор
Сообщения: 1114
Зарегистрирован: 05 окт 2010, 19:55
Контактная информация:

Re: Фильтр Гаусса от А до Я

Сообщение Бахурин Сергей »

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

Ответить