Форум программистов, компьютерный форум, киберфорум
Matlab
Форум программистов и сисадминов Киберфорум > Научный форум > Математические программы > Matlab
Войти
Регистрация
Восстановить пароль
Блоги Сообщество Поиск Заказать работу  
 
Рейтинг 4.57/35: Рейтинг темы: голосов - 35, средняя оценка - 4.57
9 / 9 / 1
Регистрация: 25.10.2009
Сообщений: 152

Смысл функции filter. Перевести код функции filter в C++

16.05.2012, 22:37. Показов 7406. Ответов 7
Метки нет (Все метки)
Студворк — интернет-сервис помощи студентам
Всем привет! у меня такая проблема. Нужно перевести код функции filter в C++.
для функции
Matlab M
1

y = filter (B, A, x)
я нашла код в интернете:
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

void filter(int ord, float *a, float *b, int np, float *x, float *y)
{
        int i,j;
    y[0]=b[0]*x[0];
    for (i=1;i<ord+1;i++)
    {
        y[i]=0.0;
        for (j=0;j<i+1;j++)
            y[i]=y[i]+b[j]*x[i-j];
        for (j=0;j<i;j++)
            y[i]=y[i]-a[j+1]*y[i-j-1];
    }
/* end of initial part */
for (i=ord+1;i<np+1;i++)
{
    y[i]=0.0;
        for (j=0;j<ord+1;j++)
        y[i]=y[i]+b[j]*x[i-j];
        for (j=0;j<ord;j++)
        y[i]=y[i]-a[j+1]*y[i-j-1];
}
} /* end of filter */
вроде с ним проблем нет. Не понятно как реализовать такую функцию:
Matlab M
1

[y, Sf] = filter (B, A, x, Si)
в частности не понятно как меняется аргумент Si. единственное что нашла, так это то что Si - начальное значение, а Sf - конечное. подскажите, пожалуйста, кто знает, что происходит с этим аргументом???
Заранее спасибо!!!
0
Лучшие ответы (1)
Сообщение: #3047907
IT_Exp
Эксперт
34794 / 4073 / 2104
Регистрация: 17.06.2006
Сообщений: 32,602
Блог
 16.05.2012, 22:37
Ответы с готовыми решениями:

