2
$\begingroup$

I have always used high-performance 32-bit microcontrollers with hardware FPU support. To calculate the coefficients of digital filters, I used the fdatool package in Matlab. I used single-precision floating point coefficients for digital filters. Next, in accordance with this figure, I implemented digital filters: enter image description here Moreover, sometimes it was possible to do a little optimization, for example with the LP filter I did this:

void lp_filter (float *x ,float *y,float *aa ,float g ,float vh)
{
    x[2]=x[1];
    x[1]=x[0];
    x[0]=vh*g;
    y[2]=y[1];
    y[1]=y[0];
    y[0]=(x[0] + x[1] + x[1] + x[2] - aa[0]*y[1] - aa[1]*y[2]); //b0=1 b1=2 b2=1
}

I knew that the b coefficients would always be (1,2,1), so I simply could not do unnecessary multiplication operations. Perhaps the compiler itself could do this, but I'm not sure about that...

Now I have a low-performance 32-bit microcontroller without an FPU. I wanted to implement digital filters, but with integer mathematics. The fdatool package also allows you to calculate integer filter coefficients, but I don’t understand how to use them... For example, here are the coefficients for a 500Hz 6th order low pass filter:

#define MWSPT_NSEC 7
const int NL[MWSPT_NSEC] = { 1,3,1,3,1,3,1 };
const int32_T NUM[MWSPT_NSEC][3] = {
  {
     265839664,           0,           0 
  },
  {
    2147483647,  2147483647,  2147483647 
  },
  {
     209661133,           0,           0 
  },
  {
    2147483647,  2147483647,  2147483647 
  },
  {
     186862344,           0,           0 
  },
  {
    2147483647,  2147483647,  2147483647 
  },
  {
    2147483647,           0,           0 
  }
};
const int DL[MWSPT_NSEC] = { 1,3,1,3,1,3,1 };
const int32_T DEN[MWSPT_NSEC][3] = {
  {
    2147483647,           0,           0 
  },
  {
    2147483647, -2147483648,  1483049895 
  },
  {
    2147483647,           0,           0 
  },
  {
    2147483647, -2024667000,   715827883 
  },
  {
    2147483647,           0,           0 
  },
  {
    2147483647, -1804502424,   404468153 
  },
  {
    2147483647,           0,           0 
  }
};

I choose the format of the filter coefficients as 32 bit signed integer. If I follow the same implementation I used for floating point, the input value may overflow 32 bits... Perhaps I need to use the coefficients as 16 bit signed integers or store the results in 64 bit registers? I tried using the int64_t format, but the filters don't work correctly. I can't understand why... Also, in the integer implementation of digital filters, I don’t understand how I should work with the final result? For example, how can I understand where the comma would be if it were a float number?


For digital filters I used single precision floating point coefficients. Here is an example of such coefficients that I obtained in Matlab:

#define MWSPT_NSEC 3
const int NL[MWSPT_NSEC] = { 1,3,1 };
const real32_T NUM[MWSPT_NSEC][3] = {
  {
    0.01185768284,              0,              0 
  },
  {
                1,              2,              1 
  },
  {
                1,              0,              0 
  }
};
const int DL[MWSPT_NSEC] = { 1,3,1 };
const real32_T DEN[MWSPT_NSEC][3] = {
  {
                1,              0,              0 
  },
  {
                1,   -1.669203162,   0.7166338563 
  },
  {
                1,              0,              0 
  }
};

Based on these coefficients, I implemented a digital filter like this:

float sek1_x[3],sek1_y[3],sek1_a[3]={1,   -1.669203162,   0.7166338563},sek1_b[3]={1,              2,              1};
void lp_filter (float *x ,float *y,float *aa ,float *bb ,float g ,float in)
{
    x[2]=x[1];
    x[1]=x[0];
    x[0]=in*g;
    y[2]=y[1];
    y[1]=y[0];
    y[0]=(x[0]*bb[0] + x[1]*bb[1] + x[2]*bb[2] - aa[1]*y[1] - aa[2]*y[2]);
} 
//Calling the filter function:
lp_filter(sek1_x,sek1_y,sek1_a,sek1_b,0.01185768284,adc_res);

The fdatool package also allows you to calculate integer filter coefficients. Here is an example of such coefficients for exactly the same filter, which I obtained in Matlab:

