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');

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.

enter image description here

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.

enter image description here

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;
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?

  • $\begingroup$ as an initial step check if you can implement a simple lowpass or highpass filter using your technique of integer implementation, if you can do so, then it's probably your notch implementation. Be warned that a notch with high selectivity requires high precision arithmetic, coefficient quantization or rounoff errors can cause problems otherwise. $\endgroup$ – Fat32 Jan 19 '17 at 11:06
  • $\begingroup$ In the first plot (floating point) you have an oscillation twice the frequency of the input signal. Why is this? $\endgroup$ – Juancho Jan 19 '17 at 11:53

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).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.