1
\$\begingroup\$

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.

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

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

\$\endgroup\$
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ – Brian Drummond Jan 19 '17 at 0:21
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ – Dirceu Rodrigues Jr Jan 19 '17 at 0:35
  • \$\begingroup\$ @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. \$\endgroup\$ – Marko Gulin Jan 19 '17 at 0:50
  • \$\begingroup\$ 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. \$\endgroup\$ – Dirceu Rodrigues Jr Jan 19 '17 at 1:51
  • \$\begingroup\$ 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. \$\endgroup\$ – Marko Gulin Jan 19 '17 at 1:56
1
\$\begingroup\$

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)

\$\endgroup\$
  • \$\begingroup\$ I've actually emulated the integer arithmetic implementation in MATLAB by using floating-point arithmetic with rounding. Because of that, I can rule out for possible overflows. I can't increase the size to 64b, as my compiler (for the MCU) doesn't support that. What I wanted to know - can I use the moving average filter, i.e., is this common practice to remove a certain frequency component? Thank you for the provided resource, I'll check it out! \$\endgroup\$ – Marko Gulin Jan 19 '17 at 1:14
  • \$\begingroup\$ By moving average, you mean FIR filters then yes. There are many different realizations of filters. Depending on how you design your filter (you can design filters with more coefficients, the more coefficients you have the better it approximates the filter) This is a good resource There are different structures, each structure has its pro's and cons (some being more computationally expensive, some being less precise on their filtering) \$\endgroup\$ – Voltage Spike Jan 19 '17 at 7:57
  • \$\begingroup\$ I've updated my post to also include a moving average filter implementation. You can find my implementation just before the questions (at the end of the post). \$\endgroup\$ – Marko Gulin Jan 19 '17 at 8:05
  • 1
    \$\begingroup\$ If you don't need the higher frequencies a moving average filter will work fine, they have similar frequency responses to a low pass filter. gaussianwaves.com/2010/11/moving-average-filter-ma-filter-2 \$\endgroup\$ – Voltage Spike Jan 19 '17 at 21:16

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.