#define MWSPT_NSEC 3
const int NL[MWSPT_NSEC] = { 1,3,1 };
const int32_T NUM[MWSPT_NSEC][3] = {
  {
      25464180,           0,           0 
  },
  {
    2147483647,  2147483647,  2147483647 
  },
  {
    2147483647,           0,           0 
  }
};
const int DL[MWSPT_NSEC] = { 1,3,1 };
const int32_T DEN[MWSPT_NSEC][3] = {
  {
    2147483647,           0,           0 
  },
  {
    2147483647, -2147483648,  1538959525 
  },
  {
    2147483647,           0,           0 
  }
};

Here is the code for implementing a digital filter with integer coefficients based on my code (as well as the figure I provided above) for implementing floating point filters:

int32_t sek1_x[3], sek1_y[3];
int32_t sek1_a[3] = { 2147483647, -2147483648, 1538959525 };
int32_t sek1_b[3] = { 2147483647, 2147483647, 2147483647 };

void lp_filter(int32_t *x, int32_t *y, int32_t *aa, int32_t *bb, int32_t g, int32_t in) {
    x[2] = x[1];
    x[1] = x[0];
    x[0] = ((int64_t)in * g) >> 31; // Scaling by shifting

    y[2] = y[1];
    y[1] = y[0];
    
    // Perform the filtering with scaling
    y[0] = (((int64_t)x[0] * bb[0] + (int64_t)x[1] * bb[1] + (int64_t)x[2] * bb[2] - (int64_t)aa[1] * y[1] - (int64_t)aa[2] * y[2])*2147483647) >> 31;
}

// Calling the filter function:
lp_filter(sek1_x, sek1_y, sek1_a, sek1_b, 25464180, adc_res);

I'm not sure that this will work, since I don't understand what part of the data is needed for, namely these lines:

const int DL[MWSPT_NSEC] = { 1,3,1 };
const int32_T DEN[MWSPT_NSEC][3] = {
  {
    2147483647,           0,           0 <- this line
  },
  {
    2147483647, -2147483648,  1538959525 
  },
  {
    2147483647,           0,           0 <- this line
  }
};

I also don’t understand how I can properly monitor overflow. I did shifts to try to avoid overflow, but I'm not sure if it's correct... Does anyone have any ideas on how to correctly use integer coefficients generated by matlab to implement digital filters?


I wrote an example from the information that I understood, but I made a mistake somewhere, but I don’t understand where... Source code with float filter implementation (6th order low pass filter and 2nd order high pass filter):

float   lp500_sek1_x[3],lp500_sek1_y[3],lp500_sek1_a[2]={-1.833933353,   0.8860449791},
        lp500_sek2_x[3],lp500_sek2_y[3],lp500_sek2_a[2]={-1.669203162,   0.7166338563},
        lp500_sek3_x[3],lp500_sek3_y[3],lp500_sek3_a[2]={-1.586906791,   0.6319990754},
        hp5_sek1_x[3],hp5_sek1_y[3],hp5_sek1_a[2]={-1.998889327,   0.9988899231};

//Calling filters
lp_filter(lp500_sek1_x,lp500_sek1_y,lp500_sek1_a,0.01302789338,((float)res_adc));
lp_filter(lp500_sek2_x,lp500_sek2_y,lp500_sek2_a,0.01185768284,lp500_sek1_y[0]);
lp_filter(lp500_sek3_x,lp500_sek3_y,lp500_sek3_a,0.01127306558,lp500_sek2_y[0]);
hp_filter(hp5_sek1_x,hp5_sek1_y,hp5_sek1_a,0.9994447827,lp500_sek3_y[0]);

//Filter functions
void lp_filter (float *x ,float *y,float *aa ,float g ,float vh)
{
    x[2]=x[1];
    x[1]=x[0];
    x[0]=vh*g;
    y[2]=y[1];
    y[1]=y[0];
    y[0]=(x[0] + x[1]+x[1] + x[2] - aa[0]*y[1] - aa[1]*y[2]); //b0=1 b1=2 b2=1
}
void hp_filter (float *x ,float *y,float *aa ,float g ,float vh)
{
    x[2]=x[1];
    x[1]=x[0];
    x[0]=vh*g;
    y[2]=y[1];
    y[1]=y[0];
    y[0]=(x[0] - (x[1]+x[1]) + x[2] - aa[0]*y[1] - aa[1]*y[2]); //b0=1 b1=-2 b2=1
}

