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\textrm{ 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\textrm{ 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\textrm{ 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, with a $10~\text{ms}$ ($200$ samples) window.

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

The denominator (recursive coefficients Ai) look OK: the poles of your system are at 45 degree angles ($\pi/4$), with magnitude 0.68 (which is not very aggressive for a notch filter; in my opinion they should be more like 0.9).

But your numerator has its roots very near $z=1$, which corresponds to frequency 0 instead of the desired $\pi/4$ for implementing the notch.

This type of notch filter is quite straightforward: 2 zeros on the unit circle at the desired frequency, and 2 poles near each zero just inside the unit circe.

Instead of complicating matters designing in continuous time and then applying the bilineal transform, I suggest designing directly in discrete time.

Assuming your notch frequency is 1/8th the sample frequency, as suggested by the position of the poles, then the numerator should be:

$$\left(1 - e^{j\pi/4}z^{-1}\right)\left(1 - e^{-j\pi/4}z^{-1}\right) = 1 - \sqrt{2}z^{-1} + z^{-2}$$

So your coefficients for the numerator should be $[1, -\sqrt{2}, 1] \times 65536$ (instead of the $[1, -2, 1]$ you have now).