Реализация БИХ-фильтра
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;
}
} |
|
Проверил код для сигнала с частотами 20 Гц и 60 Гц. После фильтрации (ФНЧ с частотой среза 40 Гц) вместо исходных составляющих я получаю 2 пика: один на частоте 10 Гц, другой на частоте 30 Гц (и он меньше по уровню, как будто бы подавился фильтром). Т.е. почему-то исходные частоты сигнала уменьшились в 2 раза. Но помимо этого я получил еще составляющие симметричные относительно 60 Гц (90 Гц - подавленный и 110 Гц). Может кто-нибудь подсказать причину?
Добавлено через 10 минут
Итог: ошибка была в том, что я неправильно считывал данные из файла в Matlab  Последний код рабочий
0
|
Последний код рабочий
(