Фильтрация сигнала (создать аналог функции filter на С#)
Привет! Помогите, пожалуйста, создать функцию фильтра, подобную функции MatLabа = filter(b, a, x, Zi). Отфильтровать необходимо...

Использование собственной функции в BindingSource Filter
Добрый день ! Стоит задача фильтровать данные в BindingSource с использованием своей функции. Но в лоб это не получается. Подскажите...

Можно ли в функции высшего порядка filter получить номер элемента входящего списка
Можно ли в функции высшего порядка filter получить номер элемента входящего списка lst удовлетворяющему критерию lambda x: x == '#'? lst...

7
2014 / 1286 / 61
Регистрация: 05.06.2010
Сообщений: 2,213
18.05.2012, 10:40
Лучший ответ Сообщение было отмечено R2D2 как решение

Решение

Цитата Сообщение от Flamе Посмотреть сообщение
вроде с ним проблем нет.
нет? покажите как работает. Странный какой то код - используется рекурсивный фильтр 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
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

// filter_state.cpp
//
 
#include <cmath>
#include <fstream>
 
/*
void filter_(int ord, float *a, float *b, int np, float *x, float *y)
{
    int i,j;
    y[0] = b[0]*x[0];
    for (i=1; i<ord+1; i++)
    {
        y[i] = 0.0;
        for (j=0; j<i+1; j++)
            y[i] = y[i] + b[j]*x[i-j];
        for (j=0; j<i; j++)
            y[i] = y[i] - a[j+1]*y[i-j-1];
    }
    // end of initial part 
    for (i=ord+1;i<np+1;i++)
    {
        y[i] = 0.0;
        for (j=0; j<ord+1; j++)
            y[i] = y[i] + b[j]*x[i-j];
        for (j=0; j<ord; j++)
            y[i] = y[i] - a[j+1]*y[i-j-1];
    }
} // end of filter 
*/
 
void filter(int ord, const double a[], const double b[], int length, double *x, double *y)
{
    double *t = new double[ord+1];
    memset(t, 0, (ord+1)*sizeof(double));
 
    for(int i=0; i<length; i++)
    {
        double xi = x[i];
        t[0] = b[0]*xi + t[1];
        double t0 = t[0];
        for(int j=1; j<ord; j++)
        {
            t[j] = b[j]*xi - a[j]*t0 + t[j+1];
        }
        t[ord] = b[ord]*xi - a[ord]*t0;
        y[i] = a[0]*t0;
    }
 
    delete [] t;
}
 
enum freqs
{
    f1 = 500, f2 = 1000, f3 = 2000
};
 
int main()
{
    const int Fs = 10000;
    const double pi = 3.14159;
    unsigned length = 1000;
    double temp;
    double *x = new double[length];
 
    for(unsigned i=0; i<length; i++)
    {
        temp = sin(2*pi*freqs::f1*i/Fs);
        temp += sin(2*pi*freqs::f2*i/Fs);
        temp += sin(2*pi*freqs::f3*i/Fs);
        x[i] = temp;
    }
 
    // Butterworth filter
    double b[] = {1.68358140723680e-06, 1.68358140723680e-05, 7.57611633256558e-05,
        0.000202029768868415, 0.000353552095519727, 0.000424262514623672, 0.000353552095519727,
        0.000202029768868415, 7.57611633256558e-05, 1.68358140723680e-05, 1.68358140723680e-06};
 
    double a[] = {1, -5.98758962981667, 16.6721933230027, -28.2587879002006,
        32.1597564876946, -25.6017495970534, 14.4056874262078, -5.64707434413249,
        1.47372793697391, -0.230919345862029, 0.0164796305471309};
 
    int order = 10;
    double *y = new double[length];
 
    //filter_(order, a, b, length, x, y);
    filter(order, a, b, length, x, y);
 
    std::ofstream f("filtred.txt");
    for(int i=0; i<length; i++)
        f << y[i] << ' ';
    f.close();
 
    delete [] x;
    delete [] y;
 
    return 0;
}
И код матлаб для проверки и сравнения:
Matlab M
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

clear all; close all; clc;
 
Fs = 10000;
t = (0:1000)/Fs;
 
x = sin(2*pi*500*t) + sin(2*pi*1000*t) + sin(2*pi*2000*t);
figure(1)
plot(x);
 
n = 10;
[b, a] = butter(n, .2, 'low');
 
% так не надо делать
y1 = filter(b, a, x(1:500));
y2 = filter(b, a, x(501:end));
y = [y1 y2];
figure(2)
plot(y);
 
% а так надо
[y1, zf] = filter(b, a, x(1:500));
y2 = filter(b, a, x(501:end), zf);
y = [y1 y2];
figure(3)
plot(y);
 
figure(4)
plot(abs(fft(y))); axis tight;
 
ss = load('D:\test\_\filter_state\filtred.txt'); % загружаем сишный результат
figure(5)
plot(abs(fft(ss))); axis tight;
Цитата Сообщение от Flamе Посмотреть сообщение
Si - начальное значение, а Sf - конечное. подскажите, пожалуйста, кто знает, что происходит с этим аргументом???
В своей функции я выделяю память для хранения текущего состояния фильтра:
C++
1

double *t = new double[ord+1];
и обнуляю эти значения:
C++
1

    memset(t, 0, (ord+1)*sizeof(double));
В случае фильтрации сигнала по частям, эти операции нужно выполнять вне функции, передавая указатель на состояние фильтра как аргумент функции. Перед первым вызовом эти значения должны быть нулями, при выходе из функции будет сохраняться состояние фильтра. Доработайте сами - это не сложно, пару строчек изменить.
3
9 / 9 / 1
Регистрация: 25.10.2009
Сообщений: 152
18.05.2012, 19:54  [ТС]
Спасибо огромное!) теперь кажется стало понятно))
0
 Аватар для starjke
