12 / 12 / 0
Регистрация: 20.11.2013
Сообщений: 167
1

Вычисление коэффициентов БИХ фильтров

10.02.2014, 09:23. Показов 17139. Ответов 26

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

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

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

http://www.dsplib.ru/content.html
http://www.mikroe.com/chapters... ters/#id32

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

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
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 = 1; j < Filter_Order+1; j++)
         {
          S[j-1,(0*Filter_Order)] = cosl(M_PI*(0.5+(2*(j)+1)/(float)(2*Filter_Order)));//левый полюс
          S[j-1,(1*Filter_Order)] = sinl(M_PI*(0.5+(2*(j)+1)/(float)(2*Filter_Order)));//правый полюс
          Sk[j-1] = exp(((-1)*M_PI)*(0.5+((2*j+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]/(float)H0p) + H0p*(S[i,(0*Filter_Order)]+(-1)*S[i,(1*Filter_Order)]);
        }
       Analog_Hs = 1/(float)Analog_Hs;
       Analog_Hs = pow(H0p,2)*Analog_Hs;
      float H0 = H0p;
//Расчет импульсной характеристики цифрового фильтра. 
     float Hz;
      Hz = 0;
      float P1,P2,P3,P4; //произведения в формуле
      //Zk - правые стороны полюсов.....вещественные их части.
      P1 = 1 - (-1)* S[0,(1*Filter_Order)];
      P2 = 1 - S[0,(0*Filter_Order)]-(-1)*S[0,(1*Filter_Order)];
      P3 = 1 - ((1+(((-1)* S[0,(1*Filter_Order)])/(float)(1-((-1)* S[0,(1*Filter_Order)])))));
      P4 = 1 - ((1+S[0,(0*Filter_Order)]+(-1)*S[0,(1*Filter_Order)])/(float)(1-S[0,(0*Filter_Order)]-(-1)*S[0,(1*Filter_Order)]));
      for (long i = 1;  i < Filter_Order; i++)
       {
        P1 *= 1 - (-1)* S[i,(1*Filter_Order)];
        P2 *= 1 - S[i,(0*Filter_Order)]-(-1)*S[i,(1*Filter_Order)];
        P3 *= 1 - ((1+(((-1)* S[i,(1*Filter_Order)])/(float)(1-((-1)* S[i,(1*Filter_Order)])))));
        P4 *= 1 - ((1+S[0,(0*Filter_Order)]+(-1)*S[0,(1*Filter_Order)])/(float)(1-S[0,(0*Filter_Order)]-(-1)*S[0,(1*Filter_Order)]));
       }
      Hz = pow((-1),(Filter_Order-M))*(P1/(float)P2)*pow(1,(Filter_Order-M))*(P3/(float)P4);
      //----------------------------------------------------
     //Закончил Фильтр Баттерворта нижних частот расчитывать
    }
  }
 
return true;
}
__________________
Помощь в написании контрольных, курсовых и дипломных работ, диссертаций здесь
2
Programming
Эксперт
94731 / 64177 / 26122
Регистрация: 12.04.2006
Сообщений: 116,782
10.02.2014, 09:23
Ответы с готовыми решениями:

Структуры КИХ и БИХ фильтров
Подскажите как мне реализовать КИХ и БИХ фильтр, например, на C++. Допустим у меня есть значения...

Нормировка коэффициентов БИХ-фильтра
Здравствуйте все! У меня есть БИХ ФВЧ. Фильтр второго порядка, коэффициенты: A1 =...

Синтез и программная реализация цифровых фильтров в классе КИХ и БИХ цепей
Синтез и программная реализация цифровых фильтров в классе КИХ и БИХ цепей. Тхнические требования:...

Вычисление биномиальных коэффициентов
честно говоря не очень поняла как это делается. прошу помощи. от чего вообще оттолкнуться?? и куда...

