Linear Regression Based Real-Time Filtering

This paper introduces real time filtering method based on linear least squares fitted line. Method can be used in case that a filtered signal is linear. This constraint narrows a band of potential applications. Advantage over Kalman filter is that it is computationally less expensive. The paper further deals with application of introduced method on filtering data used to evaluate a position of engraved material with respect to engraving machine. The filter was implemented to the CNC engraving machine control system. Experiments showing its performance are included.


Introduction
The problem of real-time signal filtering is addressed extensively in technical society.Variety of filters has been created.As an illustration, we can mention simple averaging or median filter; highpass, lowpass, bandpass filters [1] or widely used Kalman filter [2], [3].The superiority of the latter comes out of fact that it relies not only on measured data, but also on prior knowledge of a measured signal in the form of a mathematical model.Moreover, the model is iteratively updated based on previous measurements.
In this paper, we introduce filtering method, which is also using prior knowledge of the measured signal.We will show that under certain conditions linear regression can be used to filter a signal with known mathematical description.The idea of filtering lies in fitting a curve to already measured data, given the mathematical description of a curve.This basically means that if we know what kind of signal we expect we can interpolate a measured data by a curve of expected shape in real-time, thus filtering the noised data.
In the theoretical part of this paper, we describe mathematical fundamentals of linear regression (LR) filter for a line fit.
The practical part of the paper is focused on implementation of LR filter in a CNC engraving machine.For a correct functionality, a constant distance between engraving tool and the material has to be held.To accomplish a constant distance, proximity measurements are taken and the tool's position in Z axis is controlled accordingly.The proximity sensor data are corrupted by noise and need to be filtered.To filter the data and to estimate the position of material surface with reference to the machine's frame, we use LR filter.

Least Squares Filter
Least squares (LS) error fitting is widely used method in statistics and other fields of math.Usually it is used to fit a curve to already measured set of data.
Given a set of n points (data pairs) (x i ,y i ), i = 1, . . ., n, where x i is an independent variable and y i is dependent variable, we can fit a line y = ax + b, using following equations: This is a common least squares line fitting problem solution and complete derivation of above equations can be found in [4].Theoretically, it can be used to filter any mathematically described signal which is a linear combination of functions in form: where β j are coefficients and φ j are functions of independent variable xi.The signal can also be nonlinear, but there is no closed-form solution for such LR problem.Instead, numerical algorithms must be used [5].Computing LR for more complex functions requires more sophisticated algorithms such as LU (lower-upper) decomposition, or singular value decomposition, which are computationally more expensive and real time filtering would require more powerful environment.
If we know that measured signal is linear, we can use this method for real time signal filtering.With each measured point, we can compute a new fitted line.The more points there are the more accurate fit we get.Pseudo code of LS real time filter would look as follows: sumx = 0; sumy = 0; sumxx = 0; sumxy = 0; n=0 //initialize variables while (true) { x = determine_x(); y = measure_y(); n++; //iterate number of points sumx += x; //update sum of x 's sumy += y; //update sum of y 's sumxy += x * y; //update sum of (x * y) sumxx += x * x; //update sum of (x * x) a = (n * sumxy -sumx * sumy) / (n * sumxxsumx * sumx); //calculate parameter a b = (sumy -a * sumx) / n; //calculate paramenter b } Advantage of this code is, that it is very fast and consumes a small amount of memory.There is no need to remember all measured points, only sums of x and y need to be stored.Figure 1 shows filtering of noisy signal with LS filter.In Fig. 2 averaging filter was used on the same dataset.At the beginning, performance of LS filter and averaging filter is comparable.After having more than 20 data points (x = 0) measured, LS filter clearly outperforms averaging filter.
We can compare the complexity of LR filter to Kalman filter, since it also relies on a model of measured signal.The Kalman filter model assumes the true state at time k is evolved from the state at (k − 1) according to: where F k is the state transition model which is applied to the previous state x k−1 ; B k is the control input model which is applied to the control vector u k ; w k is the process noise which is assumed to be drawn from a zero mean multivarate normal distribution with covariance Q k [2].At time k an observation (or measurement) z k of the true state x k is made according to: where H k is the observation model which maps the true state space into the observed space and v k is the observation noise which is assumed to be zero mean Gaussian white noise with covariance R k [2].
The filtered value x k|k of measured signal in time k is computed in two stages: predict and update.Predict: Update: Practical implementation of Kalman filter requires a model of the system in the form of matrices F k , B k , H k , and process noise and sensor noise covariances Q k and R k .The more precise model and estimate of covariances are, the better the performance of filter is achieved.From a computing point of view, in the general case, Kalman filter requires multiplications, additions and inversion of matrices.Computationally most expensive is the inversion of matrix in Eq. ( 8).

