Форум программистов, компьютерный форум, киберфорум
Matlab
Войти
Регистрация
Восстановить пароль
 
Рейтинг 4.92/26: Рейтинг темы: голосов - 26, средняя оценка - 4.92
9 / 9 / 1
Регистрация: 25.10.2009
Сообщений: 152
1

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

16.05.2012, 22:37. Просмотров 4658. Ответов 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)
Programming
Эксперт
94731 / 64177 / 26122
Регистрация: 12.04.2006
Сообщений: 116,782
16.05.2012, 22:37
Ответы с готовыми решениями:

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

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

Можно ли в функции высшего порядка filter получить номер элемента входящего списка
Можно ли в функции высшего порядка filter получить номер элемента входящего списка lst...

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

7
2008 / 1280 / 60
Регистрация: 05.06.2010
Сообщений: 2,213
18.05.2012, 10:40 2
Лучший ответ Сообщение было отмечено 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  [ТС] 3
Спасибо огромное!) теперь кажется стало понятно))
0
3 / 3 / 1
Регистрация: 05.05.2014
Сообщений: 75
01.04.2015, 00:44 4
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
2008 / 1280 / 60
Регистрация: 05.06.2010
Сообщений: 2,213
01.04.2015, 10:01 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
...
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
3 / 3 / 1
Регистрация: 05.05.2014
Сообщений: 75
02.04.2015, 00:58 6
vital792,
Так получается если 9 коэффициентов, то так нужно переделать?
C#
1
2
int order = 2;
double[] t = new double[order + 7];
Или как? Зачем вообще эта вторая строчка? если можно сразу обозначить ордер равно 9?
0
2008 / 1280 / 60
Регистрация: 05.06.2010
Сообщений: 2,213
02.04.2015, 09:15 7
Цитата Сообщение от 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
3 / 3 / 1
Регистрация: 05.05.2014
Сообщений: 75
05.04.2015, 17:00 8
vital792, что-то я ничего уже не понимаю... Не смог разобраться со своей проблемой, хотя ордер задаю правильно, все равно ошибка. Потом попробовал Ваш пример, тоже ошибка, но уже другая.
Поможете разобраться мне? Прикрепляю скрины ошибок и кода.
0
Миниатюры
Смысл функции filter. Перевести код функции filter в C++   Смысл функции filter. Перевести код функции filter в C++  
IT_Exp
Эксперт
87844 / 49110 / 22898
Регистрация: 17.06.2006
Сообщений: 92,604
05.04.2015, 17:00

Заказываю контрольные, курсовые, дипломные и любые другие студенческие работы здесь.

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

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

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

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


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

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

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.