4
\$\begingroup\$

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?

\$\endgroup\$
  • \$\begingroup\$ Welcome to code review I hope you get some good answers. \$\endgroup\$ – pacmaninbw Sep 19 '16 at 12:42
  • \$\begingroup\$ Question, what is the order of your filter? It is much better, because of noise effects, to change your higher order filter to a number of cascaded filters with a lower order. In my company, we do (sub)nano precision control, we implement all higher order filters as cascaded filters with either order 2 or 1. \$\endgroup\$ – WG- Sep 19 '16 at 13:59
  • \$\begingroup\$ Would you like an optimized C implementation of a Filter function? \$\endgroup\$ – Royi Sep 1 '17 at 17:44
6
\$\begingroup\$

Few highlights:

std::copy(std::begin(z), std::end(z), std::begin(zOut));

To copy input values to overwrite them is just a waste of time. You will overwrite them then you do not need to copy old values, use z instead of zOut inside loops.

y.size() is possibly evaluated multiple times. It's a minor performance issue (and it may be inlined) but it depends on your actual implementation.

I don't know about valarray<T> implementation but if this is performance critical you may want to use a plain array.

I hope FilterResults constructor doesn't allocate anything (what valarray<double> default constructor does?) That said, this is C++, filter needs a state then you should design this using a Filter class with its private instance fields. You return just filtered data and you do not create/copy/move status. I imagine a Filter class like this:

class Filter
{
public:
    Filter(const valarray<double>& b, const valarray<double>& b)
    {
        // Validate parameters...
        this->a = a;
        this->b = b;

        // I omit z, I think there are very few chances you want
        // to set initial filter state...
    }

    valarray<double>& Process(const valarray<double>& input)
    {
        // Your code here, when you're done just return computed
        // result. "Final conditions" (filter status?) has not to
        // be returned, it's a private field.
    }

private:
    valarray<double> b, a, z;
};

As you can note some calculations/allocations will be performed in constructor. It won't change anything if each time you filter a different chunk of samples, using different filter, but in the other cases it will allow a big speed-up. Next step in design (to do something in C++) is to abstract away filter and its design: you may have a pure abstract base class Filter, an IirFilter implementation and one or more classes derived from FilterDesigner (let's say ChebyshevIirFilterDesigner to calculate filter coefficients). In general to mimic MATLAB approach may result in a sub-optimal implementation (because languages are different and because it will lead to a less than perfect OOP design.) If this class is a use-once function and you want to keep a function instead of a class you may drop finalConditions, do you actually use it? If you're implementing an on-line filter (for a continuous stream) then you'd better save it as private field (it's the z parameter you're using now)

You should never have using namespace std; You already use std:: consistently then you just add a possible source of conflicts.

Few more notes about performance:

In this case (I guess they do not use any esoteric math theory to perform an optimization) C++ performance should be pretty similar or even better than MATLAB. MATLAB may (I do not know what current implementation does) use SIMD/AVX instructions. Do not think you can have same throughput if you don't.

To perform a benchmark is more complex than this! At very least you should setup your machine for that and to perform calculations few thousands times to have an average and a distribution then inspect results to understand caching issues...

\$\endgroup\$
  • \$\begingroup\$ Thanks for the review. Could you elaborate on what you mean by " You return just filtered data and you do not create/copy/move status..."? \$\endgroup\$ – Phlox Midas Sep 19 '16 at 13:23
  • \$\begingroup\$ I also added few more notes, what I'd suggest is to take a look to an existing implementation from one of the DSP libraries out there. It's a great source of inspiration. If you're really go for performance you may also take a look to SIMD/AVX to speed-up calculations and, for example, OpenMP to make them parallel (if you have to filter more than one stream). \$\endgroup\$ – Adriano Repetti Sep 19 '16 at 15:01
3
\$\begingroup\$

I'm assuming you're already using the highest optimization level in your testing. If not, then you are testing apples versus oranges because MatLab functions will be optimized.

NOTE: All optimization suggestions assume that the compiler isn't doing the optimizations and may or may not improve performance, it may not be portable to other operating systems or computer architectures.

Adriano Repetti is providing some good advice about the std::copy() and about plain arrays. It may also be useful to experiment with std::array<double, SIZE> that allows you to use iterators. To make the algorithm work as it is now it only needs to copy the first 2 values of z to zOut.

C #define versus C++ const or const experssion
The code has

#define D_SCL_SECURE_NO_WARNINGS 1  

Two things about this statement

  1. It's not used anywhere in the code presented
  2. The C++ way of handling this would generally be

    constexpr int D_SCL_SECURE_NO_WARNINGS = 1;

C++ has 2 ways of defining constants, const and constexpr. Of the 2 methods constexpr is more like #define than const because constexpr is a compile time construct where const is a runtime construct. You might find more information at this stackoverflow question.

Direct Access vs Indexing
This line in the code can possibly be optimized if the compiler isn't optimizing it for you:

        y[m] = b[0] * Xm + zOut[0];

The optimization is to use an iterator or pointer to b[0] and zOut[0], or to create variable to replace b[0] and zOut[0] since b[0] never changes and zOut[0] should only change once and you can pre-calculate the value of zOut[0].

    double bZero = b[0];
    double zOutZero = bZero * x[0] + z[1] - a[1] * y[0]; // this may not be quite correct
    size_t i;
    size_t n = a.size();
    size_t mMax = y.size();
    for (size_t m = 0; m < mMax; m++)
    {
        double Xm = x[m];
        y[m] = bZero * Xm + zOutZero;
        double Ym = y[m];
        for (i = 1; i < n; i++)
        {
            zOut[i - 1] = b[i] * Xm + zOut[i] - a[i] * Ym;
        }
    }
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.