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

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?

• Integer arithmetic (or rather, fixed point) is generally preferred for DSP, floats bring some difficult to deal with problems of rounding. And you can build quite good notches using FIR filters (which are a generalisation of the moving average). My preference is to use a language that supports fixed point directly, at least for prototyping teh filter - porting to integer isn't difficult. – user_1818839 Jan 19 '17 at 0:21
• What are the values that you used for B0, B1, B2, A1 and A2 on above C code? Are you scaling these coefficients accordingly? Are you also scaling the value of the ADC (or the input)? However, programming the DSPIC in assembly has some advantages when using fixed point, such as the access to the MAC instructions and the 40-bit accumulator. You do not have to port all the code, but only the filter related functions. In addition, Microchip provides a library for DSP operations, even working in C. – Dirceu Rodrigues Jr Jan 19 '17 at 0:35
• @DirceuRodriguesJr I've updated my question with additional information. Values for parameters are given immediately after they're introduced. The ADC works on 12b resolution, which I "increase" to 15b to increase the filter resolution. – Marko Jan 19 '17 at 0:50
• So if you're using r = 15, the original coefficient for -63509 will be -1.93814 ... (not scaled). I think this would require a 32 bit Q1.14 representation. In other words: 1 bit of sign, 1 bit for the integer part and 14 bits for the fractional part (with r = 14). But remember that Q1.14 x Q1.14 leads to a result with 2 sign bits (repeated), 2 bits for the integer part and 28 bits on the fractional part; which will require appropriate shifts for correction. – Dirceu Rodrigues Jr Jan 19 '17 at 1:51
• I'm not using a fixed-point arithmetic, but integer arithmetic. In other words, I keep the position of a point in my head :) I've tried to increase the resolution by additional 5 bits, but it didn't help. – Marko Jan 19 '17 at 1:56

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.