Exactly I get a Zeitronix SDL file with Lambda values and in the Zeitronix software I can see calculated values and in the file were
only analog voltages wich dos not correspond to the voltage of the Lambda probe.
So I picked up a few points and needed to calculate a linear equation for this points.
First I searched in the net and found:
https://www.dcode.fr/function-equation-finder
But I thougt it would be nice to have this as PB code.
I used the calulation from:
https://blog.mbedded.ninja/mathematics/ ... e-fitting/
and borrowed the inverse matrix calculation from wilbert:
https://www.purebasic.fr/english/viewto ... 04#p523704
Combined, it results in LinearCurveFinder.pbi
Code: Select all
;
; General calulation of a linear curve fitting
; https://blog.mbedded.ninja/mathematics/curve-fitting/linear-curve-fitting/
;
; Calculation of an inverse matrix by wilbert
; https://www.purebasic.fr/english/viewtopic.php?p=523704#p523704
;
; Topic
; https://www.purebasic.fr/english/viewtopic.php?p=586282#p586282
;
CompilerIf #PB_Compiler_IsMainFile
EnableExplicit
CompilerEndIf
Procedure.i InverseMatrixDouble(*Matrix, *InvMatrix, n) ; n x n matrix inverse
Protected.l n1=n-1,p,k,t,v
Protected.d r,m,*s.Double,*d.Double
Protected Dim Mat.d(n1,n1)
Protected Dim Piv.u(n1)
CopyMemory(*Matrix, @Mat(), n*n*SizeOf(Double)) ; << copy *Matrix to Mat() >>
For t=0 To n1 : Piv(t)=t : Next ; << init Piv() >>
For p=0 To n1 : m=0
For t=p To n1 ; << find highest absolute value >>
r=Abs(Mat(p,Piv(t))) : If r>m : m=r : v=t : EndIf
Next
If m=0 : ProcedureReturn #False : EndIf ; << exit if all 0 >>
Swap Piv(p),Piv(v) : k=Piv(p) : r=1/Mat(p,k) : Mat(p,k)=1 ; << get k and reciprocal for Mat(p,k) >>
*d=@Mat(p,0) : For t=0 To n1 : *d\d*r : *d+8 : Next ; << multiply pivot row by reciprocal >>
*d=@Mat()
For t=0 To n1
If t<>p
*s=@Mat(p,0) : m=-Mat(t,k) : Mat(t,k)=0
For v=0 To n1 : *d\d+*s\d*m : *s+8 : *d+8 : Next ; << calculate new Mat() values >>
Else
*d+n<<3
EndIf
Next
Next
For t=0 To n1 : *d=*InvMatrix+Piv(t)*n<<3
For v=0 To n1 : *d\d=Mat(t,Piv(v)) : *d+8 : Next ; << set *InvMatrix >>
Next
ProcedureReturn #True
EndProcedure
Procedure.i LinearCurveFinder(List PointList.Point(), *a.double, *b.Double)
Protected.i Result
Protected.d A11, A12, A21, A22, B11, B21
Protected Dim DoubleArray.d(1, 1)
ForEach PointList()
A11 + PointList()\x * PointList()\x
A12 + PointList()\x
B11 + PointList()\x * PointList()\y
B21 + PointList()\y
Next
A21 = A12
A22 = ListSize(PointList())
; Debug "A11: " + Str(A11)
; Debug "A12: " + Str(A12)
; Debug "B11: " + Str(B11)
; Debug "B21: " + Str(B21)
; Debug "A22: " + Str(A22)
; [ A11 A12 ] [ a ] = [ B11 ]
; [ A21 A22 ] [ b ] = [ B21 ]
DoubleArray(0, 0) = A11
DoubleArray(0, 1) = A21
DoubleArray(1, 0) = A12
DoubleArray(1, 1) = A22
Result = InverseMatrixDouble(@DoubleArray(), @DoubleArray(), 2)
If Result
*a\d = DoubleArray(0, 0) * B11 + DoubleArray(0, 1) * B21
*b\d = DoubleArray(1, 0) * B11 + DoubleArray(1, 1) * B21
EndIf
ProcedureReturn Result
EndProcedure
CompilerIf #PB_Compiler_IsMainFile
#Width = 800
#Height = 600
Define.d a, b
NewList PointList.point()
Define.d MaxX, MaxY, MinX, MinY
Define.d XFactor, YFactor, Factor
Define.i FactorIsLargerThanOne, x, Event
Define Equation$
CompilerIf #True
AddElement(PointList())
PointList()\x = 2
PointList()\y = 1
AddElement(PointList())
PointList()\x = 3
PointList()\y = 7
AddElement(PointList())
PointList()\x = 6
PointList()\y = 5
AddElement(PointList())
PointList()\x = 7
PointList()\y = 9
CompilerElse
AddElement(PointList())
PointList()\x = 420
PointList()\y = 23
AddElement(PointList())
PointList()\x = 430
PointList()\y = 25
AddElement(PointList())
PointList()\x = 450
PointList()\y = 30
AddElement(PointList())
PointList()\x = 460
PointList()\y = 32
AddElement(PointList())
PointList()\x = 520
PointList()\y = 45
AddElement(PointList())
PointList()\x = 570
PointList()\y = 55
CompilerEndIf
If LinearCurveFinder(PointList(), @a, @b)
;Debug "y = " + StrD(a, 6) + "x + " + StrD(b, 6)
ForEach PointList()
If MaxX < PointList()\x : MaxX = PointList()\x : EndIf
If MaxY < PointList()\y : MaxY = PointList()\y : EndIf
If MinX > PointList()\x : MinX = PointList()\x : EndIf
If MinY > PointList()\y : MinY = PointList()\y : EndIf
Next
;Debug StrD(MinX, 0) + " < X < " + StrD(MaxX, 0)
;Debug StrD(MinY, 0) + " < Y < " + StrD(MaxY, 0)
Equation$ = "y = " + StrD(a, 6) + " x "
If b >= 0
Equation$ + "+"
Else
Equation$ + "-"
EndIf
Equation$ + " " + StrD(Abs(b), 6)
OpenWindow(0, 0, 0, #Width, #Height, "Line Curve Finder: " + Equation$, #PB_Window_MinimizeGadget|#PB_Window_ScreenCentered)
CanvasGadget(0, 0, 0, #Width, #Height)
XFactor = MaxX / #Width
If XFactor > 1 : FactorIsLargerThanOne = #True : EndIf
YFactor = MaxY / #Height
If XFactor > 1 : FactorIsLargerThanOne = #True : EndIf
If FactorIsLargerThanOne
If XFactor > YFactor
Factor = XFactor
Else
Factor = YFactor
EndIf
Else
If XFactor > YFactor
Factor = 1 / XFactor
Else
Factor = 1 / YFactor
EndIf
EndIf
;Debug "Factor: " + StrD(Factor)
If StartVectorDrawing(CanvasVectorOutput(0))
FlipCoordinatesY(#Height / 2)
ForEach PointList()
AddPathCircle(PointList()\x * Factor, PointList()\y * Factor, 4)
ClosePath()
Next
VectorSourceColor(RGBA(255, 0, 0, 255))
StrokePath(10)
MovePathCursor(0, b * Factor)
AddPathLine(MaxX * Factor, (a * MaxX + b) * Factor)
ClosePath()
VectorSourceColor(RGBA(0, 0, 255, 255))
StrokePath(4)
StopVectorDrawing()
EndIf
Repeat
Event = WaitWindowEvent()
Until Event = #PB_Event_CloseWindow
Else
Debug "Was not able to build the Inverse matrix"
EndIf
CompilerEndIf