1159 / 825 / 400
Регистрация: 21.10.2012
Сообщений: 2,395
1

Реализация БИХ-фильтра

06.04.2016, 16:27. Показов 1058. Ответов 1
Метки нет (Все метки)

Здравствуйте. Задача - необходимо реализовать БИХ фильтр на языке С. Нашел тему, в которой уже приводится готовый код, но он чуть-чуть не подходит, т.к. для моей задачи можно использовать только целочисленные данные. Попытался изменить код, но результат фильтрации не верен. Не могли бы пояснить ошибку?


Исходный код:
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
#include "malloc.h"
 
void butter4(double *in, double *out, int length, const double NUM[], const double DEN[])
{
    double x, y0, y1 = 0, y2 = 0, y3 = 0;
    double b0= NUM[0];
    double b1= NUM[1];
    double b2= NUM[2];
    double b3= NUM[3];
 
    double a0= DEN[0];
    double a1= DEN[1];
    double a2= DEN[2];
    double a3= DEN[3];
 
    for(int i= 0; i < length; i++)
    {
        x = *in++;
        y0 = x * b0 + y1;
        y1 = b1 * x - a1 * y0 + y2;
        y2 = b2 * x - a2 * y0 + y3;
        y3 = b3 * x - a3 * y0;
        *out++ = y0 * a0;
    }
}
 
int main(int argc, char* argv[])
{
    int Length = 1000;
    double Num[] = {6.30129962118151e-06, 1.89038988635445e-05, 1.89038988635445e-05, 6.30129962118151e-06};
    double Den[] = {1, -2.92520452880893, 2.85318010843919, -0.927925169233285};
 
    double *signal = (double*)calloc( Length, sizeof(double) );
    signal[0] = 1.0;
 
    double *filtredSignal = (double*)malloc( sizeof(double) * Length );
 
    butter4(signal, filtredSignal, Length, Num, Den);
 
    FILE *f = fopen("ImpResponse.txt","w+");
    for(int i=0; i<Length; i++)
        fprintf(f, "%.15e%c", filtredSignal[i] , ' ');
    fclose(f);
 
    free(signal);
    free(filtredSignal);
    return 0;
}

А это уже переделанный мной:
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
#include <stdio.h>
#include <stdlib.h>
#include "malloc.h"
 
void butter3(short int *in, long *out, int length, const long NUM[], const long DEN[])
{
    long x, ex, y0, y1 = 0, y2 = 0, y3 = 0;
    long b0= NUM[0];
    long b1= NUM[1];
    long b2= NUM[2];
    long b3= NUM[3];
 
    long a0= DEN[0];
    long a1= DEN[1];
    long a2= DEN[2];
    long a3= DEN[3];
 
    for(int i= 0; i < length; i++)
    {
        x = *in++;
        y0 = x * b0 + y1;
        y1 = b1 * x - a1 * y0 + y2;
        y2 = b2 * x - a2 * y0 + y3;
        y3 = b3 * x - a3 * y0;
        *out++ = y0 * a0;
    }
}
 
int main()
{
    int Length = 1200;
    long Num[] = {4654133, 13605355, 13605355, 4654133};
    long Den[] = {100000000, -138244832, 110745415, -35981606};
 
    short int *signal = (short int*)malloc( sizeof(short int) * Length  );
    for(int i = 0; i < 100; i++)
    {
        signal[12*i] = 866;
        signal[1+12*i] = -500;
        signal[2+12*i] = 0;
        signal[3+12*i] = 500;
        signal[4+12*i] = -866;
        signal[5+12*i] = -2000;
        signal[6+12*i] = -866;
        signal[7+12*i] = 500;
        signal[8+12*i] = 0;
        signal[9+12*i] = -500;
        signal[10+12*i] = 866;
        signal[11+12*i] = 2000;
    }
 
 
    long *filtredSignal = (long*)malloc( sizeof(long) * Length );
 
    butter3(signal, filtredSignal, Length, Num, Den);
 
    FILE *lol;
    lol = fopen("E:\\Poklon\\SignalInt", "wb");
    for (int q = 0; q < Length; q++) {
        fwrite(&filtredSignal[q], sizeof(long), 1, lol);
    };
    fclose(lol);
 
    free(signal);
    free(filtredSignal);
    return 0;
}
P.S. Значения сигнала, числителя и знаменателя получал из Matlab-а. Проверку производил там же.

