Understand Moving Average Filter with Python & Matlab

The moving average filter is a simple Low Pass FIR (Finite Impulse Response) filter commonly used for smoothing an array of sampled data/signal. It takes L samples of input at a time and takes the average of those L-samples and produces a single output point. It is a very simple LPF (Low Pass Filter) structure that comes handy for scientists and engineers to filter unwanted noisy component from the intended data.

As the filter length increases (the parameter L) the smoothness of the output increases, whereas the sharp transitions in the data are made increasingly blunt. This implies that this filter has excellent time domain response but a poor frequency response.

The MA filter performs three important functions:

1) It takes L input points, computes the average of those L-points and produces a single output point
2) Due to the computation/calculations involved, the filter introduces a definite amount of delay
3) The filter acts as a Low Pass Filter (with poor frequency domain response and a good time domain response).

Implementation

The difference equation for a L-point discrete-time moving average filter with input represented by the vector \mathbf{x} and the averaged output vector \mathbf{y}, is 

\[y[n] = \displaystyle{\frac{1}{L} \sum_{k=0}^{L-1}x[n-k]} \quad\quad (1) \]

For example, a 5-point Moving Average FIR filter takes the current and previous four samples of input and calculates the average. This operation is represented as shown in the Figure 1 with the following difference equation for the input output relationship in discrete-time.

\[\begin{aligned} y[n] &= \displaystyle{\frac{1}{5} \left(x[n] + x[n-1] + x[n-2] + x[n-3] + x[n-4] \right) } \\ &= 0.2 \left(x[n] + x[n-1] + x[n-2] + x[n-3] + x[n-4] \right) \end{aligned} \quad\quad (2) \]
Figure 1: Discrete-time 5-point Moving Average FIR filter

The unit delay shown in Figure 1 is realized by either of the two options:

  • Representing the input samples as an array in the computer memory and processing them
  • Using D-Flip flop shift registers for digital hardware implementation. If each discrete value of the input x is represented as a 12-bit signal line from ADC (analog to digital converter), then we would require 4 sets of 12-Flip flops to implement the 5-point moving average filter shown in Figure 1.

Z-Transform and Transfer function

In signal processing, delaying a signal x[n] by k sample period (unit delay) is equivalent to multiplying the Z-transform X[z] by z^{-k}. By applying this idea, we can find the Z-transform of the 5-point moving average filter in equation (2) as

\[ \begin{aligned} y[n] &= 0.2 \left(x[n] + x[n-1] + x[n-2] + x[n-3] + x[n-4] \right) \\ Y[z] &= 0.2 \left(z^0 + z^{-1} + z^{-2} + z^{-3} + z^{-4} \right) X[z] \\ Y[z] &= 0.2 \left(1 + z^{-1} + z^{-2} + z^{-3} + z^{-4} \right) X[z] \quad\quad (3) \end{aligned}\]

Similarly, the Z-transform of the generic L-sample Moving Average filter of equation (1) is 

\[x Y(z) = \left(\displaystyle{\frac{1}{L} \sum_{k=0}^{L-1} z^{-k}} \right) X(z) \quad\quad (4) \]

The transfer function describes the input-output relationship of the system and for the L-point Moving Average filter, the transfer function is given by

\[H(z) = \displaystyle{\frac{Y(z)}{X(z)}} =\displaystyle{\frac{1}{L} \sum_{k=0}^{L-1} z^{-k}}  \quad\quad (5)\]

Simulating the filter in Matlab and Python

In Matlab, we can use the filter function or conv (convolution) to implement the moving average FIR filter.

In general, the Z-transform Y(z) of a discrete-time filter’s output y(n) is related to the Z-transform X(z) of the input x(n) by

\[H(z) = \displaystyle{\frac{Y(z)}{X(z)}} = \displaystyle{\frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + \cdots + b_n z^{-n}}{a_0 + a_1 z^{-1} + a_2 z^{-2} + \cdots + a_m z^{-m}}} \quad\quad (6) \]