3 / 3 / 1
Регистрация: 05.05.2014
Сообщений: 75
01.04.2015, 00:44
Flamе, а как понять, какой нужен order? Альфа ритм выдает ошибку..
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

            //50hz
            int order = 2;
            double[] t = new double[order + 1];
            double[] b = {0.98300511284837244, 0.24971833246806957, 0.98300511284837244};
            double[] a = {1, 0.24971833246806957, 0.96601022569674488};
            double[] y = new double[num];
            filter(order, a, b, num, L1, y, t); 
            chart4.Series[0].Points.DataBindXY(Lx,y);
            chart4.ChartAreas[0].AxisX.Minimum = 0;
            chart4.ChartAreas[0].AxisX.Maximum = 925;
            //lowpass
            double[] t2 = new double[order + 1];
            double[] b2 = { 0.059560596234149672,0.23824238493659869,0.35736357740489805,0.23824238493659869,0.059560596234149672 };
            double[] a2 = {  1,-0.52804595798372822,0.57428885374252248,-0.11645534850163169,0.023181992489232316 };
            double[] y2 = new double[num];
            filter(order, a2, b2, num, L1, y2, t2);
            chart5.Series[0].Points.DataBindXY(Lx, y2);
            chart5.ChartAreas[0].AxisX.Minimum = 0;
            chart5.ChartAreas[0].AxisX.Maximum = 925;
            //alpha-rhythm
            double[] t3 = new double[order + 1];
            double[] b3 = {0.000042022232781020762, 0, -0.00016808893112408305, 0, 0.00025213339668612457, 0, -0.00016808893112408305, 0, 0.000042022232781020762};
            double[] a3 = {1, -7.1066864371216631, 22.515295948777023, -41.498099095062784, 48.645183603782698, -37.132329737715914, 18.027451528765006, -5.0918788345024755, 0.64122389961493442};
            double[] y3 = new double[num];
            filter(order, a3, b3, num, L1, y3, t3);
            chart7.Series[0].Points.DataBindXY(Lx, y3);
            chart7.ChartAreas[0].AxisX.Minimum = 0;
            chart7.ChartAreas[0].AxisX.Maximum = 925;
         }
 private void filter(int ord, double[] a, double[] b, int length, double[] x, double[] y, double[] t)
    {
 
        for (int i = 0; i < length; i++)
        {
            double xi = x[i];
            t[0] = b[0] * xi + t[0];
            double t0 = t[0];
            for (int j = 0; j < ord-1; j++)
            {
                t[j] = b[j+1] * xi - a[j+1] * t0 + t[j + 1];
            }
            t[ord-1] = b[ord] * xi - a[ord] * t0;
            y[i] = a[0] * t0;
        }
        t = null;
    }