Добавлено через 3 часа 22 минуты
Переделал для целочисленных величин вот в такой вид, но опять что-то не так, где-то есть ошибка

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
/* Numerator coefficients */
#define B0 0
#define B1 1
#define B2 2
 
/* Denominator coefficients */
#define A0 3
#define A1 4
#define A2 5
 
void my( const signed int * coefficients, signed int * input, signed int * output, int length)
{
  long temp;
  static signed int x[3] = { 0, 0, 0 };  /* x(n), x(n-1), x(n-2) */
  static signed int y[3] = { 0, 0, 0 };  /* y(n), y(n-1), y(n-2) */
 
    for(int i = 0; i < length; i++)
    {
     x[0] = input[i]; /* Copy input to x[0] */
 
     temp =  ( (long) coefficients[B0] * x[0]) ;   /* B0 * x(n)     */
     temp += ( (long) coefficients[B1] * x[1]);    /* B1 * x(n-1) */
     temp += ( (long) coefficients[B2] * x[2]);    /* B2 * x(n-2)   */
     temp -= ( (long) coefficients[A1] * y[1]);    /* A1 * y(n-1) */
     temp -= ( (long) coefficients[A2] * y[2]);    /* A2 * y(n-2)   */
 
     /* Здесь происходит деление на А0, которое равно 32767, поэтому просто сдвиг вправо на 15 бит*/
     temp >>= 15;
 
     y[0] = (short int) ( temp );
 
     /* Сдвиг результатов */
 
     y[2] = y[1];   /* y(n-2) = y(n-1) */
     y[1] = y[0];   /* y(n-1) = y(n)   */
 
     x[2] = x[1];   /* x(n-2) = x(n-1) */
     x[1] = x[0];   /* x(n-1) = x(n)   */
 
     output[i] = (short int) temp;
    }
}
Добавлено через 21 минуту
Проверил код для сигнала с частотами 20 Гц и 60 Гц. После фильтрации (ФНЧ с частотой среза 40 Гц) вместо исходных составляющих я получаю 2 пика: один на частоте 10 Гц, другой на частоте 30 Гц (и он меньше по уровню, как будто бы подавился фильтром). Т.е. почему-то исходные частоты сигнала уменьшились в 2 раза. Но помимо этого я получил еще составляющие симметричные относительно 60 Гц (90 Гц - подавленный и 110 Гц). Может кто-нибудь подсказать причину?

Добавлено через 10 минут
Итог: ошибка была в том, что я неправильно считывал данные из файла в Matlab Последний код рабочий
0
Programming
Эксперт
94731 / 64177 / 26122
Регистрация: 12.04.2006
Сообщений: 116,782
06.04.2016, 16:27
Ответы с готовыми решениями:

Построение АЧХ цифрового БИХ фильтра и его реализация
Добрый день! Заранее извиняюсь, если задаю откровенно глупые вопросы, но не могу никак въехать как...

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

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

Реализация CIC-фильтра
Собственно, схема: Интегратор состоит из накапливающих сумматоров. А сумматорам надо...

1
0 / 0 / 0
Регистрация: 31.03.2016
Сообщений: 29
22.06.2022, 14:34 2
Добрый день.
Подскажите пожалуйста, правильно ли я понимаю, что в исходном коде с double точностью рассмотрен пример БИХ фильтра третьего порядка? Вроде как должны быть ещё коэффициенты b4 и a4 (для фильтра 4 порядка)?
0
IT_Exp
Эксперт
87844 / 49110 / 22898
Регистрация: 17.06.2006
Сообщений: 92,604
22.06.2022, 14:34
Помогаю со студенческими работами здесь

Реализация КИХ фильтра(теория)
http://www.dsplib.ru/content/filters/butterex/butterex.html вот тут про Батерворта. Про Ких ...

Реализация полосового фильтра Баттерворта 3 порядка
Использование z-преобразования Ннч(S)--&gt; Hпф(z) с получением соответствующей матрицы, позволяет...

Kit-модуль MSP-EXP430G2. Реализация цифрового фильтра
Здравствуйте знатоки. Кто-нибудь из вас пытался реализовать на базе контроллера MSP-EXP430G2...

Реализация преобразования Гильберта с использованием КИХ-фильтра первого порядка
У меня такой вопрос, скажите можно ли реализовать преобразование Гильберта с помощью КИХ-фильтра...


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

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

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