where, \{b_i\} and \{a_i\} are the filter coefficients and the order of the filter is the maximum of n and m

For implementing equation (6) using the filter function, the Matlab function is called as

B = [b_0, b_1, b_2, ..., b_n] %numerator coefficients
A = [a_0, a_1, a_2, ..., a_m] %denominator coefficients
y = filter(B,A,x) %filter input x and get result in y

We can note from the difference equation and transfer function of the L-point moving average filter, that following values for the numerator coefficients \{b_i\} and denominator coefficients \{a_i\}.

\[\begin{aligned} \{b_i\} &= \{1/L, 1/L, …., 1/L\} &,i=0,1,\cdots,L-1 \\ \{a_i\} &= {1} &,i=0 \end{aligned} \quad\quad (7)\]

Therefore, the 5-point moving average filter can be coded as

B = [0.2, 0.2, 0.2, 0.2, 0.2] %numerator coefficients
A = [1] %denominator coefficients
y = filter(B,A,x) %filter input x and get result in y

The numerator coefficients for the moving average filter can be conveniently expressed in short notion as shown below

 L = 5
B = ones(1,L)/L %numerator coefficients
A = [1] %denominator coefficients
x = rand(1,10) %random samples for x
y = filter(B,A,x) %filter input x and get result in y

When using conv function to implement the moving average filter, the following code can be used

L = 5;
x = rand(1,10) %random samples for x;
y = conv(x,ones(1,L)/L)

There exists a difference between using conv function and filter function for implementing an FIR filter. The conv function gives the result of complete convolution and the length of the result is length(x)+ L -1. Whereas, the filter function gives the output that is of same length as that of the input x.

In python, the filtering operation can be performed using the lfilter and convolve functions available in the scipy signal processing package. The equivalent python code is shown below.

import numpy as np
from scipy import signal
L=5 #L-point filter
b = (np.ones(L))/L #numerator co-effs of filter transfer function
a = np.ones(1)  #denominator co-effs of filter transfer function
x = np.random.randn(10) #10 random samples for x
y = signal.convolve(x,b) #filter output using convolution

y = signal.lfilter(b,a,x) #filter output using lfilter function

Pole-zero plot and frequency response

A pole-zero plot for a filter transfer function H(z), displays the pole and zero locations in the z-plane. In the pole-zero plot, the zeros occur at locations (frequencies) where H(z) \rightarrow 0 and the poles occur at locations (frequencies) where H(z) \rightarrow \infty. Equivalently, zeros occurs at frequencies for which the numerator of the transfer function in equation 6 becomes zero and the poles occurs at frequencies for which the denominator of the transfer function becomes zero.

In a pole-zero plot, the locations of the poles are usually marked by cross (\times) and the zeros are marked as circles (\circ). The poles and zeros of a transfer function effectively define the system response and determines the stability and performance of the filtering system.

In Matlab, the pole-zero plot and the frequency response of the L-point moving average can be obtained as follows.

L=11
zplane([ones(1,L)]/L,1); %pole-zero plot
w = -pi:(pi/100):pi; %to plot frequency response
freqz([ones(1,L)]/L,1,w); %plot magnitude and phase response
Figure 2: Pole-Zero plot for L=11-point Moving Average filter
Figure 3: Magnitude and phase response of L=11-point Moving Average filter

The magnitude and phase frequency responses can be coded in Python as follows

[updated] See my Google colab for full Python code

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

L=11 #L-point filter
b = (np.ones(L))/L #numerator co-effs of filter transfer function
a = np.ones(1)  #denominator co-effs of filter transfer function

w, h = signal.freqz(b,a)
plt.subplot(2, 1, 1)
plt.plot(w, 20 * np.log10(abs(h)))
plt.ylabel('Magnitude [dB]')
plt.xlabel('Frequency [rad/sample]')

plt.subplot(2, 1, 2)
angles = np.unwrap(np.angle(h))
plt.plot(w, angles)
plt.ylabel('Angle (radians)')
plt.xlabel('Frequency [rad/sample]')
plt.show()
 Impulse response, Pole-zero plot, magnitude and phase response of L=11 moving average filter