26
0 / 0 / 0
Регистрация: 17.05.2016
Сообщений: 5
12.02.2014, 10:51 2
Есть функциональные зависимости АЧХ и ФЧХ на определенное число точек (частота от минус пи до пи),как по ним найти коэффициенты фильтра различного порядка (или определенного, если задать изменение по N)?
по какой формуле или подскажите где об этом почитать? везде пишут, что это ОДПФ..но не пойму как считается((
0
10226 / 6606 / 496
Регистрация: 28.12.2010
Сообщений: 21,160
Записей в блоге: 1
12.02.2014, 20:47 3
http://www.dsplib.ru/content/filters/fir/fir.html
0
12 / 12 / 0
Регистрация: 20.11.2013
Сообщений: 167
15.02.2014, 16:36  [ТС] 4
Дядьки, без обид.......Но очень разочаровываюсь в форуме....Я надеялся, что поотзывчевей все тут как то......

В общем, сам спросил - сам начинаю отвечать.

Для расчета коэффициентов цифрового фильтра нужно сначала расчитать коэффициенты передаточной функции аналогового прототипа. затем нормировать их на частоту среза..........
Я пока расчитываю фильтр Баттерворта нижних частот по мурзилкам..........Вот список:
http://www.dsplib.ru/content/f... terex.html
http://www.dsplib.ru/content/filters/ch3/ch3.html
http://www.dsplib.ru/content/filters/ch1/ch1.html

Просто ОГРОМНОЕ человеческое СПАСИБО хозяину сайта за такое доступное представление материала!!!!

Всего то надо было выполнить расчет основных параметров фильтра и его коэффициентов в знаменателе (полиномы N степени, коэффициенты которых, как оказалось, можно при помощи линейной свертки выдернуть).
Для меня все затормозили пробелы в математике.......Полупрограмми ст я - мед. кибернетик.

Ниже привожу код расчета.....Буду ОЧЕНЬ признателен, если вдруг кто косяки в нем найдет и мне на них укажет - исправлюсь

Код навырост - функция будет отрабатывать в зависимости от желаемого типа фильтра и его подтипа:

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
float Analog_CoeFFs(long Filter_Type,long Filter,long Fs,long Fp, float Rp,float Rs)
//Функция рассчета IIR фильтров
//http://www.dsplib.ru/content/filters/butterex/butterex.html
{
if (Filter == 0)
 {
  if (Filter_Type == 0)
   {
    float OmegaP = Fp/(float)(2*M_PI); //Нормированная частота пропускания
    float OmegaS = Fs/(float)(2*M_PI); //Нормальзованная частота заграждения
    float TettaP = tan(OmegaP*M_PI/(float)2);
     if (TettaP < 0) TettaP = (-1)*TettaP;
    float TettaS = tan(OmegaS*M_PI/(float)2);
     if (TettaS < 0) TettaS = (-1)*TettaS;
    float EpsP = sqrt(pow(10,Rp/(float)10)-1);
    float EpsS = sqrt(pow(10,Rs/(float)10)-1);
    long N_Grad = log(EpsS/(float)EpsP)/log(TettaS/(float)TettaP);//порядок фильтра
     if (N_Grad < 0) N_Grad =N_Grad*(-1);
    //Если порядок фильтра четный , то L = N/2
    //Нечетный - L = (N - 1)/2
    long L;
     if (N_Grad%2 == 0) L = N_Grad/(float)2;
      else L = (N_Grad-1)/(float)2;
    float Alfa = 1/pow(EpsP,1/(float)N_Grad);
    //Углы Этта считаю. количество зависит от порядка фильтра. Поэтому придется массивчик динамический создать.
     float *Etta = (float*)malloc(sizeof(float)*L);
    for (long i = 0; i < L; i++)
     {
      Etta[i] = (2*(i+1)-1)/(float)(2*N_Grad);
      Etta[i] = Etta[i]*M_PI;
     }
    //Считаю значение перед S второго слагаемого полинома
     float *Second_Zna4 = (float*)malloc(sizeof(float)*L);
     for (long i = 0; i < L; i++)
     {
      Second_Zna4[i] = 2*Alfa*sin(Etta[i]);
     }
    //Числитель (Num) для аналогового фильтра считаю по формуле 1/EpsP.
     float Num_Analog = 1/(float)EpsP;
    //Считаю коэффициенты полинома.
    float Num_CoeffS;
    if (N_Grad == 1) Num_CoeffS = 2;
    if (N_Grad == 2) Num_CoeffS = 3;
    if (N_Grad > 2) Num_CoeffS = N_Grad*2 - 2;
    float *CoeffS = (float*)malloc(sizeof(float)*Num_CoeffS);
    for (long i = 0; i < Num_CoeffS; i++)
     {
      CoeffS[i] = 0;
     }
    Linear_Convolution(0,N_Grad,L,Alfa,Second_Zna4,CoeffS);
    //Нормируем все значения коэффициентов на частоту среза EpsP. Нормирую только коэффициенты при S.
    //т. е. все элементы с начала массива кроме последнего
    for (long i = 0; i < Num_CoeffS-1; i++)
     {
      CoeffS[i] = CoeffS[i]/pow(TettaP,(N_Grad)-i);
     }
    //теперь надо сократить значения всех коэффициентов знаменателя на числитель
    for (long i = 0; i < Num_CoeffS; i++)
     {
      CoeffS[i] = CoeffS[i]/(float)Num_Analog;
      ShowMessage(CoeffS[i]);
     }
    free(Etta);
    free(Second_Zna4);
    return true;
    //Конец рассчета аналогового фильтра Баттерворта
   }
 }
}

Ну и функция для вычисления коэффициентов полинома.....Тут условий всяких продумать много пришлось....

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
float Linear_Convolution(long Polinom_Type,long N,long L,float Alfa,float *Sec_Zna4,float *CoeffS)
//линейна свертка(получение коэффициентов полиномов)
//Polinom_Type - тип полинома(Баттерворт, Чебышев)
//N - порядок фильтра
//Sec_Zna4 - массив со значенями при синусе в полиноме
//CoeffS - результирующий массив коэффициентов полинома
{
 if (Polinom_Type == 0) //полином Баттерворта
  {
    long r; //степень для (S+Alfa)
     if (N%2 == 0) r = 0;
      else r = 1;
    //Сами коэффициенты считаю. сразу делаем код независимым от порядка фильтра.
    //по сути - надо сначала рассчитать коэффициенты полинома под П и, если r позволяет, рассчитать результат его домножение на получившийся результат
    if (N == 1) //считать нехер - (S+Alfa) результат,коэффициенты = [1,Alfa];
     {
      long CoeffS_Kol = 2;
      CoeffS[0] = 1;
      CoeffS[1] = Alfa;
     return true;
     }
     if (N == 2) //та же ситуация - (S^2+Second_Zna4*S+Alfa^2)
     {
      long CoeffS_Kol = 3;
      CoeffS[0] = 1;
      CoeffS[1] = Sec_Zna4[0];
      CoeffS[2] = Alfa*Alfa;
     return true;
     }
    //Дальше уже сложнее, но мы не ищем легких путей - нам надо все унифицировать ;))
    if (N > 2)
     {
      float *Second_Mirrow = (float*)malloc(sizeof(float)*3); //зеркальный вектор (S^2+Second_Zna4*S+Alfa^2). Для линейной свертки готовлю
      long Size;
       Size = 0;
      for (long i = 0; i < L; i++)
       {
        if (i == 0)
         {
          Size = 3;
           CoeffS[0] = 1;
           CoeffS[1] = Sec_Zna4[i];
           CoeffS[2] = Alfa*Alfa;
         }
          else
           //свертка с целью получения коэффициентов результирующих
           {
            free(Second_Mirrow);
            float *Second_Mirrow = (float*)malloc(sizeof(float)*Size);
            for (long p = 0; p < Size; p++)
             {
              Second_Mirrow[p] = 0;
             }
            //Массив Second_Mirrow заполняю. Постоянно одинаковый для каждого последующего L будет.
             Second_Mirrow[0] = Alfa*Alfa;
             Second_Mirrow[1] = Sec_Zna4[i];
             Second_Mirrow[2] = 1;
 
            float *First = (float*)malloc(sizeof(float)*Size);
             //заполняю массив First
            for (long j = 0; j < Size; j++)
             {
              First[j] = CoeffS[j];
             }
            //Свертка. Тут уже должен последовательно перемножить и сложить элементы массива First с элементами массива Second_Mirrow.
              for (long j = 0; j < Size; ++j)
               {
                for (long t = 0; t < 3; ++t)
                 {
                  CoeffS[j + t] += First[j] * Second_Mirrow[t];
                 }
               }
            free(First);
            Size = Size + 1;
           }
       }
       //Теперь, если r == 1 надо домножить получившуюся последовательность на (S+Alfa)
        free(Second_Mirrow);
        if (r == 1)
         {
          Second_Mirrow = (float*)malloc(sizeof(float)*Size);
          Second_Mirrow[0] = 0;
          Second_Mirrow[1] = 1;
          for (long p = 2; p < Size; p++)
           {
            Second_Mirrow[p] = 0;
           }
 
          float *First = (float*)malloc(sizeof(float)*Size);
             //заполняю массив First
            for (long j = 0; j < Size; j++)
             {
              First[j] = CoeffS[j];
             }
            //Свертка. Тут уже должен последовательно перемножить и сложить элементы массива First с элементами массива Second_Mirrow.
              for (long j = 0; j < Size; ++j)
               {
                for (long t = 0; t < Size; ++t)
                 {
                  CoeffS[j + t] += First[j] * Second_Mirrow[t];
                 }
               }
            free(First);
            free(Second_Mirrow);
       }
     return true;
     }
 //Закончил считать коэффициенты полинома Баттерворта
  }
return true;
}
Воот.....надеюсь, что кому то даже все это полезно будет, т. к. я для себя выявил неприятную тенденцию - Всем проще Хедер вытащить с Матлаба с определенными коэффициентами под конкретный фильтр, чем что то универсальное сварганить.....
1
12 / 12 / 0
Регистрация: 20.11.2013
Сообщений: 167
17.02.2014, 12:13  [ТС] 5
Таак - накосячил чуток - неправильно все работало.....Сейчас вроде все должно нормально быть:

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
void Analog_CoeFFs(long Filter_Type,long Filter,long Fs,long Fp, float Rp,float Rs,float *CoeffS)
//Функция рассчета IIR фильтров
//http://www.dsplib.ru/content/filters/butterex/butterex.html
{
if (Filter == 0)
 {
  if (Filter_Type == 0)
   {
    float OmegaP = Fp/(float)(2*M_PI); //Нормированная частота пропускания
    float OmegaS = Fs/(float)(2*M_PI); //Нормальзованная частота заграждения
    float TettaP = tan(OmegaP*M_PI/(float)2);
     if (TettaP < 0) TettaP = (-1)*TettaP;
    float TettaS = tan(OmegaS*M_PI/(float)2);
     if (TettaS < 0) TettaS = (-1)*TettaS;
    float EpsP = sqrt(pow(10,Rp/(float)10)-1);
    float EpsS = sqrt(pow(10,Rs/(float)10)-1);
    long N_Grad = log(EpsS/(float)EpsP)/log(TettaS/(float)TettaP);//порядок фильтра
     if (N_Grad < 0) N_Grad =N_Grad*(-1);
    //Если порядок фильтра четный , то L = N/2
    //Нечетный - L = (N - 1)/2
    long L;
     if (N_Grad%2 == 0) L = N_Grad/(float)2;
      else L = (N_Grad-1)/(float)2;
    float Alfa = 1/pow(EpsP,1/(float)N_Grad);
    //Углы Этта считаю. количество зависит от порядка фильтра. Поэтому придется массивчик динамический создать.
     float *Etta = (float*)malloc(sizeof(float)*L);
    for (long i = 0; i < L; i++)
     {
      Etta[i] = (2*(i+1)-1)/(float)(2*N_Grad);
      Etta[i] = Etta[i]*M_PI;
     }
    //Считаю значение перед S второго слагаемого полинома
     float *Second_Zna4 = (float*)malloc(sizeof(float)*L);
     for (long i = 0; i < L; i++)
     {
      Second_Zna4[i] = 2*Alfa*sin(Etta[i]);
     }
    //Числитель (Num) для аналогового фильтра считаю по формуле 1/EpsP.
     float Num_Analog = 1/(float)EpsP;
    //Считаю коэффициенты полинома.
    long Num_CoeffS;
    if (N_Grad == 1) Num_CoeffS = 2;
    if (N_Grad == 2) Num_CoeffS = 3;
    if (N_Grad > 2) Num_CoeffS = N_Grad*2 - 2;
    for (long i = 0; i < Num_CoeffS; i++)
     {
      CoeffS[i] = 0;
     }
//считаю коэффициенты полинома
//Сначала под произведением коэффициенты считаю.
//Формирую треугольную матрицу размером L+1 строк и L+1 столбцов.
float *Analog_Den = (float*)malloc(sizeof(float)*Num_CoeffS*Num_CoeffS);
  //Первый элемент = 1 (полином в нулевой степени), остальные обнуляю.
  for (long i = 0; i < Num_CoeffS; i++)
   {
    for (long j = 0; j < Num_CoeffS; j++)
     {
      if ((i == 0) && (j == 0))    Analog_Den[j+(i*Num_CoeffS)] = 1;
       else  Analog_Den[j+(i*Num_CoeffS)] = 0;
     }
   }
float *Vector = (float*)malloc(sizeof(float)*Num_CoeffS);
for (long i = 0; i < Num_CoeffS; i++) 
 {
  Vector[i] = 0;    
 }
for (long i = 0; i < L; i++) 
 {
  Vector[0] = 1;
  Vector[1] = Second_Zna4[i];
  Vector[2] = Alfa*Alfa;
  Matrix_Linear_Conv_First(Num_CoeffS,Analog_Den,Vector);   
 }
 if (N_Grad %2 !=0) 
  {
   for (long i = 0; i < Num_CoeffS; i++)
   {
    Vector[i] = 0;  
   }
   Vector[0] = 1;
   Vector[1] = Alfa;
   Matrix_Linear_Conv_First(Num_CoeffS,Analog_Den,Vector);   
  } 
//Теперь вытаскиваю последнюю строку матрицы. Это и будут ненормированные коэффициенты фильтра аналогового.
for (long i = 0; i < Num_CoeffS; i++) 
 {
  CoeffS[i] = Analog_Den[(Num_CoeffS-1)+(i*Num_CoeffS)];
 } 
//Теперь вытаскиваю последнюю строку матрицы, которая и является решением.
    //Нормируем все значения коэффициентов на частоту среза EpsP. Нормирую только коэффициенты при S.
    //т. е. все элементы с начала массива кроме последнего
    for (long i = 0; i < Num_CoeffS-1; i++)
     {
      CoeffS[i] = CoeffS[i]/pow(TettaP,(Num_CoeffS-1)-i);
     }
    //теперь надо сократить значения всех коэффициентов знаменателя на числитель
    for (long i = 0; i < Num_CoeffS; i++)
     {
      CoeffS[i] = CoeffS[i]/(float)Num_Analog;
     }
    free(Etta);
    free(Second_Zna4);
    free(Vector);
    free(Analog_Den);
    //Конец рассчета аналогового фильтра Баттерворта
   }
 }
}


C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
void Matrix_Linear_Conv_First(long Num_El,float *Matrix,float *Vector)
//расчет матриц Matrix1, которые являются результатом перемножения каждого из столюцоы на вектор Matrix2
{
 float *First = (float*)malloc(sizeof(float)*Num_El); //Собственно - значения строки матрицы Matrix1, которую домножаю на Second
 float *Result_Conv = (float*)malloc(sizeof(float)*Num_El); //промeжуточный массив в памяти. Чтобы не запутаться с двумерным массивом......
 for (long j = 0; j < Num_El; j++)
  {
   Result_Conv[j] = 0;
  }
 for (long i = 0; i < Num_El-1; i++)
  {
   for (long j = 0; j < Num_El; j++)
    {
     First[j] = 0;
     First[j] = Matrix[i+(j*Num_El)];
    }
   //Свертка
   Result_Conv[0] = First[0]*Vector[0];
   for (long h = 1; h < Num_El; h++)
    {
     for (long t = h; t > 0; t--)
      {
       Result_Conv[h] += First[abs(t-h)]*Vector[t];
      }
    }
   for (long k = 0; k < Num_El; k++)
    {
     Matrix[(i+1)+(k*Num_El)] = Result_Conv[k];
    }
  }
  free(First);
  free(Result_Conv);
}
1
алкокодер
157 / 153 / 41
Регистрация: 27.12.2012
Сообщений: 550
17.02.2014, 13:29 6
Цитата Сообщение от DimKaKiber Посмотреть сообщение
Дядьки, без обид.......Но очень разочаровываюсь в форуме....Я надеялся, что поотзывчевей все тут как то......
Форум уже не торт и тема специфична.
0
12 / 12 / 0
Регистрация: 20.11.2013
Сообщений: 167
18.02.2014, 14:39  [ТС] 7
Ну да - тема специфичная......но я почему то думал, что тема натоптанная........
Коэффициенты цифрового фильтра в итоге я расчитал так:
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
void Digital_Filter_Calc(long Polinom_Type,long N,float *A1,float *B1,float *Digital_Coeffs_Num,float *Digital_Coeffs_DeNum)
//Расчет цифрового фильтра по его аналоговому прототипу
//Polinom_Type - тип полинома(Баттерворт, Чебышев)
//N - порядок фильтра
//*А - вектор коэффициентов числителя аналогового фильтра
//*В - вектор коэффициентов знаменателя аналогового фильтра
//*Digital_Coeffs_Num - коэффициенты числителя
//*Digital_Coeffs_DeNum - коэффициенты знаменателя
{
if (Polinom_Type == 0)//фильтр Баттерворта
 {
  long Num_El = N+1; //количество строк и столбцов матриц подстановки
  float *NUM = (float*)malloc(sizeof(float)*Num_El*Num_El);//Num_El строк, Num_El столбцов
 
  //Первый элемент матриц Denum и NUM = 1 (полином в нулевой степени), остальные обнуляю.
  for (long i = 0; i < Num_El; i++)
   {
    for (long j = 0; j < Num_El; j++)
     {
      if ((i == 0) && (j == 0))    NUM[j+(i*Num_El)] = 1;
       else  NUM[j+(i*Num_El)] = 0;
     }
   }
   //Получаю матрицу Num. Надо домножить каждую строку матрицы на вектор Бетта = [1,1].
   float *Betta = (float*)malloc(sizeof(float)*Num_El);
   for (long i = 0; i < Num_El; i++)
    {
     Betta[i] = 0;
    }
   Betta[1] = 1;
   Betta[0] = 1;
 Matrix_Linear_Conv_First(Num_El,NUM,Betta);
   //Получаю матрицу DeNum. Надо домножить каждую строку матрицы на вектор Альфа = [1,-1].
   float *DeN = (float*)malloc(sizeof(float)*Num_El*Num_El);//Num_El строк, Num_El столбцов
   for (long i = 0; i < Num_El; i++)
   {
    for (long j = 0; j < Num_El; j++)
     {
      if ((i == 0) && (j == 0))    DeN[j+(i*Num_El)] = 1;
       else  DeN[j+(i*Num_El)] = 0;
     }
   }
   float *Den = (float*)malloc(sizeof(float)*Num_El);
   for (long i = 0; i < Num_El; i++)
    {
     Den[i] = 0;
    }
   Den[0] = 1;
   Den[1] = -1;
   Matrix_Linear_Conv_First(Num_El,DeN,Den);
//Получаю матрицу ND.
//Перемножать надо: ND(i) = Den(N-i+2)*Num(i)
float *ND = (float*)malloc(sizeof(float)*Num_El*Num_El);
for (long i = 0; i < Num_El; i++)
 {
  for (long j = 0; j < Num_El; j++)
   {
    ND[j+(i*Num_El)] = 0;
   }
 }
float *A = (float*)malloc(sizeof(float)*Num_El);//строка Den(N-i+2)
float *B = (float*)malloc(sizeof(float)*Num_El);//строка Num(i)
float *Out = (float*)malloc(sizeof(float)*Num_El);//промежуточный массив результатов
for (long i = 0; i < Num_El; i++)
 {
  for (long j = 0; j < Num_El; j++)
   {
    A[j] = 0;
    B[j] = 0;
    Out[j] = 0;
   }
  for (long j = 0; j < Num_El; j++)
   {
    A[j] = DeN[(Num_El-1-i)+(j*Num_El)];
    B[j] = NUM[i+(j*Num_El)];
   }
  Vector_Conv(Num_El, A,B,Out);
  for (long j = 0; j < Num_El; j++)
   {
    ND[i+(j*Num_El)] = Out[j];
   }
 }
//Теперь получившиеся значения надо помножить на коэффициенты числителя и знаменателя передаточной функции аналогового фильтра.
float *NDN = (float*)malloc(sizeof(float)*Num_El*Num_El); //Результирующая матрица
for (long i = 0; i < Num_El; i++)
 {
  for (long j = 0; j < Num_El; j++)
   {
    NDN[j+(i*Num_El)] = 0;
   }
 }
//Тупо каждый элемент матрицы ND домножаю на i-й элемент вектора A
for (long i = 0; i < Num_El; i++)
 {
  for (long j = 0; j < Num_El; j++)
   {
    NDN[i+(j*Num_El)] = ND[i+(j*Num_El)]*A1[i];
   }
 }
float *NDD = (float*)malloc(sizeof(float)*Num_El*Num_El); //Результирующая матрица
for (long i = 0; i < Num_El; i++)
 {
  for (long j = 0; j < Num_El; j++)
   {
    NDD[j+(i*Num_El)] = 0;
   }
 }
//Тупо каждый элемент матрицы ND домножаю на i-й элемент векторов с коэффициентамичислителя и знаменателя
for (long i = 0; i < Num_El; i++)
 {
  for (long j = 0; j < Num_El; j++)
   {
    NDD[i+(j*Num_El)] = ND[i+(j*Num_El)]*B1[i];
   }
 }
//теперь получаю коэффициенты. Они рассчитываются как сумма коэффициентов по столбцам матриц  NDD и NDN.
//Коэффициены числителя.
for (long i = 0; i < Num_El; i++)
 {
  for (long j = 0; j < Num_El; j++)
   {
    Digital_Coeffs_Num[i] += NDN[j+(i*Num_El)];
   }
 }
//Коэффициенты знаменателя.
for (long i = 0; i < Num_El; i++)
 {
  for (long j = 0; j < Num_El; j++)
   {
    Digital_Coeffs_DeNum[i] += NDD[j+(i*Num_El)];
   }
 }
float Diver;
 Diver = Digital_Coeffs_DeNum[0];
//Теперь приводим все расчитанные значения к коэффициенту знаменателя при нулевой степени S.
for (long i = 0; i < Num_El; i++)
 {
  Digital_Coeffs_Num[i] = Digital_Coeffs_Num[i]/Diver;
 }
for (long i = 0; i < Num_El; i++)
 {
  Digital_Coeffs_DeNum[i] = Digital_Coeffs_DeNum[i]/Diver;
 }
 
  free(Betta);
  free(Den);
  free(NUM);
  free(Betta);
  free(DeN);
  free(Den);
  free(ND);
  free(A);
  free(B);
  free(Out);
  free(NDN);
  free(NDD);
  //Закончил цифровой фильтр Баттерворта расчитывать
 }
}
C++
1
2
3
4
5
6
7
8
9
10
11
void Vector_Conv(long Size, float *A,float *B,float *Out)
{
   Out[0] = A[0]*B[0];
   for (long h = 1; h < Size; h++)
    {
     for (long t = h; t > -1; t--)
      {
       Out[h] += A[abs(t-h)]*B[t];
      }
    }
}

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

Тут я привел расчеты фильтра Баттерворта нижних частот. Сейчас никак не могу правильно реализовать саму процедуру фильтрации. В теории до конца разберусь - сделаю, думаю


Ну и ошибки, если в коде найдете (что, к сожалению, весьма вероятно) - просьба написать о них.....будем исправлять.....Но вроде операции пошагово отладил и на бумажке руками расчеты сделал....Все сошлось
0
алкокодер
157 / 153 / 41
Регистрация: 27.12.2012
Сообщений: 550
20.02.2014, 06:32 8
DimKaKiber, я думаю тема достойна оформления в блог.
1
12 / 12 / 0
Регистрация: 20.11.2013
Сообщений: 167
20.02.2014, 06:43  [ТС] 9
Пока нет времени на оформление - с головой в решении задачи, т. к. помимо реализации фильтров ндо еще и экспериментальные исследования провести в части выявления их оптимальных параметров для достижения поставленной передо мной цели, но приятно видеть такое мнение)
К сожалению почему то желаемого результата при фильтрации не достигаю........КИХ фильтр намного лучше все делает и качественней (именно то, что от него ожидается). Сейчас буду пытаться ошибки в рассчетах отлавливать (либо в реализации самой процедуры фильтрации, т. к. КИХ составляющая моего расчетного фильтра (числитель) работает верно).
0
10226 / 6606 / 496
Регистрация: 28.12.2010
Сообщений: 21,160
Записей в блоге: 1
20.02.2014, 09:01 10
Проблема не в форуме, а в отсутствии желания TC запостить специфичную тему там где положено. Кто ее тут видит? Специалисты по ЦОС? Так это не раздел ЦОС или Матлаба, шансы ее увидеть только просматривая все новые темы сразу. А ну да, пишите на С++, а вопрос-то совсем не в С++. Сами виноваты. Вот такие темы и вылавливаешь периодически.

что тема не нова совсем - в сети нет никаких внятных алгоритмоы рассчета коэффициентов фильтров. Да, и если уж совсем честно, и внятной теории маловато
ой, ну не надо плакаться. Море внятных алгоритмов, начиная с уже упомянутого ресурса dsplib, ресурса model.exponenta.ru, заканчивая спецлитературой, к примеру Эммануила С. Айфичера http://books.google.ru/books?i... B2&f=false
0
12 / 12 / 0
Регистрация: 20.11.2013
Сообщений: 167
20.02.2014, 09:17  [ТС] 11
А кто то плакался???
Когда читаешь одно и то же практически во всех учебниках, то и появляется такое мнение. Но это мнение именно, а не жалобы и мольбы о помощи.
+ на форуме подходящего раздела не обнаружил.

За книжку спасибо - почему то не натыкался на неё - погляжу.
Тема не подходит форуму - удаляйте - мне не жалко)
0
10226 / 6606 / 496
Регистрация: 28.12.2010
Сообщений: 21,160
Записей в блоге: 1
20.02.2014, 14:44 12
+ на форуме подходящего раздела не обнаружил.
для особо непонимающих:
https://www.cyberforum.ru/digi... processing
1
12 / 12 / 0
Регистрация: 20.11.2013
Сообщений: 167
21.02.2014, 10:51  [ТС] 13
Спасибо за перенос темы

