Moving average algorithm that conserves integral

Question

Suppose I have point array $y_i$ of size $N$. How to implement moving average algorithm that conserves quantity $$ I = \sum_{i=1}^{N}y_i $$ NOTE: I don't want time shift so I would prefer to use symmetric SMA in the middle: $$ \hat{y}_i = \frac{1}{n}\sum_{k = i-\frac{n-1}{2}}^{i+\frac{n-1}{2}}y_k $$ where $\hat{y}_i$ - filtered point, $n$ - window size of SMA.

It is easy to show that the above algorithm conserves $I$ "in the middle" but problem is when processing edges of a signal. So question is how to process edges of a signal in order to conserve $I$?

After several failures to construct simple method to do it I am here.

Answer 1

Why do you want to preserve the sum? Is it to make the filter with the same gain everywhere?

If you want to keep the sum the same, you need to pad the beginning and the end of the signal as: $$ .... ,y_2,y_1,y_0 ---- y_0,y_1, y_2,...., y_n ---- y_n, y_{n-1},y_{n-2}... $$ Example: suppose that you have 6 samples

$$ y_0, y_1, y_2, y_3,y_4, y_5 $$

and you want to to do a moving average with length 5:

Special treatment at the beginning:

$$ \hat y_0 = (y_1 + y_0 + y_0 + y_1 + y_2)/5, $$ $$ \hat y_1 = (y_0 + y_0 + y_1 + y_2 + y_3)/5 $$

Normal filtering:

$$ \hat y_2 = (y_0 + y_1 + y_2 + y_3 + y_4)/5, $$ $$ \hat y_3 = (y_1 + y_2 + y_3 + y_4 + y_5)/5, $$

Special treatment at the end:

$$ \hat y_4 = (y_2 + y_3 + y_4 + y_5 + y_5)/5, $$ $$ \hat y_5 = (y_3 + y_4 + y_5 + y_5 + y_4)/5, $$

Answer 2

Adding to other comments: I believe that you are describing a filtering operation that avoids tapering at the edges, similar to what is often desired in image filtering.

There are built-in MATLAB functions that aim to «match initial condition» or »avoid edge bias», eg filtfilt and xcorr, but I gravitate towards two solutions. One using a straight forward convolution and compensation gain at the edges:

%% task specification
x = pi*ones(10,1);
ma_length = 4;

k = ones(ma_length,1)/ma_length;
y = conv(x, k,'full');
eq_gain = ma_length./(1:ma_length)';
y2 = y;
y2(1:ma_length) = y(1:ma_length) .* eq_gain;
y2(end-ma_length+1:end) = y(end-ma_length+1:end) .* eq_gain(end:-1:1)

Or, alternatively, using a CIC-type accumulator/comb filter:

y_tmp = 0;
for n = 1:ma_length
    y_tmp = y_tmp + x(n);
    y4(n) = y_tmp/n;
end
for n = (ma_length+1) : length(x)
    y_tmp = y_tmp - y4(n-ma_length) + x(n);
    y4(n) = y_tmp/ma_length;
end
for n = (length(x)+1) : (length(x)+ma_length-1)
    y_tmp = y_tmp - y4(n-ma_length);
    y4(n) = y_tmp/((length(x)+ma_length-n));
end

Code have not been tested other than observing that I got 13 ones. May well contain bugs