Figure 4: Impulse response, Pole-zero plot, magnitude and phase response of L=11 moving average filter

Case study:

Following figures depict the time domain & frequency domain responses of a L-point Moving Average filter. A noisy square wave signal is driven through the filter and the time domain response is obtained.

Processing a signal through Moving average filter of various lengths
Figure 5: Processing a signal through Moving average filter of various lengths

On the first plot, we have the noisy square wave signal that is going into the moving average filter. The input is noisy and our objective is to reduce the noise as much as possible. The next figure is the output response of a 3-point Moving Average filter. It can be deduced from the figure that the 3-point Moving Average filter has not done much in filtering out the noise. We increase the filter taps to 10-points and we can see that the noise in the output has reduced a lot, which is depicted in next figure.

With L=51 tap filter, though the noise is almost zero, the transitions are blunted out drastically (observe the slope on the either side of the signal and compare them with the ideal brick wall transitions in the input signal).

Frequency response of moving average filters of various lengths
Figure 6: Frequency response of moving average filters of various lengths

From the frequency response of lower order filters (L=3, L=10), it can be asserted that the roll-off is very slow and the stop band attenuation is not good. Given this stop band attenuation, clearly, the moving average filter cannot separate one band of frequencies from another. As we increase the filter order to 51, the transitions are not preserved in time domain. Good performance in the time domain results in poor performance in the frequency domain, and vice versa. Compromise need for optimal filter design.

Rate this article: PoorBelow averageAverageGoodExcellent (52 votes, average: 4.33 out of 5)

References:

[1] Filtering – OpenCourseWare – MIT.↗
[2] HC Chen et al.,”Moving Average Filter with its application to QRS detection”,IEEE Computers in Cardiology, 2003.↗

Similar topics

Essentials of Signal Processing
● Generating standard test signals
 □ Sinusoidal signals
 □ Square wave
 □ Rectangular pulse
 □ Gaussian pulse
 □ Chirp signal
Interpreting FFT results - complex DFT, frequency bins and FFTShift
 □ Real and complex DFT
 □ Fast Fourier Transform (FFT)
 □ Interpreting the FFT results
 □ FFTShift
 □ IFFTShift
Obtaining magnitude and phase information from FFT
 □ Discrete-time domain representation
 □ Representing the signal in frequency domain using FFT
 □ Reconstructing the time domain signal from the frequency domain samples
● Power spectral density
Power and energy of a signal
 □ Energy of a signal
 □ Power of a signal
 □ Classification of signals
 □ Computation of power of a signal - simulation and verification
Polynomials, convolution and Toeplitz matrices
 □ Polynomial functions
 □ Representing single variable polynomial functions
 □ Multiplication of polynomials and linear convolution
 □ Toeplitz matrix and convolution
Methods to compute convolution
 □ Method 1: Brute-force method
 □ Method 2: Using Toeplitz matrix
 □ Method 3: Using FFT to compute convolution
 □ Miscellaneous methods
Analytic signal and its applications
 □ Analytic signal and Fourier transform
 □ Extracting instantaneous amplitude, phase, frequency
 □ Phase demodulation using Hilbert transform
Choosing a filter : FIR or IIR : understanding the design perspective
 □ Design specification
 □ General considerations in design

Books by the author

Wireless Communication Systems in Matlab
Wireless Communication Systems in Matlab
Second Edition(PDF)

PoorBelow averageAverageGoodExcellent (159 votes, average: 3.81 out of 5)

Digital modulations using Python
Digital Modulations using Python
(PDF ebook)

PoorBelow averageAverageGoodExcellent (122 votes, average: 3.60 out of 5)

digital_modulations_using_matlab_book_cover
Digital Modulations using Matlab
(PDF ebook)

PoorBelow averageAverageGoodExcellent (125 votes, average: 3.69 out of 5)

Hand-picked Best books on Communication Engineering
Best books on Signal Processing

Post your valuable comments !!!