Добавлено через 12 минут
vital792, можете просмотреть код? В чем ошибка? Почему не работают банд-пасс фильтры для выделения ритмов? Фильтры 50гц и фнч работают.
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

     //50hz
            int order = 2;
            double[] t = new double[order + 1];
            double[] b = {0.98300511284837244, 0.24971833246806957, 0.98300511284837244};
            double[] a = {1, 0.24971833246806957, 0.96601022569674488};
            double[] y = new double[num];
            filter(order, a, b, num, L1, y, t); 
            chart4.Series[0].Points.DataBindXY(Lx,y);
            chart4.ChartAreas[0].AxisX.Minimum = 0;
            chart4.ChartAreas[0].AxisX.Maximum = 925;
            //lowpass
            double[] t2 = new double[order + 1];
            double[] b2 = { 0.059560596234149672,0.23824238493659869,0.35736357740489805,0.23824238493659869,0.059560596234149672 };
            double[] a2 = {  1,-0.52804595798372822,0.57428885374252248,-0.11645534850163169,0.023181992489232316 };
            double[] y2 = new double[num];
            filter(order, a2, b2, num, L1, y2, t2);
            chart5.Series[0].Points.DataBindXY(Lx, y2);
            chart5.ChartAreas[0].AxisX.Minimum = 0;
            chart5.ChartAreas[0].AxisX.Maximum = 925;
            //alpha-rhythm
            double[] t3 = new double[order + 1];
            double[] b3 = {0.000042022232781020762, 0, -0.00016808893112408305, 0, 0.00025213339668612457, 0, -0.00016808893112408305, 0, 0.000042022232781020762};
            double[] a3 = {1, -7.1066864371216631, 22.515295948777023, -41.498099095062784, 48.645183603782698, -37.132329737715914, 18.027451528765006, -5.0918788345024755, 0.64122389961493442};
            double[] y3 = new double[num];
            filter(order, a3, b3, num, L1, y3, t3);
            chart7.Series[0].Points.DataBindXY(Lx, y3);
            //chart7.ChartAreas[0].AxisX.Minimum = 0;
            //chart7.ChartAreas[0].AxisX.Maximum = 925;
            //beta-rhythm
            double[] t4 = new double[order + 1];
            double[] b4 = { 0.015014740623995418, 0 -0.060058962495981673, 0, 0.090088443743972513, 0, -0.060058962495981673, 0, 0.015014740623995418};
            double[] a4 = {1, -3.8483506604531263, 7.4219375225148241, -9.2412355487591089, 8.1272387364941743, -5.1014790793630498, 2.2315198084593453, -0.62514756293261808, 0.0907551178832674};
            double[] y4 = new double[num];
            filter(order, a4, b4, num, L1, y4, t4);
            chart8.Series[0].Points.DataBindXY(Lx, y4);
            chart8.ChartAreas[0].AxisX.Minimum = 0;
            chart8.ChartAreas[0].AxisX.Maximum = 925;
            //delta-rhythm
            double[] t5 = new double[order + 1];
            double[] b5 = { 0.0000029142055881038877, 0, -0.000011656822352415551, 0, 0.000017485233528623325, 0, -0.000011656822352415551, 0, 0.0000029142055881038877};
            double[] a5 = {1, -7.7714250554471498, 26.432110976797333, -51.389986320680549, 62.46829391139503, -48.615650506783084, 23.655312859527022, -6.5796087124733509, 0.80095284767278929};
            double[] y5 = new double[num];
            filter(order, a5, b5, num, L1, y5, t5);
            chart6.Series[0].Points.DataBindXY(Lx, y5);
            chart6.ChartAreas[0].AxisX.Minimum = 0;
            chart6.ChartAreas[0].AxisX.Maximum = 925;
            //teta-rhythm
            double[] t6 = new double[order + 1];
            double[] b6 = {0.0000012194060670318125, 0, -0.0000048776242681272498, 0, 0.0000073164364021908743, 0, -0.0000048776242681272498, 0, 0.0000012194060670318125};
            double[] a6 = {1, -7.7144455673812846, 26.143432615848749, -50.832116666968687, 62.020956046008855, -48.625039856626231, 23.92252106319896, -6.7526345681114215, 0.83732746738682529};
            double[] y6 = new double[num];
            filter(order, a6, b6, num, L1, y6, t6);
            chart9.Series[0].Points.DataBindXY(Lx, y6);
            chart9.ChartAreas[0].AxisX.Minimum = 0;
            chart9.ChartAreas[0].AxisX.Maximum = 925;
         }
 private void filter(int ord, double[] a, double[] b, int length, double[] x, double[] y, double[] t)
    {
 
        for (int i = 0; i < length; i++)
        {
            double xi = x[i];
            t[0] = b[0] * xi + t[0];
            double t0 = t[0];
            for (int j = 0; j < ord-1; j++)
            {
                t[j] = b[j+1] * xi - a[j+1] * t0 + t[j + 1];
            }
            t[ord-1] = b[ord] * xi - a[ord] * t0;
            y[i] = a[0] * t0;
        }
        t = null;
    }