По ходу дела нашел ошибки в своей реализации.

теперь как в статье все коэффициенты получаются. Но только для фильтров до 3-его порядка (Фильтр Баттерворта).
Появилось достаточно много вопросов:
1. Почему при большем порядке фильтра коэффициент числителя для цифрового фильтра обнуляется?
2. Как верно тогда реализовать расчет при таком раскладе?
3. Как правильно нормировать коэффициенты (в моем случае преобразование ФНЧ-ФНЧ) под мою частоту дискретизации и частоту среза?
4. Не могу понять как правильно в коде реализовать разностную схему БИХ фильтра.

Ниже привожу свой исправленный код расчета для ФНЧ Баттерворта.

Расчет порядка фильтра:

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
float IIR_Gradation(int Filter_Type,long Filter,long Fs,long Fp, float Rp,float Rs)
 {
 if (Filter == 0)
 {
  if (Filter_Type == 0)
   {
   float OmegaP = Fp/(float)(2*M_PI); //Нормированная частота пропускания
     OmegaP = round(OmegaP,0.1);
   float OmegaS = Fs/(float)(2*M_PI); //Нормальзованная частота заграждения
    OmegaS = round(OmegaS,0.1);
    float TettaP = tan((OmegaP*M_PI)/(float)2);
     if (TettaP < 0) TettaP = (-1)*TettaP;
    TettaP = round(TettaP,0.0001);
    float TettaS = tan((OmegaS*M_PI)/(float)2);
     if (TettaS < 0) TettaS = (-1)*TettaS;
    TettaS = round(TettaS,0.0001);
    float EpsP = sqrt(pow(10,Rp/(float)10)-1);
    EpsP = round(EpsP,0.0001);
    float EpsS = sqrt(pow(10,Rs/(float)10)-1);
    EpsS = round(EpsS,0.0001);
    float N = log10(EpsS/(float)EpsP)/log10(OmegaS/(float)OmegaP);
     N *= 10;
      N = floor(N+0.5);
       N /= 10;
    long N_Grad;
     N_Grad = 0;
      N_Grad = N;
     if (N_Grad < 0) N_Grad =N_Grad*(-1);
   return N_Grad;
   }
 }
 }