Application of LS Filter
In our application, LS filter is used to filter a signal from a proximity sensor which is measuring a distance of material that is engraved by CNC machine.The machine was built for engraving planar solid surfaces, such as stone blocks.The engraved picture is not 3D, but only planar (similar to drawing).Engraving is performed by hammering on the material with pointed diamond.This would theoretically require X − Y plotter which has a working plane in parallel with material.This parallelism must be rather accurate, so that a constant distance between the tool and material is held.In our case, required distance is 1 ± 0, 5 mm.Bedding the material with such accuracy is often a difficult task.For this reason, a degree of freedom in Z-axis was added to the machine.The only purpose of Z-axis movement is to compensate for non-parallelism of material surface with respect to working plane of the machine.A deflection from parallelism is evaluated by measuring the distance between the diamond tool tip and the material surface in certain points.Measurements are performed by rather unconventional proximity sensor which uses touching of a tool tip with material to evaluate the distance.A detailed description of the sensor can be found in [6].From the nature of touch sensor principle, it is obvious that the bigger the distance, the longer it takes to measure it.In our machine small distances up to 3 mm can be measured in tens of milliseconds, but a 15 mm distance can take up to 3 seconds (!) to measure.This constraint comes out of fact, that small distances are measured only by moving the tool tip which has low mass.In contrast, for longer distances, whole engraving tool with considerably bigger mass has to be moved.Taking in account above constraints, sensor allows to measure distance with each and every stroke of hammering tool while engraving (the engraving distance is 1 mm).In addition, the measurement is taken during the time the tool is performing its stroke, so it does not introduce any additional time delay.
On the other hand, initial evaluating of material surface plane might take longer, because every point in which we measure distance might take 3 seconds to measure.That is why amount of initial measurement points is limited.

Estimation of Material Surface Plane Position
To evaluate the parallelism deflection of material surface and machine working plane, we measure distance in a few representative points and fit a plane on these points.Based on this equation, a Z-axis would be driven, so that a constant distance between the tool and material is kept.To define a plane, at least 3 points are necessary.To make sure, that engraved picture lies completely on a material surface, we picked from 4 to 6 representative measurement points depending on outer shape of the picture.Since the sensor has some error, measured points do not lie exactly on a plane.So we came to over-determined task of fitting a plane to 4 ÷ 6 points.To do this, we used a multiple linear least squares regression [7], [8], [9], [10].The plane equation is: z can be calculated as: Since there is an additional sensor error in z direction we can write: where As we have more measured points, we can write the Eq.(13) in matrix form: where n is the number of measured points and k is a vector of unknown parameters.
The square residual error E 2 is calculated by rearranging and then squaring E as is shown below: At the minimum square residual error, the derivatives of E 2 with respect to the unknown k are equal to zero: This gives the equation: By solving the equation we get the vector of parameters The equation of fitted plane ρ is then: Estimation of material surface position is depicted in Fig. 3.The origin of coordinate system is placed into the tip of the tool.After measuring five representative points (pink) a fitted plane ρ is calculated.

