Фильтрация сигнала (численные методы, C++)

[colding]

Здравствуйте.
Выполняю фильтрацию сигнала мгновенной мощности (т.е. имеются датчики тока и напряжения, с них синхронно поступают данные с определенной частотой). Эта фильтрация позволяет выделить постоянную поставляющую мощности, а именно активную составляющую мощности.
На рисунке представлена схема в matlab, реализующая данную фильтрацию ( ссылка на яндекс диск https://yadi.sk/i/IntrXGamriNhG , не получилось залить картинку сюда).

Задача состоит в том, чтобы реализовать данную фильтрацию на С++. Для интегрирования выбран метод Эйлера (одношаговый).
Попробовал реализовать, но что-то не так сделал. Подскажите, что не так или вообще может другой способ реализации посоветуете.

Используется синтаксис C++

// double Current[] - массив значений мгновенных токов
// double Voltage[] - массив значений мгновенных напряжений
// double Step - шаг опроса
// int N - количество значений в массиве токов (напряжения)

double Filter (double Current[], double Voltage[], double Step, int N) {
        double y1Past=0; // предыдущее значение y1
        double yPast=0;  // предыдущее значение y
        double xPast=0;  // предыдущее значение х
        double y1=0;     // текущее значение y1
        double y=0;      // текущее значение y
        double x=0;      // текущее значение х
        double ksi=0.707;
        double T=1/(5*M_PI);
        for (int i = 0; i < N; i++) {
                x= (Current[i] * Voltage[i] - y1Past - 2*ksi*T*yPast)/(T*T);
                y= yPast + (x-xPast)*Step;
                xPast = x;
                y1 = y1Past + (y-yPast) * Step;
                yPast = y1 ;
        }
                return y1;
}
}

[svv]

(Неверная схема удалена.)

Распилите Ваш вопрос пополам: запишите в $\TeX$ математическую формулу для текущего значения $y_1$, и получатся два отдельных более простых вопроса:
$\bullet$ Правильная ли формула?
$\bullet$ При условии, что формула верна, правильно ли Вы реализовали её в программе?

-- Вс май 15, 2016 12:33:22 --

Судя по схеме, использовать надо предыдущее значение $y$, а Вы оператором yPast=y1; запоминаете значение $y_1$.

[colding]

Извините, скинул не тот скрин.
Вот нужный скрин

svv в сообщении #1123667 писал(а):

Судя по схеме, использовать надо предыдущее значение $y$, а Вы оператором yPast=y1; запоминаете значение $y_1$.

Это тоже я заметил.

1. после преобразования передаточная функция примет вид W(p)=1/(T^2\cdot p^2+2\cdot\xi\cdot T\cdot p+1)
$W(p)=y1/u$
получаем $T^2 \cdot y1''+2\cdot\xi\cdot T\cdot y1'+y1=u$
относительно старшей производной $y1''=(u-2\cdot\xi\cdot T\cdot y1'-y1)/T^2$

Код выглядит так

Используется синтаксис C++

// double Current[] - массив значений мгновенных токов
// double Voltage[] - массив значений мгновенных напряжений
// double Step - шаг опроса
// int N - количество значений в массиве токов (напряжения)

double Filter (double Current[], double Voltage[], double Step, int N) {
        double y1Past=0; // предыдущее значение y1
        double yPast=0;  // предыдущее значение y
        double xPast=0;  // предыдущее значение х
        double y1=0;     // текущее значение y1
        double y=0;      // текущее значение y
        double x=0;      // текущее значение х
        double ksi=0.707;
        double T=1/(5*M_PI);
        for (int i = 0; i < N; i++) {
                x= (Current[i] * Voltage[i] - y1Past - 2*ksi*T*yPast)/(T*T);
                y= yPast + (x-xPast)*Step;
                xPast = x;
                y1 = y1Past + (y-yPast) * Step;
                yPast=y;
                y1Past = y1 ;
        }
                return y1;
}

Реализация все равно не корректна :(

[Joker_vD]

Ну тогда ручками считаете шага три или даже четыре, потом вставляете внутрь цикла отладочную печать, запускаете, читаете вывод, сравниваете с посчитанным вручную, находите, где начинаются расхождения и исправляете.

[svv]

У Вас неправильно организовано интегрирование. Не нужно вычитать xPast и yPast.
Думаю, это не единственная ошибка.

[colding]

svv, Joker_vD спасибо за помощь. Результат программы на С++ и моделирования в матлаб на разных участках (от 0 до 1 секунды) очень близки. Сказывается маленький шаг дискретизации.
Еще раз спасибо, без Вас ничего бы не получилось.

[svv]

colding
Спасибо за благодарность , но всё-таки имейте в виду, что вот такое:
$s_i=s_{i-1}+(y_i-y_{i-1})h$
— это не интегрирование $y$ ни в каком приближении. (обозначения мои)

[colding]

svv в сообщении #1123916 писал(а):

colding
Спасибо за благодарность , но всё-таки имейте в виду, что вот такое:
$s_i=s_{i-1}+(y_i-y_{i-1})h$
— это не интегрирование $y$ ни в каком приближении. (обозначения мои)

Правда Ваша

Может кому-нибудь пригодится (но думаю, что этот материал мало ценен):

Используется синтаксис C++

// double Current[] - массив значений мгновенных токов
// double Voltage[] - массив значений мгновенных напряжений
// double Step - шаг опроса
// int N - количество значений в массиве токов (напряжения)

double Filter (double Current[], double Voltage[], double Step, int N) {
        double y1Past=0; // предыдущее значение y1
        double yPast=0;  // предыдущее значение y
        double xPast=0;  // предыдущее значение x
        double y1;       // текущее значение y1
        double y;        // текущее значение y
        double x;        // текущее значение x
        double ksi=0.707;
        double T=1/(5*M_PI);
        for (int i = 0; i < N; i++) {
                x = (Current[i] * Voltage[i] - y1Past - 2*ksi*T*yPast)/(T*T);
                y = yPast + (x)*Step;
                y1 = y1Past + (y) * Step;
                yPast= y;
                y1Past = y1 ;

        }
                return y1;
}

Результат отрисовки:


Для сравнения, приложу результаты матлаба

[svv]

Отлично!