Source code with an integer implementation of the same filters (6th order low-pass filter and 2nd order high-pass filter):

int32_t ilp500_sek1_a[2]={-1969170942, 951383551},
        ilp500_sek2_a[2]={-1792293247, 769479743},
        ilp500_sek3_a[2]={-1703928191, 678603839},
        ihp5_sek1_a[2]={-2146291070,  1072549887},
        ilp6_sek1_a[2]={-2145098367, 1071359231};

int64_t ilp500_sek1_x[3],ilp500_sek1_y[3],
        ilp500_sek2_x[3],ilp500_sek2_y[3],
        ilp500_sek3_x[3],ilp500_sek3_y[3],
        ihp5_sek1_x[3],ihp5_sek1_y[3],
        ilp6_sek1_x[3],ilp6_sek1_y[3];

//Calling filters
lp_filter_int(ilp500_sek1_x,ilp500_sek1_y,ilp500_sek1_a,13988593,((int64_t)res_adc));
lp_filter_int(ilp500_sek2_x,ilp500_sek2_y,ilp500_sek2_a,12732089,ilp500_sek1_y[0]);
lp_filter_int(ilp500_sek3_x,ilp500_sek3_y,ilp500_sek3_a,12104361,ilp500_sek2_y[0]);
hp_filter_int(ihp5_sek1_x,ihp5_sek1_y,ihp5_sek1_a,1073145662,ilp500_sek3_y[0]);

//Filter functions
void lp_filter_int (int64_t *x ,int64_t *y,int32_t *aa ,int32_t g ,int64_t vh)
{
    x[2]=x[1];
    x[1]=x[0];
    x[0]=(vh*(int64_t)g)>>30;
    y[2]=y[1];
    y[1]=y[0];
    y[0]=((x[0]*(int64_t)0x3FFFFFFF) + (x[1]*(int64_t)0x7FFFFFFF) + (x[2]*(int64_t)0x3FFFFFFF) - ((int64_t)aa[0]*y[1]) - ((int64_t)aa[1]*y[2]))>>30; //b0=1 b1=2 b2=1
}
void hp_filter_int (int64_t *x ,int64_t *y,int32_t *aa ,int32_t g ,int64_t vh)
{
    x[2]=x[1];
    x[1]=x[0];
    x[0]=(vh*(int64_t)g)>>30;
    y[2]=y[1];
    y[1]=y[0];
    y[0]=((x[0]*(int64_t)0x3FFFFFFF) - (x[1]*(int64_t)0x7FFFFFFF) + (x[2]*(int64_t)0x3FFFFFFF) - ((int64_t)aa[0]*y[1]) - ((int64_t)aa[1]*y[2]))>>30; //b0=1 b1=-2 b2=1
}

After the integer low-pass filter, I see a similar signal as after the float low-pass filter, but after the high-pass filter everything breaks down... I can’t understand what I’m doing wrong?

$\endgroup$
7
  • $\begingroup$ Fixed point processing is complicated. Most people use the so-called "Q" notation to track scale factors and avoid overflows. See: en.wikipedia.org/wiki/Q_(number_format) . You need to track the scale for EVERY state variable and intermediate result. Stage ordering and pole zero pairing can make a significant difference for SNR. Stay away from Direct Form II and Transposed Form I . Make sure you implement proper rounding. $\endgroup$
    – Hilmar
    Nov 24, 2023 at 12:58
  • 1
    $\begingroup$ Looking a little closer, I can see one problem. You're using Q1.31 and you should be using Q2.30 instead. Your $a_1$ coefficient (which gets scaled to be in aa[1]) can range from -2 to +2 in a biquad. So your integer value of your coefficients should be scaled up by $2^{30}$ and not by $2^{31}$. Do you understand this? $\endgroup$ Nov 26, 2023 at 5:54
  • $\begingroup$ My goodness, I knew that MATLAB had a C code generator, but I didn't know the quality of the code was so poor. They could do far better. $\endgroup$ Nov 26, 2023 at 6:00
  • $\begingroup$ @robert bristow-johnson As you say, the value of my float coefficients, which was generated by mallab, is in the range -2 to +2. If I want to convert these coefficients to integer, then I allocate 1 bit for the sign, and float=2 becomes equal to 2^31, which will give me coverage of the entire range, right? That is, for my example, the float array a={1,-1.669203162,0.7166a338563} will become the integer array a={1073741824,-1792293247,769479744}. $\endgroup$
    – red15530
    Nov 26, 2023 at 12:44
  • $\begingroup$ @robert bristow-johnson But there are two problems that I can't seem to figure out: 1. How do I control overflow? If my range is -2..2, overflows will happen all the time... 2. What should I do with integer coefficients for example b={1,2,1}, I will be forced to convert them to integer and do operations that I can avoid in float for speed (for example, multiplying by 1). $\endgroup$
    – red15530
    Nov 26, 2023 at 12:45

