Introduction

You have this Data Acquisition System that coughs out data in real time, and you want to add a digital filter that takes each filter sample, throws out the filtered output. That’s the task of this post.

Let’s say the digital data is over your computer’s serial port, and you have set up some Matlab code to plot it already, now we need to filter each sample before storing it for plotting.

First of all you need to decide what kind of filter you want to use – low pass/bandpass etc, butterworth/elliptic etc, the filters order, cutoff frequency etc. And depending on all that use the matlab functions like butter(), ellip() to get the coefficients of the filter you want to use. You could calculate the coefficients by any other method also. This post assumes that you already have the filter coefficients and you just need some simple code to implement the filter.

Now that you have the coefficients usually denoted by the vectors a and b for the denominator and numerator respectively, let’s look at the general form of an IIR filter

Make sure that in your coefficients, $a(0) = 1$.

Use the following code to implement the above filter:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

[b_lp,a_lp] = butter(8,40/(Fs/2),'low');
global dbuffer ;
dbuffer = zeros(length(b_lp),1); % buffer for LPF
%....... Read the ADC value from the serialport...
% Pass it to the filter subfunction
filtered_value = NewFilter(b_lp, a_lp, value);
%..... Plot data....

% NewFilter subfunction
function [data] = NewFilter(b, a,value)
k = 1;
global dbuffer ;
while(k<length(b))
dbuffer(k) = dbuffer(k+1);
k=k+1;
end
dbuffer(length(b)) = 0;
k = 1;
while(k<(length(b)+1))
dbuffer(k) = dbuffer(k) + value * b(k);
k=k+1;
end

k = 1;
while(k<length(b))
dbuffer(k+1) = dbuffer(k+1) - dbuffer(1) * a(k+1);
k=k+1;
end

data = dbuffer(1);

That’s all to it! Now consider an FIR filter, this is easier because now we do not have the vector a , so just delete the Multiply and accumulate loop for the a terms in the above code.

A moving average filter can also be implemented in a similar fashion. The code looks like this:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

function [data] = MavgFilter(value)
 
k = 1;
global Dbuffer ;
global N;
% 1. shift by one element (discard the left most element)
while(k<N)
Dbuffer(k) = Dbuffer(k+1);
k=k+1;
end
Dbuffer(N) = 0;

k = 1;
while(k<N+1)
Dbuffer(k) = Dbuffer(k) + value;
k=k+1;
end

data = Dbuffer(1)/N;

Obviously, the same code can be used in any other language like C, or C# if you want to do analysis using them. If this implementation does not seem to be obvious from the equation, it’s because this is a Direct-Form Transpose II implementation. An advantage with this is that you only need to use one buffer to store intermediate results. In the more-straightforward Direct-Form Transpose I implementation, you would need two buffers - one for the numerator terms and another for the denominator terms.