Gaußfit eingebaut

1 files changed, 179 insertions, 0 deletions

diff --git a/werteunit.pas b/werteunit.pas
index e35b54a..a643a7b 100644
--- a/werteunit.pas
+++ b/werteunit.pas
@@ -62,6 +62,7 @@ type
     procedure integriereSingle(qu: pTLLWerteSingle; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung);
     procedure integriereDouble(qu: pTLLWerteDouble; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung);
     procedure integriereExtended(qu: pTLLWerteDouble; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung);
+    procedure gauszFit(amplituden,breiten,positionen,ueberlappe,hintergruende: pTLLWerteExtended; von,bis: longint; senkrecht: boolean; fensterBreite,maxBreite,maxVerschiebung: extended; positionsMitten: tExtendedArray);
   end;
   tLLWerteSingle = specialize tLLWerte<single>;
   tLLWerteDouble = specialize tLLWerte<double>;
@@ -1254,6 +1255,184 @@ begin
   end{of case};
 end;
 
+procedure tLLWerte.gauszFit(amplituden,breiten,positionen,ueberlappe,hintergruende: pTLLWerteExtended; von,bis: longint; senkrecht: boolean; fensterBreite,maxBreite,maxVerschiebung: extended; positionsMitten: tExtendedArray);
+var
+  i,j,ii,zaehl,                    // Zähler
+    qpSchritt,qsSchritt,           // Schrittlängen in der Quelle
+    zpSchritt,zdSchritt,           // Schrittlängen im Ziel
+    pLen,offset,                   // Datenlänge und Offset in der Quelle
+    wiMin,wiMax,                   // momentane Fensterränder (für den Zählindex) in der Quelle
+    verbesserung,                  // <0: schlechter, 1: besser, 2: viel besser
+    aParams,nParams:     longint;  // Index der alten und neuen Parameter
+  parameter:             array[0..1,boolean,0..3] of extended;
+    // dim0 nummeriert Parametersätze
+    // dim1: Ableitung (true) oder Werte (false)
+    // dim2: position, 1/breite, amplitude, hintergrund
+  fehlers:               array[0..1] of extended;
+  t0,t1,t2,t3,                     // Zwischenergebnisse
+    qdrSumm,                       // Quadratesumme im betrachteten Fenster
+    pMi,pMa,                       // Achsenenden in Datenrichtung
+    schrittFaktor:       extended;
+  ignBr,ignVersch:       boolean;  // Breite/Verschiebung am Anschlag!
+begin
+  if senkrecht then begin
+    qpSchritt:=params.xsteps; // Schritt in der Quelle parallel zur Fit-Richtung
+    qsSchritt:=1;             // Schritt in der Quelle senkrecht zur Fit-Richtung
+    zpSchritt:=params.xsteps; // Schritt im Ziel in positionsMitten-Richtung
+    zdSchritt:=1;             // Schritt im Ziel in Daten-Richtung (= senkrecht zur Fit-Richtung)
+    pMi:=params.tstart;
+    pMa:=params.tstop;
+    pLen:=params.tsiz;
+  end
+  else begin
+    qsSchritt:=params.xsteps;
+    qpSchritt:=1;
+    zdSchritt:=length(positionsMitten);
+    zpSchritt:=1;
+    pMi:=params.xstart;
+    pMa:=params.xstop;
+    pLen:=params.xsteps;
+  end;
+  maxBreite:=1/maxBreite;
+
+  for i:=von to bis do
+    for j:=0 to length(positionsMitten)-1 do begin
+      offset:=i*qsSchritt;
+      aParams:=0;
+
+      wiMin:=0;
+      wiMax:=pLen-1;
+      if fensterBreite>0 then begin
+        wiMin:=max(wiMin,round(positionsMitten[j] - fensterBreite / 2));
+        wiMax:=min(wiMax,round(positionsMitten[j] + fensterBreite / 2));
+      end;
+
+      qdrSumm:=0;
+      t0:=0;
+      t1:=0;
+      t2:=0;
+      schrittFaktor:=0;
+      for ii:=wiMin to wiMax do begin
+        t3:=werte[offset+qpSchritt*ii];
+        t3:=abs(t3);
+        qdrSumm:=qdrSumm + sqr(t3);
+        t0:=t0 + t3;
+        t1:=t1 + ii*t3;
+        t2:=t2 + sqr(ii)*t3;
+        if t3>schrittFaktor then schrittFaktor:=t3;
+      end;
+
+      t1:=t1/t0;
+      t2:=t2/t0;
+
+      parameter[aParams,false,0]:=t1;                     // Erwartungswert
+      parameter[aParams,false,2]:=schrittFaktor;          // Maximalwert
+      parameter[aParams,false,1]:=parameter[aParams,false,2]/t0*sqrt(pi);
+//      parameter[aParams,false,1]:=1/sqrt(2*(t2-sqr(t1))); // beinhalted Standardabweichung  <-- schlechter!!
+      parameter[aParams,false,3]:=0;
+
+      if (maxBreite>=0) and (parameter[aParams,false,1]<2*maxBreite) then
+        parameter[aParams,false,1]:=maxBreite*2;
+
+      nParams:=aParams;
+      {$DEFINE gauszFitBerechneWerte}
+      {$I gauszFit.inc} // Werte + Gradienten berechnen
+      {$UNDEF}
+
+      schrittFaktor:=1;
+      if maxVerschiebung>=0 then
+        schrittFaktor:=min(schrittFaktor,maxVerschiebung/max(abs(parameter[nParams,true,0]),1e-50));
+      if maxBreite>=0 then
+        schrittFaktor:=min(schrittFaktor,abs(parameter[nParams,false,2]-maxBreite)/max(abs(parameter[nParams,true,1]),1e-50));
+
+      zaehl:=0;
+      repeat
+        nParams:=(aParams+1) mod length(parameter);
+        for ii:=0 to length(parameter[aParams,false])-1 do
+          parameter[nParams,false,ii]:=
+            parameter[aParams,false,ii] - schrittFaktor * parameter[aParams,true,ii]; // * byte(ii<3);
+
+        ignVersch:=false;
+        if maxVerschiebung>0 then begin
+          if parameter[nParams,false,0]-positionsMitten[j] > maxVerschiebung then begin
+            parameter[nParams,false,0]:=positionsMitten[j] + maxVerschiebung;
+            ignVersch:=true;
+          end;
+          if positionsMitten[j]-parameter[nParams,false,0] > maxVerschiebung then begin
+            parameter[nParams,false,0]:=positionsMitten[j] - maxVerschiebung;
+            ignVersch:=true;
+          end;
+        end;
+
+        ignBr:=false;
+        if maxBreite>0 then
+          if parameter[nParams,false,1] < maxBreite then begin
+            parameter[nParams,false,1]:=maxBreite;
+            ignBr:=true;
+          end;
+
+        {$DEFINE gauszFitBerechneWerte}
+        {$I gauszFit.inc} // Werte + Gradienten berechnen
+        {$UNDEF}
+
+        if fehlers[aParams]>2*fehlers[nParams] then begin
+          verbesserung:=2;
+          schrittFaktor:=schrittFaktor * 2;
+          aParams:=nParams;
+        end
+        else if fehlers[aParams]>fehlers[nParams] then begin
+          verbesserung:=1;
+          schrittFaktor:=schrittFaktor * 1.3;
+          aParams:=nParams;
+        end
+        else begin
+          if verbesserung>0 then
+            verbesserung:=-1
+          else
+            dec(verbesserung);
+          schrittFaktor:=schrittFaktor * 0.5;
+        end;
+
+        if schrittFaktor<1e-50 then
+          fehler('Sehr kleiner Schrittfaktor ('+floattostr(schrittFaktor)+') beim Fitten!');
+
+        inc(zaehl);
+        {$IFDEF gauszFitStatus}
+        if (zaehl mod 10000)=0 then
+          gibAus(
+            floattostr(fehlers[aParams])+' '+
+            floattostr(qdrSumm)+' '+
+            inttostr(byte((verbesserung<-10) or (fehlers[aParams]*100<qdrSumm)))+' '+
+            inttostr(byte(ignVersch or (abs(schrittFaktor*parameter[aParams,true,0])<1)))+' '+
+            inttostr(byte(ignBr or (abs(schrittFaktor*parameter[aParams,true,1]/parameter[aParams,false,1])<0.01)))+' '+
+            inttostr(byte(abs(schrittFaktor*parameter[aParams,true,2]/parameter[aParams,false,2])<0.01))+' '+
+            inttostr(byte(abs(schrittFaktor*parameter[aParams,true,3]/max(abs(parameter[aParams,false,3]),1e-10))<0.01))+' '+
+            inttostr(verbesserung),
+            3);
+        {$ENDIF}
+
+      until ((zaehl>=100000) or (verbesserung<-10) or (fehlers[aParams]*10<qdrSumm)) and
+            (ignVersch or (abs(schrittFaktor*parameter[aParams,true,0])<1)) and
+            (ignBr or (abs(schrittFaktor*parameter[aParams,true,1]/parameter[aParams,false,1])<0.01)) and
+            (abs(schrittFaktor*parameter[aParams,true,2]/parameter[aParams,false,2])<0.01) and
+            (abs(schrittFaktor*parameter[aParams,true,3]/max(abs(parameter[aParams,false,3]),1e-10))<0.01);
+      {$IFDEF gauszFitStatus}
+      gibAus('',3);
+      {$ENDIF}
+
+      if assigned(positionen) then
+        positionen^.werte[j*zpSchritt + i*zdSchritt]:=pMi+(pMa-pMi)/(pLen-1)*parameter[aParams,false,0];
+      if assigned(breiten) then
+        breiten^.werte[j*zpSchritt + i*zdSchritt]:=(pMa-pMi)/(pLen-1)/parameter[aParams,false,1];
+      if assigned(amplituden) then
+        amplituden^.werte[j*zpSchritt + i*zdSchritt]:=parameter[aParams,false,2];
+      if assigned(hintergruende) then
+        hintergruende^.werte[j*zpSchritt + i*zdSchritt]:=parameter[aParams,false,3];
+      if assigned(ueberlappe) then
+        ueberlappe^.werte[j*zpSchritt + i*zdSchritt]:=fehlers[aParams]/qdrSumm;
+    end;
+end;
+
 // tWavelet ********************************************************************
 
 function tWavelet.setzeTyp(s: string): boolean;