1 Answer 1

1
$\begingroup$

There is a lotta detail in the question.

I have done a hella lotta fixed-point DSP in my day. Both with a fixed-point DSP (56K) and with an integer CPU such as the 68K but also in generic C code. I posted an answer about specific C code that does a biquad in fixed point and uses noise-shaping at the sole quantization point. You should look at it.

Here is the main principle. You saw a comment refer to the "Q" notation. E.g. Q2.14 means a 16-bit fixed-point number with 2 bits to the left of the binary point and 14 to the right. What this really means is that the 16-bit integer has an integer value (this is the value that the C program sees) that is scaled up by a factor of $2^{14}$. This doesn't really matter for the signal; whatever scaling of the signal in the input is the same scaling for the output. But it does matter for the coefficients. Every coefficient is represented as an integer that is $2^{14}$ times larger than the actual numerical value of the coefficient.

So, when you integer multiply that coefficient times another value (like a signal or a state), the result will be $2^{14}$ times bigger than it should be. So you'll have to shift it 14 bits to the right to scale it with $2^{-14}$ and undo the scaling of the coefficient necessary to make it an integer. This shifting right causes 14 bits to fall offa the edge, which is quantization error. Noise shaping can help minimize the consequences of that quantization.

With 32 bits, the quantization error will be much smaller in effect.

So you have to deal with the range of your coefficients, so you know how much to scale them. For a biquad, we know the denominator coefficients will have absolute value of less than 2. So we can represent those coefficients well as Q2.30 fixed-point values. So using int32_t as your coefficients and int64_t as your accumulator, then after multiplying, you will have to shift the 64-bit value 30 bits to the right.

Besides noise-shaping (in this specific case, "fraction saving" is the noise shaping), there are some other useful tricks in fixed-point DSP. Another trick is to use the Direct Form 1 and do not use Direct Form II. The reasons why are two-fold: First the DF1 has only one point of quantization. And second, the DF1 puts the effect of the zeros before the poles, so in a high-Q filter, there is attenuation (from a zero) that will occur before there is gain (from the pole). This means that unnecessary internal overflow is avoided. And keeping a 64-bit accumulator and adding up all 5 values before quantization allows you to quantize after the addition.

$\endgroup$
24
  • $\begingroup$ Thank you very much for such a detailed explanation. I just don't understand some points. Why do I need to select the number form as Q2.30? I allocate one bit for the sign, and the remaining 31 I can use for my range of numbers. For example, if I want to express the numbers from -2 to 2. I can use like 2^31 = 2. $\endgroup$
    – red15530
    Nov 26, 2023 at 12:57
  • $\begingroup$ I also don’t understand what the right thing to do is: convert my float coefficients myself or have Matlab do this work? For example, if my original float array was a={1,-1.669203162,0.7166338563}, then the integer array would become a={1073741824,-1792293247,769479744}. But Matlab gives me a={2147483647,-2147483648, 1538959525}. $\endgroup$
    – red15530
    Nov 26, 2023 at 13:01
  • $\begingroup$ //"if I want to express the numbers from -2 to 2. I can use like 2^31 = 2."// - - - - - - really? You better look closely into that. $\endgroup$ Nov 26, 2023 at 13:01
  • $\begingroup$ If I understand you correctly, and I want to use the Q2.30 format, then I give the two most significant bits for the integer part and the sign, and the remaining 30 for the fractional part. So? My float array a={1,-1.669203162,0.7166338563}, would it be a={0x4000 0000,0xEAD4 397F,0x2DDD 543F}? $\endgroup$
    – red15530
    Nov 26, 2023 at 13:26
  • $\begingroup$ I think so. It looks correct. I am not converting the hex, but it looks right. Now, when you shift right, instead of shifting right 31 bits (>>31), you will be shifting right by 30 bits. $\endgroup$ Nov 26, 2023 at 13:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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