Problems with implementation of a band-stop filter on an MCU (dsPIC) using fixed-point arithmetic

Question

I'm trying to implement a band-stop (notch) filter on a microcontroller dsPIC33EP64GS506 to filter out a 100 Hz component from an input signal. The problem here is that my microcontroller doesn't have a floating-point unit, so I have to use a fixed-point arithmetic, but to be more precise, I use integer arithmetic. The system sample time is \$T_s=50~\mu\text{s}\$. Here I explain my problem in detail, and the questions are at the very end of the post. Please find attached a MATLAB code to run simulations if you want: download link on my Dropbox. Note that you don't have to be a registered Dropbox user in order to be able to download the file.

The ADC resolution is 12b, which I "increase" to 15b, not to increase the masurement resolution (which cannot be done), but to increase the filter resolution:

v_in = ADCBUF0<<3;

The notch filter transfer function in s-domain is given as:

$$G_s(s) = \frac{s^2 + 2\zeta_1\omega_n s + \omega_n^2}{s^2 + 2\zeta_2\omega_ns + \omega_n^2},$$

where filter parameters \$\zeta_1\$, \$\zeta_2\$, and \$\omega_n\$ fully determine the filter. If we want to filter out a 100 Hz component, then we set \$\omega_n=2\pi\cdot100~\text{s}^{-1}\$. As for other parameters, without detailed explanations, I use \$\zeta_1=0.001\$ and \$\zeta_2=1\$. To implement this filter on a digital system, we need to discretize the transfer function in s-domain, and for that, I use a Tustin discretization method.

$$G_z(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}}.$$

The corresponding recursive equation, as implemented on a digital system, is:

$$y(k) = b_0 u(k) + b_1 u(k-1) + b_2 u(k-2) - a_1 y(k-1) - a_2 y(k-2).$$

Here is the MATLAB code to get the transfer function in z-domain:

Ts = 50e-6;
zeta1 = 0.001;
zeta2 = 1;
wn = 100*(2*pi);
Gs = tf([1 2*zeta1*wn wn^2], [1 2*zeta2*wn wn^2]);
Gz = c2d(Gs, Ts, 'tustin');
bodeplot(Gz);

Please see below a frequency characteristics of the filter in z-domain generated by MATLAB. As one can see, this type of a filter is very sensitive in terms of a frequency response. Even smallest changes in parameters could cause completely different behavior, e.g., an unexpected gain, phase shift etc.

Now, since I don't have a floating-point unit at disposal, I use a well-known method called binary scaling. Any decimal number \$ d \in \mathbb{R} \$ can be represented as \$\frac{\lfloor d\cdot2^r \rceil}{2^r}\$, where \$r\in\mathbb{N}\$, and \$\lfloor\cdot\rceil\$ is a round-to-nearest-integer function. Binary scaling is a method used when we want to avoid using a division, which is very expensive operation in terms of required CPU cycles (typically 18 cycles), since binary shift can produce the same results in much less cycles (typically 1 cycle). For this example I used \$r=15\$, which is the highest precision that I could use considering available bits to do a multiplication (a result should be within 32 bits). The corresponding recursive equation using integer arithmetic is implemented as follows:

yk0 = ( B0*uk0+B1*uk1+B2*uk2 - A1*yk1-A2*yk2 + (1<<14) ) >> 15;
uk2 = uk1; uk1 = uk0;
yk2 = yk1; yk1 = yk0;

The term (1<<14) is used to ensure that the result is rounded to nearest integer after the bit shift operation. Note that >>15 is practically an integer division by 32768, but we aware - it is not completely equivalent! Bit shift always rounds to minus infinity, while integer division always rounds to zero. The filter parameters values are as follows: B0=31771, B1=-63509m B2=31769, A1=-63509, A2=30772. This can also be found in a provided MATLAB code.

I was really shocked when I realized that the integer variant of the filter doesn't work at all. Please see below the responses of three different filters used to remove a 100 Hz component from the input signal. From top to bottom: (1) notch filter in floating point implementation, (2) notch filter in integer arithmetic implementation using \$r=15\$, (3) a simple moving average filter.

Here is how I've implemented a moving average filter in C:

uk200 = window[ind];
window[ind] = uk0;
win_sum = win_sum - uk200 + uk0;
yk0 = (((win_sum+(1<<3))>>4)*5243+(1<<15))>>16;
ind++;
if (ind==200) ind=0;

Questions to the community.

Since I don't have that much of an experience in digital filtering, can someone confirm to me are these results expected. Is really that problematic to implement a notch filter on a digital system using integer arithmetic only? What digital filter is typically used to remove a certain frequency? As I see from these results, the moving average filter outperforms both of notch filter implementations. Are there any downsides for the moving average filter that I should be aware of, except for obviously increased memory demand?

Answer 1

There are two problems with digital filters, one is dynamic range (ability to represent a range of numbers), the other is quantization noise from rounding.

Integer filters have a much reduced dynamic range vs floating point filters and also have more quantization noise.

The filter needs to be check to ensure its not saturating registers. (The filter is always multiplying a big number by a big constant, the mulitplier will overflow and this will cause instability in the filter or noise). There are ways to overcome this, the coefficents can be adjusted to get the same response but avoiding.

See this resource for more information.

One thing to try (if its possible with your compiler\hardware) is to increase the size of the variables from 32 to 64 and see if that helps.

Another thing you will want to check (and I haven't done this yet, I need to figure out what form your filter is before I can convert it to a different form) is the filter order and form. Different forms have different results. The filter order you have is a second order filter. The form looks like a second order section of some kind but I'm not sure (right now)

Digital filters can have a higher order, a notch filter can be designed (or highpass or lowpass, whatever) filter with a second order filter or more, the higher the order the more constants, and multiplies and adds (more computational resources) but the better the filter is approximated.

Its probably better to design the filter via digital methods, there are papers describing how to calculate the coefficents, but there are also calculators online and of course matlab, but use the filter design tool. Make sure that whatever method you use your calcualted constants match the form (bi-quad second order sections, Direct form 1, Direct form II ect)