Implementation notes of least square fitting conic in stm32

        The sensor shall be calibrated before use , Curve fitting must be carried out during calibration . It's easy to deal with it by computer , But in practice, the microcontroller is used to bid timing , Only general algebraic operations can be performed , No matrix operation , It is very inconvenient to handle . The algebraic formula for calculating the undetermined coefficients of quadratic polynomial curve fitting is derived by the least square method , It is very convenient to use these formulas to process data .

        There is a set of measured data (x i , yi)i =1 , 2 , … , n , The fitting function is

                                                                                  

• One 、MATLAB Verification in China

Assume that the standard input and output curve of a sensor is as shown in Figure 1( Red curve ) Shown , The input and output curves of the sensor deviate due to the environment and external interference , Pictured 1( Blue curve ). If you want to get the correct data from the sensor, you need a correction curve , Correct the input curve with error to a red standard curve . First we have to confirm the error curve , Pictured 2 Shown . Red is the smoothed curve .

chart 1： Standard curve and actual curve

chart 2： Error curve

 

MATLAB The correction process is briefly described as follows ：

1. Import standard curves into MATLAB
2. Import the error curve collected on site into MATLAB
3. Take the standard curve and the curve with error 7 Data
4. Calculation error ( Output of fitting curve )
5. Set error curve by least square method
6. F( standard output )  = F( Field input )-F( error )

The fitted curve is shown in the figure 3 Shown ：

Fit the curve

 

MATLAB The code is as follows ( Direct operation )：

%============== Define the standard ====================%
a = [310.0 305.9 302.6 300.2 298.1 296.0 294.1 292.4 290.8 288.8 287.4 286.1 284.8 283.7 282.7 281.9 281.0 280.2 279.4 278.8 278.2 277.6 277.0 276.4 275.9 275.5 275.0 274.3 273.8 273.4];
a = a/10.0;
%============== Define the actual input ====================%
b = [308.7 304.3 301.0 298.4 296.2 294.1 291.9 290.0 288.2 286.6 285.2 283.7 282.5 281.5 280.5 279.5 278.7 277.9 277.2 276.5 275.9 275.2 274.7 274.0 273.5 273.0 272.5 272.0 271.5 271.1];
b = b/10.0;
bb = [b(1),b(5),b(10),b(15),b(20),b(25),b(30)];% Data to be fitted , Take one of the input 7 Data 
%============= Define standard output ( Is the error )=================%
out = [a(1)-b(1),a(5)-b(5),a(10)-b(10),a(15)-b(15),a(20)-b(20),a(25)-b(25),a(30)-b(30)];
%================= The following is the calculation parameter ===================%
b1 = sum(bb);
b2 = sum(bb.^2);
b3 = sum(bb.^3);
b4 = sum(bb.^4);

c1 = sum(out);
c2 = sum(out.*bb);
c3 = sum(out.*bb.^2);

n=7;% The length of fitting data is  7 
k = n*b2*b4 + 2*b1*b2*b3 - n*b3^2 - b1^2*b4 - b2^3;%% Calculation K
a0 = ((b2*b4-b3^2)*c1 + (b2*b3-b1*b4)*c2 + (b1*b3-b2^2)*c3 )/k;
a1 = ((b2*b3-b1*b4)*c1 + (n*b4-b2^2)*c2 + (b1*b2-n*b3)*c3 )/k;
a2 = ((b1*b3-b2^2)*c1 + (b1*b2-n*b3)*c2 + (n*b2-b1^2)*c3 )/k;
%% Fit function %%%%%
syms funb x0
funb = a0 + a1*x0 + a2*x0*x0;%% function 
%================ Error curve and smoothed curve ============%
error = a-b;% Get the data error 
t = 1:1:30;% Define abscissa 
error_p = zeros(1,30);% Define an empty vector 
for i=2:30
    error_p(i) = error_p(i-1)*0.8 + error(i)*(1-0.8);% First order low-pass filter 
end
plot(t,error);% Draw an error curve 
hold on;
plot(t,error_p);% Draw a smoothing error curve 
figure;
ezplot(funb,[27,38]);% Draw a fit curve 

 

• Two 、 stay STM32 To realize

since STM32F4 After the launch of the series ,STM32 Already with DSP Arithmetic unit , Even if you need real-time curve fitting ,STM32 Can also be fully competent .

The calculation code is as follows ：

/*
 Show function 
 Error curve fitting 
 among  n  To fit the data length 
 among g_dData  Is the curve data to be fitted (x)
 among g_dOut  Is the output of the fitting curve (Y)
*/
void disposeHMI()
{
	double b1=0;
	double b2=0;
	double b3=0;
	double b4=0;
	double c1=0;
	double c2=0;
	double c3=0;
	double n=0;
	double k=0;
	unsigned char i;
	unsigned int j;
	if(g_cUartData == 0xfe)// If a measurement instruction is issued 
	{
		g_cUartData = 0;// Empty 
		Bi = 0;
		delay_ms(200);
		Bi = 1;
		//============================= The following is the fitting ==========================//
		//========================= The fitting data is  1 5 10 15 20 25 30 ====7 Number ======//
		{
			n=7;
			//================ Let's calculate BCk coefficient 
			for(i=0;i<7;i++)
			{
				b1 += g_dData[i];
			}
			for(i=0;i<7;i++)
			{
				b2 += pow(g_dData[i],2);
			}
			for(i=0;i<7;i++)
			{
				b3 += pow(g_dData[i],3);
			}
			for(i=0;i<7;i++)
			{
				b4 += pow(g_dData[i],4);
			}
			
			for(i=0;i<7;i++)
			{
				c1 += g_dOut[i];
			}
			for(i=0;i<7;i++)
			{
				c2 += g_dOut[i]*g_dData[i];
			}
			for(i=0;i<7;i++)
			{
				c3 += g_dOut[i]*pow(g_dData,2);
			}
			k = n*b2*b4 + 2*b1*b2*b3 - n*b3*b3 - b1*b1*b4 - b2*b2*b2;
			//============== Let's calculate the coefficient ============================//
			a0 = ((b2*b4-b3*b3)*c1 + (b2*b3-b1*b4)*c2 + (b1*b3-b2*b2)*c3 )/k;
			a1 = ((b2*b3-b1*b4)*c1 + (n*b4-b2*b2)*c2 + (b1*b2-n*b3)*c3 )/k;
			a2 = ((b1*b3-b2*b2)*c1 + (b1*b2-n*b3)*c2 + (n*b2-b1*b1)*c3 )/k;
		}
	}
}

3、 ... and 、 summary

　 It is convenient to calculate the undetermined coefficient by using the algebraic calculation formula of fitting the undetermined coefficient with the quadratic term curve derived in this paper , then , Then curve fitting , Its fitting accuracy is also very high , In particular, the advantages of calibration with single-chip microcomputer are more significant . The author's ability is limited , If there is any error, please prompt for correction .