Minimum working example below. The Matlab filter function filters the input data x
using a rational transfer function defined by the numerator and denominator coefficients b
and a
and initial conditions z
. My C++ version of the code is way slower than the Matlab one. My signal x
is a huge valarray, a
and b
are three elements long each and the initial conditions z
is just three zeros.
#include <iostream>
#include<ctime>
#include <cstdlib>
using std::rand;
#include <valarray>
using std::valarray;
#include <ctime>
#define D_SCL_SECURE_NO_WARNINGS 1
struct FilterResults
{
valarray<double> filteredSignal;
valarray<double> finalConditions;
};
FilterResults Filter(const valarray<double> &b,
const valarray<double> &a, const valarray<double> &x, const valarray<double> &z)
{
valarray<double> y(x.size());
// Pad z
valarray<double> zOut(a.size());
std::copy(std::begin(z), std::end(z), std::begin(zOut));
double Xm, Ym;
size_t i, n = a.size();
for (size_t m = 0; m < y.size(); m++)
{
Xm = x[m];
y[m] = b[0] * Xm + zOut[0];
Ym = y[m];
for (i = 1; i < n; i++)
{
zOut[i - 1] = b[i] * Xm + zOut[i] - a[i] * Ym;
}
}
valarray<double> s = zOut[std::slice(0, zOut.size() - 1, 1)];
FilterResults r;
r.filteredSignal = std::move(y);
r.finalConditions = std::move(s);
return r;
}
int main() {
std::srand(std::time(0));
valarray<double> b{ { (double)rand(), (double)rand(), (double)rand() } };
valarray<double> a{ { (double)rand(), (double)rand(), (double)rand() } };
valarray<double> z{ { 0, 0, 0 } };
valarray<double> x( 500000 );
for (size_t i = 0; i < 500000; i++)
{
x[i] = (double)rand();
}
clock_t begin = clock();
auto r = Filter(b, a, x, z);
clock_t end = clock();
double elapsed_secs = double(end - begin) / CLOCKS_PER_SEC;
std::cout << "elapsed_secs: " << elapsed_secs << std::endl;
system("pause"); // C++ is 0.184 seconds on my computer, Matlab is 0.006060.
}
Any idea on how to speed it up?