Расчет аналогового прототипа:

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
void Analog_CoeFFs(long N_Grad,long Filter_Type,long Filter,long Fs,long Fp, float Rp,float Rs,double *CoeffS)
//Функция рассчета IIR фильтров
//http://www.dsplib.ru/content/filters/butterex/butterex.html
{
if (Filter == 0)
 {
  if (Filter_Type == 0)
   {
   double OmegaP = Fp/(float)(2*M_PI); //Нормированная частота пропускания
     OmegaP = round(OmegaP,0.1);
   double OmegaS = Fs/(float)(2*M_PI); //Нормальзованная частота заграждения
    OmegaS = round(OmegaS,0.1);
    double TettaP = tan((OmegaP*M_PI)/(float)2);
     if (TettaP < 0) TettaP = (-1)*TettaP;
    TettaP = round(TettaP,0.0001);
    double TettaS = tan((OmegaS*M_PI)/(float)2);
     if (TettaS < 0) TettaS = (-1)*TettaS;
    TettaS = round(TettaS,0.0001);
    double EpsP = sqrt(pow(10,Rp/(float)10)-1);
    EpsP = round(EpsP,0.0001);
    double EpsS = sqrt(pow(10,Rs/(float)10)-1);
    EpsS = round(EpsS,0.0001);
    //Если порядок фильтра четный , то L = N/2
    //Нечетный - L = (N - 1)/2
    long L;
     if (N_Grad%2 == 0) L = N_Grad/(float)2;
      else L = (N_Grad-1)/(float)2;
    double Alfa = 1/pow(EpsP,(1/(float)N_Grad));
    Alfa = round(Alfa,0.0001);
    //Углы Этта считаю. количество зависит от порядка фильтра. Поэтому придется массивчик динамический создать.
     float *Etta = (float*)malloc(sizeof(float)*L);
    for (long i = 0; i < L; i++)
     {
      Etta[i] = 0;
       Etta[i] = (2*(i+1))-1;
        Etta[i] =  Etta[i] * M_PI;
         Etta[i] = Etta[i]/(float)(2*N_Grad);
     }
    //Считаю значение перед S второго слагаемого полинома
     double *Second_Zna4 = (double*)malloc(sizeof(double)*L);
     for (long i = 0; i < L; i++)
      {
       Second_Zna4[i] = 0;
      }
     for (long i = 0; i < L; i++)
     {
      Second_Zna4[i] = round(2*Alfa*sin(Etta[i]),0.0001);
     }
    //Числитель (Num) для аналогового фильтра считаю по формуле 1/EpsP.
     double Num_Analog = 1/(float)EpsP;
    //Считаю коэффициенты полинома.
    long Num_CoeffS;
    if (N_Grad == 1) Num_CoeffS = 2;
    if (N_Grad == 2) Num_CoeffS = 3;
    if (N_Grad > 2) Num_CoeffS = N_Grad*2 - 2;
    for (long i = 0; i < Num_CoeffS; i++)
     {
      CoeffS[i] = 0;
     }
//Теперь надо рассчитать значения полинома в знаменателе передаточной характеристики аналогового фильтра.
//Тут просто все довольно таки - перемножаю два вектора между собой (линейная свертка) (S^2+Second_Zna4*S+Alfa^2)
//В зависмости от L меняю значения Second_Zna4 вектора.
double *Vector = (double*)malloc(sizeof(double)*Num_CoeffS);
double *Bufer = (double*)malloc(sizeof(double)*Num_CoeffS);//Буфер для расчета коэффициентов
for (long i = 0; i < Num_CoeffS; i++)
//инициализирую массив Vector, чтобы шняги лишней не понасчитать
 {
  Vector[i] = 0;
  Bufer[i] = 0;
 }
for (long i = 0; i < L; i++)
 {
  if (i == 0)
  //Тут по любому всего три коэффициента получается
   {
    Vector[2] = pow(Alfa,2);
    Vector[1] = Second_Zna4[i];
    Vector[0] = 1;
    CoeffS[2] = Vector[2];
    CoeffS[1] = Vector[1];
    CoeffS[0] = Vector[0];
   }
   else
   //Домножаем предыдущий результат на обновленный вектор
    {
     Vector[1] = 0;
     Vector[1] = Second_Zna4[i];
     Vector_Conv(Num_CoeffS,CoeffS,Vector,Bufer);
     //Пишем получившиеся значения в массив CoeffS
     for (long j = 0; j < Num_CoeffS; j++)
      {
       CoeffS[j] = Bufer[j];
      }
    }
 }
 //Теперь, в зависимости от порядка фильтра (четность его или нечетность проверяем),
 //нужно провести еще одну операцию - домножение на (S+Alfa)
if (N_Grad%2 != 0)
 {
  //Обновляем массив Vector.
   for (long i = 0; i < Num_CoeffS; i++)
   //инициализирую массив Vector, чтобы шняги лишней не понасчитать
    {
     Vector[i] = 0;
    }
   Vector[0] = 1;
   Vector[1] = Alfa;
  //домножаем на выражение:
   Vector_Conv(Num_CoeffS,CoeffS,Vector,Bufer);
  //Пишем коэффициенты в массив CoeffS
   for (long i = 0; i < Num_CoeffS; i++)
    {
     CoeffS[i] = Bufer[i];
    }
 }
free(Etta);
free(Second_Zna4);
free(Vector);
free(Bufer);
//Нормируем все значения коэффициентов на частоту среза EpsP. Нормирую только коэффициенты при S.
    //т. е. все элементы с начала массива кроме последнего
    for (long i = 1; i < Num_CoeffS; i++)
     {
      CoeffS[i-1] = CoeffS[i-1]/(float)pow(TettaP,(Num_CoeffS)-i);
     }
    //теперь надо сократить значения всех коэффициентов знаменателя на числитель
    for (long i = 0; i < Num_CoeffS; i++)
     {
      CoeffS[i] = CoeffS[i]/(float)Num_Analog;
     }
    //Конец рассчета аналогового фильтра Баттерворта
   }
 }
}
Расчет цифрового фильтра:

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
void Digital_Filter_Calc(long Polinom_Type,long N,float *A1,double *B1,double *Digital_Coeffs_Num,double *Digital_Coeffs_DeNum)
//Расчет цифрового фильтра по его аналоговому прототипу
//Polinom_Type - тип полинома(Баттерворт, Чебышев)
//N - порядок фильтра
//*А - вектор коэффициентов числителя аналогового фильтра
//*В - вектор коэффициентов знаменателя аналогового фильтра
//*Digital_Coeffs_Num - коэффициенты числителя
//*Digital_Coeffs_DeNum - коэффициенты знаменателя
{
if (Polinom_Type == 0)//фильтр Баттерворта
 {
  long Num_El = N+1; //количество строк и столбцов матриц подстановки
  long *NUM = (long*)malloc(sizeof(long)*Num_El*Num_El);//Num_El строк, Num_El столбцов
 
  //Первый элемент матриц Denum и NUM = 1 (полином в нулевой степени), остальные обнуляю.
  for (long i = 0; i < Num_El; i++)
   {
    for (long j = 0; j < Num_El; j++)
     {
      if ((i == 0) && (j == 0))    NUM[j+(i*Num_El)] = 1;
       else  NUM[j+(i*Num_El)] = 0;
     }
   }
   //Получаю матрицу Num. Надо домножить каждую строку матрицы на вектор Бетта = [1,1].
   long *Betta = (long*)malloc(sizeof(long)*Num_El);
   for (long i = 0; i < Num_El; i++)
    {
     Betta[i] = 0;
    }
   Betta[1] = 1;
   Betta[0] = 1;
 Matrix_Linear_Conv_FirstL(Num_El,NUM,Betta);
   //Получаю матрицу DeNum. Надо домножить каждую строку матрицы на вектор Альфа = [1,-1].
   long *DeN = (long*)malloc(sizeof(long)*Num_El*Num_El);//Num_El строк, Num_El столбцов
   for (long i = 0; i < Num_El; i++)
   {
    for (long j = 0; j < Num_El; j++)
     {
      if ((i == 0) && (j == 0))    DeN[j+(i*Num_El)] = 1;
       else  DeN[j+(i*Num_El)] = 0;
     }
   }
   long *Den = (long*)malloc(sizeof(long)*Num_El);
   for (long i = 0; i < Num_El; i++)
    {
     Den[i] = 0;
    }
   Den[0] = 1;
   Den[1] = -1;
   Matrix_Linear_Conv_FirstL(Num_El,DeN,Den);
//Получаю матрицу ND.
//Перемножать надо: ND(i) = Den(N-i+2)*Num(i)
long *ND = (long*)malloc(sizeof(long)*Num_El*Num_El);
for (long i = 0; i < Num_El; i++)
 {
  for (long j = 0; j < Num_El; j++)
   {
    ND[j+(i*Num_El)] = 0;
   }
 }
long *A = (long*)malloc(sizeof(long)*Num_El);//строка Den(N-i+2)
long *B = (long*)malloc(sizeof(long)*Num_El);//строка Num(i)
long *Out = (long*)malloc(sizeof(long)*Num_El);//промежуточный массив результатов
for (long i = 0; i < Num_El; i++)
 {
  for (long j = 0; j < Num_El; j++)
   {
    A[j] = 0;
    B[j] = 0;
    Out[j] = 0;
   }
  for (long j = 0; j < Num_El; j++)
   {
    A[j] = DeN[(Num_El-1-i)+(j*Num_El)];
    B[j] = NUM[i+(j*Num_El)];
   }
  Vector_ConvL(Num_El, A,B,Out);
  for (long j = 0; j < Num_El; j++)
   {
    ND[i+(j*Num_El)] = Out[j];
   }
 }
//Теперь получившиеся значения надо помножить на коэффициенты числителя и знаменателя передаточной функции аналогового фильтра.
double *NDN = (double*)malloc(sizeof(double)*Num_El*Num_El); //Результирующая матрица
for (long i = 0; i < Num_El; i++)
 {
  for (long j = 0; j < Num_El; j++)
   {
    NDN[j+(i*Num_El)] = 0;
   }
 }
//Тупо каждый элемент матрицы ND домножаю на i-й элемент вектора A
for (long i = 0; i < Num_El; i++)
 {
  for (long j = 0; j < Num_El; j++)
   {
    NDN[i+(j*Num_El)] = ND[i+(j*Num_El)]*A1[i];
   }
 }
double *NDD = (double*)malloc(sizeof(double)*Num_El*Num_El); //Результирующая матрица
for (long i = 0; i < Num_El; i++)
 {
  for (long j = 0; j < Num_El; j++)
   {
    NDD[j+(i*Num_El)] = 0;
   }
 }
//Тупо каждый элемент матрицы ND домножаю на i-й элемент векторов с коэффициентамичислителя и знаменателя
for (long i = 0; i < Num_El; i++)
 {
  for (long j = 0; j < Num_El; j++)
   {
    NDD[i+(j*Num_El)] = ND[i+(j*Num_El)]*B1[i];
   }
 }
//теперь получаю коэффициенты. Они рассчитываются как сумма коэффициентов по столбцам матриц  NDD и NDN.
//Коэффициены числителя.
for (long i = 0; i < Num_El; i++)
 {
  for (long j = 0; j < Num_El; j++)
   {
    Digital_Coeffs_Num[i] += NDN[j+(i*Num_El)];
   }
 }
//Коэффициенты знаменателя.
for (long i = 0; i < Num_El; i++)
 {
  for (long j = 0; j < Num_El; j++)
   {
    Digital_Coeffs_DeNum[i] += NDD[j+(i*Num_El)];
   }
 }
double Diver;
 Diver = Digital_Coeffs_DeNum[0];
//Теперь приводим все расчитанные значения к коэффициенту знаменателя при нулевой степени S.
for (long i = 0; i < Num_El; i++)
 {
  if ((Diver != 0) && (Diver >  -0.01)) Digital_Coeffs_Num[i] = Digital_Coeffs_Num[i]/Diver;
   else break;
 }
for (long i = 0; i < Num_El; i++)
 {
  if ((Diver != 0) && (Diver >  -0.001)) Digital_Coeffs_DeNum[i] = Digital_Coeffs_DeNum[i]/Diver;
   else break;
 }
  free(Betta);
  free(Den);
  free(NUM);
  free(Betta);
  free(DeN);
  free(Den);
  free(ND);
  free(A);
  free(B);
  free(Out);
  free(NDN);
  free(NDD);
  //Закончил цифровой фильтр Баттерворта расчитывать
 }
}

