Фильтр Гаусса от А до Я
Фильтр Гаусса от А до Я
Давно хочу реализовать фильтр Гаусса при помощи фильтра с БИХ. Содрать я могу, это не проблема. А вот досконально разобраться для меня важнее. Так как по мимо Гауссовского фильтра надо будет проектировать и другие фильтры. Из всего ЦОС я специализируюсь на обработки изображений.
Так теперь к вопросам.
Вопрос первый.
(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. Следующий вопрос будет завтра. А то этот вопрос я и так устал набирать.
Так теперь к вопросам.
Вопрос первый.
(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. Следующий вопрос будет завтра. А то этот вопрос я и так устал набирать.
- Бахурин Сергей
- Администратор
- Сообщения: 1116
- Зарегистрирован: 05 окт 2010, 19:55
- Контактная информация:
Re: Фильтр Гаусса от А до Я
Насколько я понял при беглом ознакомлении со статьями, то фильтр Гаусса имеет передаточную характеристику . Тогда все полюса расположены в левой полуплоскости плоскости s и фильтр устойчив. Соответственно для перехода в z плоскость надо произвести билинейное преобразование . В этом смысле никакая теория не ломается все фильтры так и рассчитываются только выделяя полюса из левой полуплоскости (подробно в этом разделе). Если же вы в фильтр подставите и полюса в правой полуплоскости, то фильтр будет неустойчивым. Поэтому факторизация передаточной характеристики на и необходима для выделения устойчивого фильтра.
Re: Фильтр Гаусса от А до Я
Вопрос второй.
Решил не затягивать. Есть реализация Гаусса через БИХ.
К примеру от Intel
http://origin-software.intel.com/en-us/ ... xtensions/
Так вот к статье прилагается код.
http://software.intel.com/file/34862
Поясните откуда берутся формулы cprev cnext? А да что-то я таких коэффициентов ни в одной книги их не видел.
А также откуда взялась константа 0.726 ?
А да если кто может, то напишите по какому методу тут идет расчёт?
Есть еще от ATI/AMD коэффициенты другие, но принцип похожий.
Решил не затягивать. Есть реализация Гаусса через БИХ.
К примеру от 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);
}
А также откуда взялась константа 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);
}
- Бахурин Сергей
- Администратор
- Сообщения: 1116
- Зарегистрирован: 05 окт 2010, 19:55
- Контактная информация:
Re: Фильтр Гаусса от А до Я
К сожалению, пока не могу вам ответить на этот вопрос надо разбираться что это за магические числа в коде.
Re: Фильтр Гаусса от А до Я
А у Вас нет исходников расчета коэффициентов фильтра Гаусса для КИХ структуры? Спасибо.Pavia писал(а):Давно хочу реализовать фильтр Гаусса при помощи фильтра с БИХ.
- Бахурин Сергей
- Администратор
- Сообщения: 1116
- Зарегистрирован: 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;
}
Re: Фильтр Гаусса от А до Я
Я лично таким пользуюсь.
Код: Выделить всё
// Расчет коэффициентов Гауссовского фильтра для реализации через КИХ фильтр
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;
Re: Фильтр Гаусса от А до Я
Бахурин Сергей
Pavia
Спасибо большое
Pavia
Спасибо большое
Re: Фильтр Гаусса от А до Я
Добрый день!
Вопрос следующий: в dsp.dll библиотеке которая свободно распространяется с этого сайта так же имеется реализация расчета коэффициентов фильтра Гауса. Замечу еще раз коэффициентов фильтра! Которые в графическом виде представляют собой импульсную характеристику фильтра Гауса (Гаусова кривая - колокол)... А где же сам процесс фильтрации сигнала? Значит, так такового фильтра Гауса в библиотеке - НЕТ! Как же программно организовать процесс фильтрации?!!!
Вопрос следующий: в dsp.dll библиотеке которая свободно распространяется с этого сайта так же имеется реализация расчета коэффициентов фильтра Гауса. Замечу еще раз коэффициентов фильтра! Которые в графическом виде представляют собой импульсную характеристику фильтра Гауса (Гаусова кривая - колокол)... А где же сам процесс фильтрации сигнала? Значит, так такового фильтра Гауса в библиотеке - НЕТ! Как же программно организовать процесс фильтрации?!!!
- Бахурин Сергей
- Администратор
- Сообщения: 1116
- Зарегистрирован: 05 окт 2010, 19:55
- Контактная информация:
Re: Фильтр Гаусса от А до Я
что значит нет фильтра гаусса? К-ты фильтра есть, значит и фильтр есть. Программно фильтрация осуществляется функцией filterFIR, где в качестве импульсной характеристики нужно подставить те к-ты что вернула функция расчета фильтра гаусса.