Document fait avec Nvu Document made with Nvu

Technical pages
Numerical filtering



The goal of this page is not to make a lesson on the numerical filtering but to show how to do it in the FORTH system described in this web site.

Numerical filtering is one of the main applications of a real time kernel. To optimize the computing times, anything of such as to work with integers rather than floating numbers especially when we want to use a processor in the moderate price. This implies to use a filtering cell "biquad" allowing to realize all the second order filters and and being able to be put in cascade with other cells for higher order filters:

for a biquad cell, we have

H(z) = (a0 + a1*z^-1 + a2*z^-2)/(1 + b1*z^-1 + b2*z^-2)

corresponding to the algorithm

y(n) = BIQUAD(x(n)) with the successive operations

w(n) = x(n) + b1*w(n-1) + b2*w(n-2)
y(n) = a0*w(n) + a1*w(n-1) + a2*w(n-2)
w(n-2) = w(n-1)
w(n-1) = w(n)
for a more than order 2 filter, we have

H'(z) = H1(z)*H2(z)*...*Hm(z)

corresponding to the algorithm

y(n) = (BIQUADm(...(BIQUAD2(BIQUAD1(x(n)))...))

It is necessary to define now the biquad filter parameters with optimization of the computing precision and time. There are in fact 5 coefficients a0, a1, a2, b1 ,b2 and 2 variables w(n-1), w(n-2). Each coefficient is multiplied by a 2^N factor in order to the biggest absolute value is the most possible closed bu less than 32768 because this number is equal to -32768 in the 16 bits signed format. For a gived N, 2^N is corresponding to the value unity. All those parameters are furnished through the address of a 8 16 bits integers table:

adr_biquad --> N, a2<<N, b2<<N, w(n-2), a1<<N, b1<<N, w(n-1), a0<<N
w(n-1) e w(n-2) are initialized to 0, the computing algorithm is the following:

w(n) = (x(n)<<N + b1<<N*w(n-1) + b2<<N*w(n-2))>>N
y(n) = (a0<<N*w(n) + a1<<N*w(n-1) + a2<<N*w(n-2))>>N
w(n-2) = w(n-1)
w(n-1) = w(n)

To determine the value of N, we have to search the coefficien with the great absolute value then multiply it by 2 while this value is less than 16383 or divide it by 2 while this value is more than 32767. We add to N, initialized to 0, 1 for each multiplication by 2 or -1 for each division by 2.

Each biquad filter used in an application is so defined by 8 integers table (6 16 bits words + 2 32 bits words) and only one instruction is necessary. A FORTH source code of BIQUAD instruction is avalaible in the file biquad.txt.


Only one instruction

x[n],adr_biquad BIQUAD y[n]

Computing of the output y[n] of the biquad filter described by adr_biquad and with the input value x[n]:
- x[n] is the input sample,
- adr_biquad is the address of the filter parameters,
adr_biquad --> N a2<<N b2<<N w[n-2] a1<<N b1<<N w[n-1] a0<<N
- y[n] is the output sample computed as following,
w[n] = [x[n]<<N + b1<<N*w[n-1] + b2<<N*w[n-2]]>>N
y[n] = [a0<<N*w[n] + a1<<N*w[n-1] + a2<<N*w[n-2]]>>N
w[n-2] = w[n-1]
w[n-1] = w[n]


We can start from the transfert function of a Butterworth second order low pass filter for example:

H(p) = 1/(p^2 + 1,414*p + 1)

following the bi-linear transformation

p = K*(1 - z^-1)/(1 + z^-1) et K = 1/(tang(PI*fc/fe))

with the cutting frequency fc and the sampling frequency fe, we obtains

H(z) = (1 + 2*z^-1 + z^-2)/((K^2 + 1,414*K + 1)+(2*(1 - K^2)*z^-1 + (K^2 - 1,414*K + 1)*z^-2)

then the algorithm

w(n) = x(n) - [(2*(1 - K^2)*w(n-1) + (K^2 - 1,414*K + 1)*w(n-2))/(K^2 + 1,414*K + 1)]
y(n) = (w(n) + 2*w(n-1) + w(n-2))/(K^2 + 1,414*K + 1)
w(n-2) = w(n-1)
w(n-1) = w(n)

If fc = 50 Hz and fe = 1 KHz, then K = 6,314 so:

w(n) = x(n) + 1,561*w(n-1) - 0,639*w(n-2)
y(n) = 0,02*w(n) + 0,04*w(n-1) + 0,02*w(n-2)

this is corresponding to the following table: 14, 329, -10508, 0, 658, 25576, 0, 329.

The table is created with the following FORTH sequence:

14 ,
329 , -10508 , 0 2,
658 , 25576 , 0 2,
329 ,

w(n-2) et w(n-1) are stored with 32 bits words in order to optimize the computing precision.
To compute a sample value:


With fe = 1 KHz, this action mus be realized each milliseconde.


To illustrate this page, I have realized the software biquad.txt and used the system window of the FORTH core emulator. and the curves displaying application.

There are 5 filters with a cutting frequency of 50 Hz and a sampling frequency of 1 KHz. The detailed calculates are avalaibles in the software source file.

First order low pass:

First order high pass:

Band pass with quality coefficien of 1:

Band stop with quality coefficien of 1:

Second order low pass:

Second order high pass: