Demonstration for a naive implementation of a moving average filterΒΆ

../../_images/sphx_glr_fig-examp-moving-average-stepresponse_001.png
import numpy as np
import matplotlib.pyplot as plt



def exponential_moving_average_filter_recursive(y, a, b=0, gamma=1):
    """
    Moving average filter with sliding window borders a and b.

    Note: this is a naive (slow) implementation and for demonstration purpose only.

    :param y: signal
    :param a: left window boundary
    :param b: right window boundary (optional)
    :param gamma: exponential decay factor 0..1
    :return: filtered signal
    """

    y_avg = np.zeros_like(y) # create empty array
    K = len(y)

    if (b>0): # correct initialization of y_sum for non-causal windows
        pass # ... here whe should initialize y_avg[0] correctly ...

    mu = (gamma**(b)-gamma**(a-1))/(1-1/gamma)

    for k in range(K):
        y_sum = 0

        if (k+a >= 0) and (k+a < K):
             y_sum = y_sum - (gamma**(a-1))*y[k+a]

        if (k+1+b >= 0) and (k+1+b < K):
             y_sum = y_sum + (gamma**(b))*y[k+1+b]

        if (k+1 >= 0) and (k+1 < K):
             y_avg[k+1] =  y_avg[k]* gamma**(-1)  + y_sum / mu

    return y_avg




# Load Signal
K = 200 # number of samples to load
k = range(K)

# generate rectangular pulase
y = np.concatenate((np.zeros(int(K/4)), np.ones(int(K/2)), np.zeros(int(K/4))))


# --- moving average filter and peak detection
L=20 # averaging window size
gamma = 1.1

y_ma1 = exponential_moving_average_filter_recursive(y, a=-2*L, b=0, gamma=1.00001) # non-symmetric filter

y_ma2 = exponential_moving_average_filter_recursive(y, a=-L, b=L, gamma=1.001) # symmetric filter

y_ma3 = exponential_moving_average_filter_recursive(y, a=-200*L, b=0, gamma=gamma) # non-symmetric filter



# Plots
_, axs = plt.subplots(1, 3, sharex='all', figsize=(10, 2))
plt.subplots_adjust(hspace=0.3)

ax = axs[0]
ax.plot(k, y, lw=0.5, c='tab:gray', label='$y$')
ax.plot(k, y_ma1, lw=1.0, c='b', label='${y}_{avg}$')
ax.set(ylabel=r'$y_k$')
ax.set(xlabel='k')
ax.set_title('asymmetric, L='+str(2*L))
# ax.legend(loc=4)


ax = axs[1]
ax.plot(k, y, lw=0.5, c='tab:gray', label='$y$')
ax.plot(k, y_ma2, lw=1.0, ls='-', c='b', label='${y}_{avg}$')
ax.set(xlabel='k')
ax.set_title('symmetric, L='+str(L))
# ax.legend(loc=4)


ax = axs[2]
ax.plot(k, y, lw=0.5, c='tab:gray', label='$y$')
ax.plot(k, y_ma3, lw=1.0, ls='-', c='b', label='${y}_{avg}$')
ax.set(xlabel='k')
ax.set_title('exp. asymmetric, gamma='+str(gamma))
# ax.legend(loc=4)


plt.show()

Total running time of the script: ( 0 minutes 0.183 seconds)

Gallery generated by Sphinx-Gallery