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: 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?
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$