Z-Axis Position Control
The computed fitted plane equation is fetched to the machine control system in the form of coordinates of image corners.
Matrix C basically represents a rectangle that belongs to plane ρ.Each row of C is a corner of engraved image.We denote them C1, C2, C3 and C4 where C1 = (x 1 , y 1 , z 1 ) and so on.Subsequently, two lines l and r, identical with left and right side of the rectangle are computed using parametric equations: where parameter t ∈< 0, 1 >.Finally, equation of the middle line m, intersecting left and right line, perpendicular to both lines is computed: where parameter s ∈< 0, 1 >.The engraved image is the result of thousands of strokes of the tool.When projected to X − Y plane (Fig. 4), the tool trajectory starts at the top of the image and continues row by row, until it reaches the bottom of the image.We calculate the Z-axis tool position in any 2D point p with coordinates (x, y) as follows: 1. Calculate points L L and L R letting t = y/img height.
2. Calculate point L M letting s = x/img width.
Z coordinate of point L M is a desired Z-axis tool position in point p.Since coordinates (x, y) are known, it is sufficient to calculate only z-coordinate of vectors L L , L R and L M .

On the Fly Parameters Tuning
A few experiments have shown that driving a Z-axis tool position based only on the fitted plane ρ is not accurate enough.Every measured point involves some error.Fitted plane minimizes squared sum of errors.
In adverse cases it might happen that summed error of fitted plane is much smaller than measurement error, which concludes to rather big error in coincidence of a real and fitted plane.In worst cases, the Z axis error was bigger than 1 mm.However, the fitted plane gives a good estimate of where the real material surface is.
Computing z-coordinate of the tool is based on lines l, r and m.By adjusting the parameters of these lines, more accurate Z-axis control could be achieved.Refinement of parameters of these lines is accomplished by carrying out further measurements of surface distance during engraving.Proximity sensor allows measuring distance with each stroke of hammering tool.Regular image consists of tens to hundreds of strokes in each row (Fig. 4).By fitting a line m f it to points measured in a row, we get more accurate estimate of surface in given row (Fig. 5).In order to handle cases with a few measurement points (< 20), we also include marginal points (intersections of lines m, l and m, r) to set of points fitted by m f it .The equation of m f it is calculated at the end of each row.Based on this equation, we calculate a z-coordinate of new point L L tuned (0, y curr , z) and point L (R tuned (img width, y curr , z), where y curr is y-coordinate of the current row.Points L L tuned and L R tuned serve for tuning parameters of lines l and m, respectively.This is done by calculating new lines l f it and r f it .Line l f it is fitted on set {C1, C3, L L tuned i }, line r f it is a fitted on set {C2, C4, L R tuned i }, where i = 0, . . ., y curr (Fig. 6).Finally, we update a zcoordinate of corners C1, C2, C3 and C4 by calculating it using starting and ending points of l f it and r f it .In the next image row, new values of corners are used to calculate lines l, r and m.In Fig. 8, there is a X − Z projection of engraved image from Fig. 4. The measurement shows that the distance is kept on 1 mm with an error < 0, 5 mm which comprises with initial requirements of engraving machine.

Conclusion
In this paper the real time filtering method based on linear regression has been introduced.The filter is based on the idea of fitting a curve on a set of already measured data, supposing that a mathematical model of the signal is available.The experimental measurements were focused on practical application of filter in CNC machine, to estimate a position of machined material, where a filtered signal has the shape of a line.Therefore, computing of fit is simple and fast.LR filter in comparison to Kalman filter is less expensive from a computing point of view because in the general case Kalman filter in practical implementation requires multiplications, additions and inversion of matrices.
In this application, a main goal was to keep a constant distance of engraving tool from engraved material.Numerous portraits have been engraved using proposed system, while measuring and saving the tool position data, which confirm a correct functionality of the system, i.e. that the system is capable of keeping the tool in required distance.The control system of the machine is based on an ATmega162, 8-bit, fixed-point microcontroller.
The advantage of our method is in computational simplicity since it requires only a few summations, multiplications and divisions of floating point variables so it is well suited for a real time use in computationally non-powerful environments e.g. with 8-bit microcontrollers.

Fig. 3 :
Fig. 3: CNC machine mb2300; Red -origin of coordinate system is placed on the tip of diamond tool.Yellow-fitted plane ρ is calculated out of five measured points (pink).

Figure 4 toFig. 8 :
Figure4to Fig.7were created by plotting a measured data from real engraved picture.In Fig.7, fitted plane ρ (yellow) computed out of five measured points (pink crosses) is depicted.Blue plane is the final update of ρ at the end of engraving.As can be seen, two planes are slightly displaced.That confirms the fact,