Теперь свертки (для нахождения коэффициентов полинома):

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
void Matrix_Linear_Conv_First(long Num_El,double *Matrix,double *Vector)
//расчет матриц Matrix1, которые являются результатом перемножения каждого из столбцов на вектор Matrix2
{
 double *First = (double*)malloc(sizeof(double)*Num_El); //Собственно - значения строки матрицы Matrix1, которую домножаю на Second
 double *Result_Conv = (double*)malloc(sizeof(double)*Num_El); //промeжуточный массив в памяти. Чтобы не запутаться с двумерным массивом......
 for (long j = 0; j < Num_El; j++)
  {
   Result_Conv[j] = 0;
  }
 for (long i = 0; i < Num_El-1; i++)
  {
   for (long j = 0; j < Num_El; j++)
    {
     First[j] = 0;
     First[j] = Matrix[i+(j*Num_El)];
    }
   //Свертка
   Result_Conv[0] = First[0]*Vector[0];
   for (long h = 1; h < Num_El; h++)
    {
     for (long t = h; t > 0; t--)
      {
       Result_Conv[h] += First[abs(t-h)]*Vector[t];
      }
    }
   for (long k = 0; k < Num_El; k++)
    {
     Matrix[(i+1)+(k*Num_El)] = Result_Conv[k];
    }
  }
  free(First);
  free(Result_Conv);
}
C++
1
2
3
4
5
6
7
8
9
10
11
void Vector_Conv(long Size, double *A,double *B,double *Out)
{
   Out[0] = A[0]*B[0];
   for (long h = 1; h < Size; h++)
    {
     for (long t = h; t > -1; t--)
      {
       Out[h] += A[abs(t-h)]*B[t];
      }
    }
}
Добавлено через 2 часа 0 минут
Реально себя тупым ощущаю (((((

ФНЧ-ФНЧ никак не могу допетрить как делать...........
Пересчитываю частоту пропускания, привязывая её к частоте дискретизации

C++
1
2
3
float Wn;//Пересчитанная частота пропускания.
     Wn = (2*Fp)/(float)SampleRate;
      Wn = tan(M_PI*(Wn/(float)2));
Дальше как корректировать коэффициенты фильтра то??? Такое ощущение, что смотрю в книгу - вижу очень большууууую фигу((((
0
12 / 12 / 0
Регистрация: 20.11.2013
Сообщений: 167
23.04.2014, 17:58  [ТС] 14
Всем доброго времени суток! Отвлекался на другие задачи......Некоторое время назад правда обратно вернулся к теме вычисления коэффициентов БИХ фильтров.

Озарение, что ли появилось - начал реализовывать метод поиска коэффициентов по полюсам и нолям фильтра.
Источник все тот же - http://www.mikroe.com/chapters/view/73/#id32

Все тот же фильтр Баттерворта мучаю
Оказывается все просто, если въехать чуток в тему:
1. Сначала высчитывается порядок фильтра (http://www.dsplib.ru/content/f... terex.html), или руками жестко задается.

2. Затем , исходя из этого порядка, высчитываются полюса фильтра. Это комплексные числа, поэтому надо сразу готовиться к реализации математических операций над ними (ну или готовые библиотеки подключить. Я свое написал, отладил и уверен в правильности работы). Так как я плохо понимаю с чем комплексные числа едят - пригодились разжеванные описания операций с ними - http://www.fxyz.ru/формулы_по_... ных_чисел/ и http://mathprofi.ru/kompleksny... nikov.html.

3. После этого находится нормальзованная частота среза фильтра нашего (TettaC = tan(M_PI*(Fc/SampleRate))). Fc - частота среза, SampleRate = частота дискретизации.

4. Затем считается Аналоговый прототип фильтра проектируемого. Прототип сразу же нормируется на (3). Свернуть надо будет в одно число все, а не хранить коэффициенты аналогового фильтра еще до кучи.

5. Ну и апофеоз - Производим Z-преобразование и расчитываем в итоге коэффициенты цифрового фильтра.


Авторы (книжки по которой делаю все) конечно как могли запутали читателей, но это только на пользу - хоть вникаешь сидишь.

Реализация метода шустрая (хоть я и не старался даже оптимизировать код), но не нравится, что приходится очень широкий коридор частот задавать.

Ниже привожу часть кода (без служебных функций) :

Поиск полюсов ФНЧ Баттерворта:

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
double *PolesOfButter_LPF(long N,long SampleRate,long BP) //Расчет полюсов ФНЧ Баттерворта
{
//Считаем нормированную полосу среза.
 double Tetta = tan(M_PI*(BP/(float)SampleRate));
 double *Poles = (double*)malloc(sizeof(double)*(N*2));//массив, содержащий значение нолей и полюсов фильтра
 long Count;
  Count = -1;
 for (long i = 0; i < N; i++)
  {
   Count++;
   Poles[Count] = Tetta*cos(M_PI*(0.5+((2*i+1)/(float)(2*N))));//Действительная часть
   Count++;
   Poles[Count] =  Tetta*sin(M_PI*(0.5+((2*i+1)/(float)(2*N))));//Мнимая часть.
  }
 return Poles;
Аналоговый прототип ФНЧ Баттерворта (Все, где встречаются слова Complex - операции с комплексными числами, либо объявление типа, который я для них создал):

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
double Analog_Butter_LPF(long N,long BP,long SampleRate,double *Pole)//Расчитываю аналоговый фильтр Баттерворта.
{
 long double Ho;
    double Tetta = tan(M_PI*(BP/(float)SampleRate));
    Ho = pow(Tetta,N)*pow(-1,N);
   //Считаю, что там у меня в аналоговом фильтре получается (Делаю его свертку до конкретного числа)
    TComplex A,B,C;
    for (long i = 0; i < (N*2); i++)
     {
      if (i == 0)
       {
        A.a = 1 - Pole[i];
        i++;
        A.b = Pole[i];
        i++;
        B.a = 1 - Pole[i];
        i++;
        B.b = Pole[i];
        C = Multiply_Complex(A,B);
       }
       else
        {
         A.a = C.a;
         A.b = C.b;
         B.a = 1 - Pole[i];
         i++;
         B.b = Pole[i];
         C = Multiply_Complex(A,B);
        }
     }
    if (fabs(C.b) < 0.000001)
     {
       Ho = Ho/C.a;
     }
     else //Тут, по идее, деление комплексных чисел мутить надо.
     {
       A.a = 0;
       A.a = Ho;
       A.b = 0;
       C = Div_Complex(A,B);
       double Mod;
        Mod = 0;
       Mod = sqrt(pow(C.a,2)+pow(C.b,2));
       Ho = Ho/Mod;
     }
 return Ho;
}
Ну и Z-преобразование (Vector_Conv() - операция линейной свертки):

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
double *Z_Transform(long N,double *Pole) //Z-преобразование
{
 double *Modul = (double*)malloc(sizeof(double)*N);
  long Count;
   Count = -1;
  TComplex A,B,C,D;
  A.a = 0;   A.a = 1.0; A.b = 0 ;
  for (long i = 0; i < N*2; i++)
   {
    //Считываю значение полюса.
    B.a = 0; B.b  = 0;
    B.a = Pole[i];
    i++;
    B.b = Pole[i];
    //Считаю результат в выражении при Z.
    C.a = 0; C.b = 0;
    C = Sum_Complex(A,B);
    D.a = 0; D.b = 0;
    D =  Raznost_Complex(A,B);
    B.a = 0; B.b  = 0;
    B = Div_Complex(C,D);
    //Считаю модуль получившегося комплексного числа.
    Count++;
    Modul[Count] = 0;
    Modul[Count] = sqrt(pow(B.a,2)+pow(B.b,2));
   }
 
 double *First = (double*)malloc(sizeof(double)*(N+1));
 double *Second = (double*)malloc(sizeof(double)*(N+1));//Массивы, которые будут перемножаться между собой.
  double *Result = (double*)malloc(sizeof(double)*(N+1));
 for (long i = 0; i < N; i++)
  {
   First[i] = 0; Second[i] = 0; Result[i] = 0;
  }
 First[0] = 1.0; First[1] = -1*Modul[0];
 Second[0] = 1.0; Second[1] = -1*Modul[1];
 Vector_Conv(N+1,First,Second,Result);
for (long i = 2; i < N; i++) 
 {
  for (long j = 0; j < i; j++) 
   {
    First[i] = 0;
    First[i] = Result[i];  
   }
  Second[0] = 1.0; Second[1] = -1*Modul[i]; 
  Vector_Conv(N+1,First,Second,Result);
 }   
 
 free(Modul);
 free(First);
 free(Second);
 return Result;
}
0
12 / 12 / 0
Регистрация: 20.11.2013
Сообщений: 167
23.04.2014, 18:24  [ТС] 15
Чтобы не быть голословным добавляю Исходный файл, фильтрованный файл и картинки, показывающие АЧХ спроектированного фильтра, вид исходного и отфильтрованного массива данных (Амплитудно-временное представление). Оговорюсь только - поставил себе задачу - выделить биения сердца из исходного сигнала.

Если честно, то пока меня больше устраивает качество работы БИХ фильтров. Да и не такая уж у них и замороченная реализация, как у БИХ. Но все равно - пробывать надо все ) В поисках лучших решений.

Пока есть один очень волнующий меня вопрос..............Реализация БИХ фильтров каким способом лучше??? Т. е. каким способом можно реализовать более эффективные фильтра такого типа??? Через нули и полюса, или же как это сделано в Матлабе (матричные операции, в основном)???
Миниатюры
Вычисление коэффициентов БИХ фильтров   Вычисление коэффициентов БИХ фильтров   Вычисление коэффициентов БИХ фильтров  

Вычисление коэффициентов БИХ фильтров  
Вложения
Тип файла: zip Звуки.zip (449.4 Кб, 49 просмотров)
0
0 / 0 / 0
Регистрация: 08.02.2016
Сообщений: 2
08.02.2016, 17:07 16
Доброго времени суток! Надеюсь тема еще актуальна... Ознакомился с Вашим изложением проблемы.
На данный момент передо мной стоит задача подобная Вашей, касаемо БИХ-фильтра, по аппроксимации Баттерворта, хоть и относится к радиотехнической области. Перед этим форумом сперва наткнулся на этот сайт http://www.dsplib.ru/content/f... terex.html и по нему так же начал считать.
Для определенности мои исходные данные: частота среза 130 кГц, частота заграждения 572 кГц, параметры Rp и Rs для меньшей путаницы взял как в статье 1дБ и 30дБ соответственно. Расчетно получилось, что необходим 3-й порядок фильтра. Реализую я все в Simulink, в общем Matlab. Т.е. хоть и не точно, но имею возможность проверить конечный результат коэффициентов по крайней мере по их порядку.
Вопрос в следующем: Вам достаточно сосчитать коэффициенты программным способом или в ручную все таки расчет получился? Мне просто необходимо рассчитать это вручную, и когда начал считать по этой методике все было вроде бы неплохо до последнего шага с билинейным преобразованием. Хотя уже после частотного преобразования я получил в знаменателе коэффициенты при s порядка -6 степени в силу моих частот. И не знаю как так вышло, но после билинейного преобразования (в его процессе же были эти мутерные разложения полиномов третьей степени, приведения к общем знаменателям и прочие прелести математики) я получил совсем неправильный результат, получился и в числителе и в знаменателе -6 порядок степени при z, хотя должен получается -3 максимум быть и коэффициенты при числителе порядка единиц, а в знаменателе тоже порядка единиц т.к. пришлось просто пренебрегать уже ничтожными множителями, которые вышли после частотного преобразования. В общем хотелось бы узнать, просчитали ли вы это вручную?

Добавлено через 14 минут
DimKaKiber, и прошу прощения, адресовано мое крайнее сообщение Вам
0
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
08.02.2016, 17:44 17
Простой способ расчета коэффициентов БИХ-фильтра, основанный на матрицах z-преобразования, описан здесь
0
0 / 0 / 0
Регистрация: 08.02.2016
Сообщений: 2
16.02.2016, 12:42 18
A_Santik, DimKaKiber, просто интересует метод, описанный именно здесь http://www.dsplib.ru/content/f... terex.html . Метод на матрицах z- преобразования (конкретно, который значится как "II. z- преобразование") выглядит конечно привлекательно из-за своей простоты, но увы или к счастью хотелось бы разобраться с более сложным.Есть несколько нюансов:
1) Пройдясь по просторам интернета заметил, что отличаются формулы по пересчету аналоговой частоты в цифровую, они однотипные, но все отличаются так, что конечный результат разный. Причем разумеется как в книжках, так и в каком то контенте на сайтах.
2) Не совсем понятно, почему в этом методе для расчетов задается нормированная частота заграждения ws для фнч Баттерворта, и уж тем более неравномерность в полосе заграждения, в то время как при задании параметров для подобного фильтра даже в готовых утилитах используется только частота среза и частота дискретизации. Тем более о какой неравномерности в полосе заграждения идет речь, если по мат.части фнч Баттерворта имеет монотонный характер в ней.
3) Данный метод видимо предназначен для очень низких частот (килогерцы уже не прокатят), как я понял. Иначе работая с килогерцами к примеру, как я попытался это сделать, даже закрыв глаза на сомнения вышеуказанные, я просто зарылся в вычислениях и в конце концов не получил ничего путного разумеется.
0
12 / 12 / 0
Регистрация: 20.11.2013
Сообщений: 167
16.02.2016, 12:57  [ТС] 19
Alex_SP, чтобы не зарываться в расчетах посмотрите, пожалуйста, эту тему: Реализация ФНЧ Баттерворта
Я в ней привел просто свою реализацию расчета ФНЧ Баттерворта по мурзилке (http://www.dsplib.ru/content/f... terex.html).
Вот только никак не разберусь как реализовать формирование блоков с коэффициентами как в Матлабе сделано. По идее - их удобнее будет использовать при высоких порядках фильтров. Если разберетесь с этим - поделитесь опытом, пожалуйста
0
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
16.02.2016, 13:56 20
Alex_SP, DimKaKiber, метод матриц z-преобразования предполагает, что для получения коэффициентов любого цифрового фильтра (ФНЧ, ФВЧ, полосового, режекторного) необходимо знать только коэффициенты аналогового ФНЧ -прототипа.
0
IT_Exp
Эксперт
87844 / 49110 / 22898
Регистрация: 17.06.2006
Сообщений: 92,604
16.02.2016, 13:56
Помогаю со студенческими работами здесь

Вычисление биномиальных коэффициентов
Задание во вложении

Вычисление биномиальных коэффициентов
Напишите рекурсивную программу, которая по заданным N и M вычисляет С(m,n)=C(m,n-1)+C(m-1,n-1), где...

Вычисление коэффициентов на стипендию
Примечание: Коэффициент на стипендию (k) зависит от значения среднего балла за семестр. Если...

Вычисление биномиальных коэффициентов
Напишите рекурсивную программу, которая по заданным N и M вычисляет Входные данные Два целых...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
20
Ответ Создать тему
Опции темы

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin
Copyright ©2000 - 2022, CyberForum.ru