Moving average algorithm that preserves the integral and number of points

Question

I need to smooth a noisy signal. I think that the moving average algorithm would be fine. However, I want the following conditions to be satisfied:

1. The total integral (the sum of all values) before and after the moving average filtering is kept unchanged
2. The number of points before and after the moving average filtering is kept unchanged

Implementation

The implementation is based on the answer to this post. The main question here is how to process the edges of a signal. The idea is to "reflect" the signal at the edges:

$$a_{-i} = a_{i-1},\, i > 0\quad \text{at the beginning of a signal},$$ $$a_{L+i} = a_{L-i-1},\, i\geq 0\quad \text{at the end of a signal},$$ where L denotes the length of the signal.

Note

To produce more or less appropriate results this algorithm should be used with the signals which start and end with relatively flat waveform (the length of "flatness" depends on the size of the window used).

Algorithm

#include <vector>
#include <iostream>

//Input :
//       in : the original signal
//       w  : the size of the window of the moving average
//Output :
//       out : the filtered signal
void MovingAverage( const std::vector<float>& in, std::vector<float>& out, int w )
{
    //****** PREPARATION ******/
    out.clear();

    //Make the window size odd
    w += (w % 2) ? 0 : 1;

    if( (w < 1)  or (in.size() < w) )
    {
        std::cerr << "Bad input!" << std::endl;
        return;
    }

    //****** ALGORITHM ******/
    int i = 0;  //the current position

    //Firstly, average the head
    for( ; i < w/2; i++ )
    {
        float averagedValue = 0.;
        for( int j = -w/2; j < (w/2 + 1); j++ )
        {
            averagedValue += (i + j >= 0) ? in[i + j] :
                                            in[-(i + j) - 1];//(i + j) is negative here
        }
        out.push_back( averagedValue / w );
    }

    //The middle is averaged normally
    for( ; i + w/2 < in.size(); i++ )
    {
        float averagedValue = 0.;
        for( int j = -w/2; j < (w/2 + 1); j++ )
        {
            averagedValue += in[i + j];
        }
        out.push_back( averagedValue / w );
    }

    //Finally, average the tail
    for( ; i < in.size(); i++ )
    {
        float averagedValue = 0.;
        for( int j = -w/2; j < (w/2 + 1); j++ )
        {
            //(i + j) is not negative here
            averagedValue += (i + j < in.size()) ? in[i + j] :
                                                   in[in.size() - (i + j - in.size()) - 1];//(i + j - in.size()) is not negative here
        }
        out.push_back( averagedValue / w );
    }
}

Test

Answer 1

• Where you make the window length odd: w += (w % 2) ? 0 : 1; I would prefer the conditional part to be something like ((w % 2) != 0). Others will no doubt disagree with me, but I personally prefer conditionals to be boolean.

  But even better, you can get rid of the conditional entirely: w |= 0x1;. Maybe this makes some kind of not-guaranteed-to-be-fully-portable assumptions? But I'd be comfortable using it.

• Using the same index i in the multiple consecutive loops makes me think a little. On the one hand, it's really good how it doesn't duplicate code that calculates where the loops start and end. On the other hand, one doesn't normally expect a loop index, particularly one named i, to be used after the end of a loop. What if some future maintainer doesn't realize this, and inserts some other loop in here that uses i as its index? Possibly you could rename i to something else, but that might make your subsequent calculations harder to follow. In the end, I'd probably at least put a comment by the declaration of i warning that it is re-used through multiple loops.

• I don't like something that is purely DSP to write to std::cerr. Either return a status value, or throw an exception when you encounter invalid parameters. I think a status return value is preferable, but an exception may be OK if this really is an exceptional case.

• You use float to hold the intermediate averagedValue. In most cases this is fine, but depending on the nature of your data and your window size, it is possible that you'd benefit from using double for the intermediate variable.

Answer 2

Unless the middle, flat, part of the window is very small, consider changing the algorithm to maintain the state of that flat part between calls. The flat part is easily updated, by adding in the next new value and subtracting the oldest value; only the shoulders of the window need updating from scratch. That would reduce that part of the algorithm from O(n✕m)` down to just O(n).

I would recommend abstracting out the element access, so we can use a single loop once we've populated the initial sum.

There's undefined behaviour if w/2 is greater than in.size() - we end up accessing outside the valid range of in.

Consider returning the output vector as a value, rather than passing in a modifiable vector. And be sure to reserve() space for the entries we'll be adding.

Use an unsigned type as the window width, to avoid having to check for negative values. And perhaps if we're passed 0 or 1, we could just return a copy of the input, unfiltered.

Try to avoid comparing signed values with unsigned ones - e.g. for( ; i < in.size(); i++ ).

---

Here's my version of the same algorithm, addressing the points above, and also those in Eric's answer. As there's no unit tests, I haven't verified this; watch out for off-by-one errors!

#include <iostream>
#include <numeric>
#include <type_traits>
#include <vector>

//Input :
//       in : the original signal
//       w  : window size for the moving average
//Output :
//       return : the filtered signal
std::vector<float> moving_average(const std::vector<float>& in, unsigned int w)
{
    using Index = std::make_signed_t<size_t>;

    if (w < 2) {
        // tiny window - no filtering
        return in;
    }

    // Make the window size odd
    w |= 1;
    // We use half the window at each end
    auto const h = w/2;

    auto const len = static_cast<Index>(in.size());
    if (h >= len) {
        // window too big for data
        throw std::invalid_argument("w");
    }

    std::vector<float> out;
    out.reserve(in.size());

    auto const get_element
        = [len,w,&in](Index i){return i < 0 ? in[-i] : i < len ? in[i] : in[2*len-1-i]; };

    double rolling_sum = in[0] + 2 * std::accumulate(in.begin(), in.begin()+h, 0.0);
    for (Index i = 0;  i < len;  ++i) {
        out.push_back(static_cast<float>(rolling_sum / w));
        rolling_sum -= get_element(i-h);
        rolling_sum += get_element(i+h);
    }

    return out;
}