Расчет коэффициентов БИХ(IIR) фильтра.

[Kider_Dim]

Добрый день!
 В настоящее время ушел в обработку дискретных сигналов (аудиозаписей).
 Никак не могу разобраться с тем, как рсчитать коэффициенты для фильтров с бесконечной импульсной характеристикой (для КИХ фильтров все просто и понятно).
 В Matlabе все просто и понятно - ввел начальные параметры (частоты пропускания, среза, дисретизации, уровень подавления/искажения и т. д.) и получил Хэдер с рассчетными нумератором и денумератором (по Мурзилкам если - коэффициентами перед z, которые находытся в числителе и знаменателе передаточной функции).
 Но мне необходимо, чтобы программа без участия Matlaba с какой-либо стороны в зависимости от входных параметров рассчитывала необходиме коэффициенты и, соответственно, фильтровала входной сигнал по разностной схеме фильтра.

 В настоящее время разобрался и реализовал рассчет аналогового фильтра Баттерворта нижних частот до момента рсчета его импульсной характеристики. Не могу никак дойти до того, как рассчитать значения для какждой из секций фильтра в виде Num, Denum. Судя по теории на каждую секцию приходится по 3 их значения (по три коэффициента).
 Сталкивался ли кто с такой задачей? Помогите, пожалуйста, понять до конца принцип расчета такого рода фильтров (реализовать расчет, если понимание будет, итак смогу).

 Литературу долго искал с вменяемым описанием математики и принципов и остановился на таких источниках:

http://www.dsplib.ru/content.html
http://www.mikroe.com/chapters/view/...-filters/#id32

Ниже привожу код функции (расчет фильтра Баттерворта), которую пока для расчета накропал. Только до рассчета передаточной функции:

Code: [Select]

double Filter_IIR_Coeffs(long Filter_Type,long Filter,long Filter_Order,long Sample_Rate,long Fs,long Fp, double Rp,double Rs)
//Функция рассчета IIR фильтров
//http://www.mikroe.com/chapters/view/73/chapter-3-iir-filters/
{
 //Индексы Filter_Type:
  //0 - Баттерворт
  //1- Чебышев первого порядка
  //2- Чебышев второго порядка
 //Индексы Filter:
  //0 - LPF
  //1 - HPF
  //2 - BPF
  //3 - BSF
 //Filter_Order - порядок фильтра
 //Sample_Rate - частота дискретизации
 //Fp - частота пропускания
 //Fs - частота затухания
 //Rp - частота искажения в полосе пропускания, Дб
 //Rs - частота подавления в полосе затухания,  Дб
 if (Filter == 0)
  //Считаю LPF фильтр
  {
   if (Filter_Type == 0)
//Фильтр Баттерворта нижних частот рассчитываю.
{
  //Расчет аналогового прототипа фильтра. Считаю константу H0.
float H0p = tan(M_PI*(Fp/(float)Sample_Rate));
float H0s = tan(M_PI*(Fs/(float)Sample_Rate));
//Рассчитываю левые полюса(S) аналогового фильтра.
  float *S = (float*)malloc(sizeof(float)*2*(Filter_Order)); //массив значений полюсов.столбик 0 - левый полюс, столбик 1 - правый столбик
  float *Sk = (float*)malloc(sizeof(float)*(Filter_Order));//Массив значений Sk левых полюсов. По сути - это значение полюса в целом
for (long j = 0; j < Filter_Order; j++)
{
  S[j,(0*Filter_Order)] = cosl(M_PI*(0.5+(2*(j)+1)/(float)(2*Filter_Order)));//левая сторона полюса
  S[j,(1*Filter_Order)] = sinl(M_PI*(0.5+(2*(j)+1)/(float)(2*Filter_Order)));//правая сторона полюса
}
for (long i = 0; i < Filter_Order; i++)
  {
   Sk[i] = exp(((-1)*M_PI)*(0.5+((2*i+1)/(float)(2*Filter_Order))));
  }
//Рассчет передаточной характеристики аналогового фильтра.
  long M;
   M = Filter_Order;//количество нулевых полюсов, рассчитанное для аналогового фильтра.
  float Analog_Hs;
   Analog_Hs = (Sk[0]/(float)H0p) + H0p*(S[0,(0*Filter_Order)]+(-1)*S[0,(1*Filter_Order)]);
   for (long i = 1; i < M; i++)
{
Analog_Hs *= Sk[i] + (S[i,(0*Filter_Order)]+(-1)*S[i,(1*Filter_Order)]);
}
   Analog_Hs = 1/(float)Analog_Hs;
   Analog_Hs = pow(H0p,2)*Analog_Hs;
  //----------------------------------------------------
//Закончил Фильтр Баттерворта нижних частот расчитывать
}
  }

return true;
}