0
2014 / 1286 / 61
Регистрация: 05.06.2010
Сообщений: 2,213
01.04.2015, 10:01
порядок фильтра определяется по максимальной длине массива коэффициентов
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

...
const double twopi = 3.14159;
int order, num = 1000;
//alpha-rhythm
double[] b3 = { 0.000042022232781020762, 0, -0.00016808893112408305, 0, 0.00025213339668612457, 0, -0.00016808893112408305, 0, 0.000042022232781020762 };
double[] a3 = { 1, -7.1066864371216631, 22.515295948777023, -41.498099095062784, 48.645183603782698, -37.132329737715914, 18.027451528765006, -5.0918788345024755, 0.64122389961493442 };
double[] y3 = new double[num];
if (b3.Length >= a3.Length) order = b3.Length - 1;
else order = a3.Length - 1;
double[] t3 = new double[order + 1];
 
double [] L1 = new double[num + 1]; // test data
for (int i = 0; i < num + 1; i++)
{
    L1[i] = Math.Sin(twopi * i * .05) + Math.Sin(twopi * i * .1) + Math.Sin(twopi * i * .3); // sin before + inside + after band
}
for (int i = 0; i < order + 1; i++)
{
    t3[i] = 0;
}
 
// t3 must be zeros!!!
filter(order, a3, b3, num, L1, y3, t3);
...
1
 Аватар для starjke
3 / 3 / 1
Регистрация: 05.05.2014
Сообщений: 75
02.04.2015, 00:58
vital792,
Так получается если 9 коэффициентов, то так нужно переделать?
C#
1
2

int order = 2;
double[] t = new double[order + 7];
Или как? Зачем вообще эта вторая строчка? если можно сразу обозначить ордер равно 9?
0
2014 / 1286 / 61
Регистрация: 05.06.2010
Сообщений: 2,213
02.04.2015, 09:15
Цитата Сообщение от starjke Посмотреть сообщение
если можно сразу обозначить ордер равно 9?
если вы заранее знаете длину фильтра, то да, вы сразу знаете и его порядок. Но в вашем коде выше, используются несколько фильтров с разной длиной и соответственно порядком. Порядок (для БИХ) выбирается как максимальная длина из числителя (b) и знаменателя(a) минус 1. В моем фрагменте кода порядок фильтра так и вычисляется.
Цитата Сообщение от starjke Посмотреть сообщение
Зачем вообще эта вторая строчка?
Массив t содержит состояние фильтра. При первом вызове функции filter он должен быть заполнен нулями. При выходе из функции в нем будет содержаться текущее состояние фильтра. Это удобно при обработке длинных сигналов (или непрерывного сигнала поступающего с внешнего устройства). Вы не фильтруете сигнал целиком, а обрабатываете его по частям, применяя функцию filter для каждой выборки:
t = 0 // начальное состояние фильтра
for ...
x = signal(idx : idx + length); // выборка сигнала
y = filter(x, ..., t); // фильтрация. На выходе сохраняется состояние фильтра.
...
idx = idx + length;
end
1
 Аватар для starjke
3 / 3 / 1
Регистрация: 05.05.2014
Сообщений: 75
05.04.2015, 17:00
vital792, что-то я ничего уже не понимаю... Не смог разобраться со своей проблемой, хотя ордер задаю правильно, все равно ошибка. Потом попробовал Ваш пример, тоже ошибка, но уже другая.
Поможете разобраться мне? Прикрепляю скрины ошибок и кода.
Миниатюры
Смысл функции filter. Перевести код функции filter в C++   Смысл функции filter. Перевести код функции filter в C++  
0
Надоела реклама? Зарегистрируйтесь и она исчезнет полностью.
BasicMan
Эксперт
29316 / 5623 / 2384
Регистрация: 17.02.2009
Сообщений: 30,364
Блог
 05.04.2015, 17:00
