diff options
Diffstat (limited to 'werteunit.pas')
| -rw-r--r-- | werteunit.pas | 179 |
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; |