Помогаю со студенческими работами здесь

Напишите функцию String->String, удаляющую из строки все гласные буквы используя функции filter и elem
Помогите, пожалуйста решить задачу. Напишите функцию String-&gt;String, удаляющую из строки все гласные буквы используя функции filter и...

Перевести код из функции в шаблоны функции
Короче, у меня есть прога написанная с помощью функций, теперь мне нужно написать её с помощью шаблонов функции... #include...

Filter
Почему нельзя сделать так: Private Sub Combo33_Change() Me.Filter = &quot;Name = Combo33.Value&quot; End SubКак правильно?

FS Filter
Нужно написать драйвер-фильтр файловой системы. Есть юзермодное приложение, в котором задаются шаблоны файлов, которые нужно скрывать...

Filter
Скажите пожалуйста как решить проблему с пробелами при фильтрации не работает и в combobox и в edit коды фильтрации ниже для combobox -...


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

Или воспользуйтесь поиском по форуму:
8
Ответ Создать тему
Новые блоги и статьи
SDL3 для Web (WebAssembly): Обработчик клика мыши в браузере ПК и касания экрана в браузере на мобильном устройстве
8Observer8 02.02.2026
Содержание блога Для начала пошагово создадим рабочий пример для подготовки к экспериментам в браузере ПК и в браузере мобильного устройства. Потом напишем обработчик клика мыши и обработчик. . .
Философия технологии
iceja 01.02.2026
На мой взгляд у человека в технических проектах остается роль генерального директора. Все остальное нейронки делают уже лучше человека. Они не могут нести предпринимательские риски, не могут. . .
SDL3 для Web (WebAssembly): Вывод текста со шрифтом TTF с помощью SDL3_ttf
8Observer8 01.02.2026
Содержание блога В этой пошаговой инструкции создадим с нуля веб-приложение, которое выводит текст в окне браузера. Запустим на Android на локальном сервере. Загрузим Release на бесплатный. . .
SDL3 для Web (WebAssembly): Сборка C/C++ проекта из консоли
8Observer8 30.01.2026
Содержание блога Если вы откроете примеры для начинающих на официальном репозитории SDL3 в папке: examples, то вы увидите, что все примеры используют следующие четыре обязательные функции, а. . .
SDL3 для Web (WebAssembly): Установка Emscripten SDK (emsdk) и CMake для сборки C и C++ приложений в Wasm
8Observer8 30.01.2026
Содержание блога Для того чтобы скачать Emscripten SDK (emsdk) необходимо сначало скачать и уставить Git: Install for Windows. Следуйте стандартной процедуре установки Git через установщик. . . .
SDL3 для Android: Подключение Box2D v3, физика и отрисовка коллайдеров
8Observer8 29.01.2026
Содержание блога Box2D - это библиотека для 2D физики для анимаций и игр. С её помощью можно определять были ли коллизии между конкретными объектами. Версия v3 была полностью переписана на Си, в. . .
Инструменты COM: Сохранение данный из VARIANT в файл и загрузка из файла в VARIANT
bedvit 28.01.2026
Сохранение базовых типов COM и массивов (одномерных или двухмерных) любой вложенности (деревья) в файл, с возможностью выбора алгоритмов сжатия и шифрования. Часть библиотеки BedvitCOM Использованы. . .
SDL3 для Android: Загрузка PNG с альфа-каналом с помощью SDL_LoadPNG (без SDL3_image)
8Observer8 28.01.2026
Содержание блога SDL3 имеет собственные средства для загрузки и отображения PNG-файлов с альфа-каналом и базовой работы с ними. В этой инструкции используется функция SDL_LoadPNG(), которая. . .
КиберФорум - форум программистов, компьютерный форум, программирование
Поддержать - Реклама - Условия использования - Обратная связь
Powered by vBulletin
Copyright ©2000 - 2026, CyberForum.ru