diff options
Diffstat (limited to 'werteunit.pas')
| -rw-r--r-- | werteunit.pas | 209 |
1 files changed, 199 insertions, 10 deletions
diff --git a/werteunit.pas b/werteunit.pas index 4d67c4c..9a068a4 100644 --- a/werteunit.pas +++ b/werteunit.pas @@ -19,6 +19,7 @@ type Auch die korrekte Parallelisierung obliegt dem übergeordneten Programmteil. } private + procedure sortiereMaxima(var maxima: tIntPointArray); public werte: array of wgen; params: tExtrainfos; @@ -46,8 +47,8 @@ type procedure kopiereLOVerzerrt(original: pTLLWerteExtended; xmin,xmax,tmin,tmax: longint; verhHo,verhVe: extended); overload; destructor destroy; override; function liesDateien(dateien: tGenerischeInputDateiInfoArray): boolean; - function fft(senkrecht,invers: boolean; algo: tFFTAlgorithmus; fen: tFenster; out pvFehler: extended): boolean; overload; inline; - function fft(smin,smax: longint; senkrecht,invers: boolean; algo: tFFTAlgorithmus; fen: tFenster; out pvFehler: extended): boolean; overload; + function fft(senkrecht,invers: boolean; algo: tFFTAlgorithmus; fen: tFenster; hg: extended; out pvFehler: extended): boolean; overload; inline; + function fft(smin,smax: longint; senkrecht,invers: boolean; algo: tFFTAlgorithmus; fen: tFenster; hg: extended; out pvFehler: extended): boolean; overload; procedure spiegle; overload; procedure spiegle(tmin,tmax: longint); overload; procedure fft2dNachbearbeitungA(nb: tFFTDatenordnung); @@ -66,6 +67,8 @@ type 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); + function ermittleHintergrund: extended; + procedure filtereHoheFrequenzen(betraege: tLLWerte; xFak,yFak: extended); end; tLLWerteSingle = specialize tLLWerte<single>; tLLWerteDouble = specialize tLLWerte<double>; @@ -90,6 +93,92 @@ implementation uses systemunit; +// tLLWerte ******************************************************************** + +procedure tLLWerte.sortiereMaxima(var maxima: tIntPointArray); +var + mins,maxs: tExtendedArray; + pivot,wert,wertLi,wertRe: extended; + vons,biss: tLongintArray; + tmp: tIntPoint; + i,li,re,cnt: int64; +begin + setlength(vons,1); + vons[0]:=0; + setlength(biss,1); + biss[0]:=length(maxima)-1; + setlength(mins,1); + setlength(maxs,1); + mins[0]:=0; + maxs[0]:=0; + for i:=vons[0] to biss[0] do begin + wert:=werte[maxima[i]['x']+maxima[i]['y']*params.xsteps]; + if (i=vons[0]) or (wert>maxs[0]) then + maxs[0]:=wert; + if (i=vons[0]) or (wert<mins[0]) then + mins[0]:=wert; + end; + cnt:=1; + + while cnt>0 do begin + li:=vons[cnt-1]; + re:=biss[cnt-1]; + + if (li>=re) or (mins[cnt-1]=maxs[cnt-1]) then begin + dec(cnt); + continue; + end; + + if cnt>=length(vons) then begin + setlength(vons,cnt+100); + setlength(biss,cnt+100); + setlength(mins,cnt+100); + setlength(maxs,cnt+100); + end; + + pivot:=sqrt(maxs[cnt-1]*mins[cnt-1]);//(maxs[cnt-1]+mins[cnt-1])/2; + + mins[cnt]:=mins[cnt-1]; + biss[cnt]:=biss[cnt-1]; + maxs[cnt]:=mins[cnt]; + mins[cnt-1]:=maxs[cnt-1]; + + while li<=re do begin + wertLi:=werte[maxima[li]['x']+maxima[li]['y']*params.xsteps]; + if wertLi>=pivot then begin + if wertLi<mins[cnt-1] then + mins[cnt-1]:=wertLi; + inc(li); + continue; + end; + wertRe:=werte[maxima[re]['x']+maxima[re]['y']*params.xsteps]; + if wertRe<=pivot then begin + if wertRe>maxs[cnt] then + maxs[cnt]:=wertRe; + dec(re); + continue; + end; + if wertLi>maxs[cnt] then + maxs[cnt]:=wertLi; + if wertRe<mins[cnt-1] then + mins[cnt-1]:=wertRe; + tmp['x']:=maxima[re]['x']; + tmp['y']:=maxima[re]['y']; + maxima[re]['x']:=maxima[li]['x']; + maxima[re]['y']:=maxima[li]['y']; + maxima[li]['x']:=tmp['x']; + maxima[li]['y']:=tmp['y']; + inc(li); + dec(re); + end; + + vons[cnt]:=li; + biss[cnt-1]:=re; + + inc(cnt); + end; +end; + constructor tLLWerte.create(ps: tExtrainfos); begin inherited create; @@ -902,19 +991,19 @@ begin end; end; -function tLLWerte.fft(senkrecht,invers: boolean; algo: tFFTAlgorithmus; fen: tFenster; out pvFehler: extended): boolean; +function tLLWerte.fft(senkrecht,invers: boolean; algo: tFFTAlgorithmus; fen: tFenster; hg: extended; out pvFehler: extended): boolean; begin if senkrecht then - result:=fft(0,params.xsteps-1,senkrecht,invers,algo,fen,pvFehler) + result:=fft(0,params.xsteps-1,senkrecht,invers,algo,fen,hg,pvFehler) else - result:=fft(0,params.tsiz-1,senkrecht,invers,algo,fen,pvFehler); + result:=fft(0,params.tsiz-1,senkrecht,invers,algo,fen,hg,pvFehler); end; -function tLLWerte.fft(smin,smax: longint; senkrecht,invers: boolean; algo: tFFTAlgorithmus; fen: tFenster; out pvFehler: extended): boolean; +function tLLWerte.fft(smin,smax: longint; senkrecht,invers: boolean; algo: tFFTAlgorithmus; fen: tFenster; hg: extended; out pvFehler: extended): boolean; var i,j,pmax,pstep,sstep: longint; in0,out0: boolean; - vorher,nachher: extended; + vorher,nachher,offset: extended; begin result:=false; if senkrecht then begin @@ -929,12 +1018,15 @@ begin end; if assigned(fen) and fen.aktiv then begin + if invers then + gibAus('fft: Warnung, hier wird bei einer inversen FFT gefenstert - soll das so sein?',1); + offset:=byte(not invers)*hg; if length(fen.Werte)<>pmax+1 then fen.berechneWerte(pmax+1); for i:=smin to smax do // Werte fenstern for j:=0 to pmax do werte[i*sstep+j*pstep]:= - werte[i*sstep+j*pstep]*fen.Werte[j]; + (werte[i*sstep+j*pstep]-offset)*fen.Werte[j]; end; vorher:=0; @@ -972,6 +1064,12 @@ begin if (nachher=0) and (vorher=0) then pvFehler:=0 else pvFehler:=abs(nachher-vorher)/(nachher+vorher); + if invers and (hg<>0) then + for i:=smin to smax do // Hintergrund addieren + for j:=0 to pmax do + werte[i*sstep+j*pstep]:= + werte[i*sstep+j*pstep]+hg; + gibAus(inttostr(byte(senkrecht))+' '+inttostr(byte(invers))+' (Parseval-Fehler = '+floattostr(pvFehler)+') ... nämlich von '+floattostr(vorher)+' zu '+floattostr(nachher),1); if in0 then gibAus('Nur Nullen im Input der FFT!',1); if out0 then gibAus('Nur Nullen im Output der FFT!',1); @@ -1137,7 +1235,7 @@ begin // bearbeitet nur den Hauptteil (außer erster und mittlerer Zeile/Spalte) for j:=1 to params.tsiz div 2 -1 do begin werte[i+j*params.xsteps]:= sqr(extended(werte[i+j*params.xsteps]-werte[params.xsteps-i+(params.tsiz-j)*params.xsteps])) // Re^2 - +sqr(extended(werte[params.xsteps-i+j*params.xsteps]+werte[i+(params.tsiz-j)*params.xsteps])); // Im^2 + +sqr(extended(werte[params.xsteps-i+j*params.xsteps]+werte[i+(params.tsiz-j)*params.xsteps])); // Im^2 werte[params.xsteps-i+j*params.xsteps]:=werte[i+j*params.xsteps]; werte[i+(params.tsiz-j)*params.xsteps]:=werte[i+j*params.xsteps]; werte[params.xsteps-i+(params.tsiz-j)*params.xsteps]:=werte[i+j*params.xsteps]; @@ -1656,6 +1754,97 @@ begin end; end; +function tLLWerte.ermittleHintergrund: extended; +var + i: int64; +begin + result:=0; + for i:=0 to params.xsteps-1 do + result:=result + werte[i]; + for i:=1 to params.tsiz-2 do + result:=result + werte[i*params.xsteps]; + for i:=0 to params.xsteps-1 do + result:=result + werte[i + (params.tsiz-1)*params.xsteps]; + for i:=1 to params.tsiz-2 do + result:=result + werte[params.xsteps-1 + i*params.xsteps]; + result:= + result/2/(params.xsteps+params.tsiz-2); +end; + +procedure tLLWerte.filtereHoheFrequenzen(betraege: tLLWerte; xFak,yFak: extended); +var + maxima: tIntPointArray; + i,im,j,jm,mCnt: int64; + wert,minWert,maxWert: extended; +begin + setlength(maxima,0); + mCnt:=0; + for j:=0 to params.tsiz div 2 + 1 do begin + jm:=j - 1 + params.tsiz*byte(j=0); + for i:=0 to params.xsteps div 2 + 1 do begin + im:=i - 1 + params.xsteps*byte(i=0); + wert:=betraege.werte[i+j*params.xsteps]; + if + (wert > betraege.werte[im+j*params.xsteps]) and + (wert > betraege.werte[i+jm*params.xsteps]) and + (wert > betraege.werte[(i+1)+j*params.xsteps]) and + (wert > betraege.werte[i+(j+1)*params.xsteps]) then begin + if length(maxima)<=mCnt then + setlength(maxima,length(maxima)+1024); + maxima[mCnt]['x']:=i; + maxima[mCnt]['y']:=j; + inc(mCnt); + end; + end; + end; + setlength(maxima,mCnt); + writeln(length(maxima),' (von ',params.xsteps*params.tsiz,')'); + + betraege.sortiereMaxima(maxima); + mCnt:=1; + maxWert:=0; + for i:=1 to length(maxima)-1 do begin + minWert:=0; + for j:=0 to i-1 do begin + wert:=sqr((maxima[j]['x']-maxima[i]['x'])*xFak)+sqr((maxima[j]['y']-maxima[i]['y'])*yFak); + if (j=0) or (wert<minWert) then begin + minWert:=wert; + if (minWert<=maxWert) and (i<>1) then + break; + end; + end; + if (i=1) or (minWert>maxWert) then begin + maxWert:=minWert; + mCnt:=i; + end; + end; + + im:=params.xsteps div 2 + 1; + jm:=params.tsiz div 2 + 1; + for j:=0 to jm do + for i:=0 to im do begin + wert:=(sqr(i)*xFak+sqr(j)*yFak)/(sqr(maxima[mCnt]['x'])*xFak+sqr(maxima[mCnt]['y'])*yFak); + if wert > 0.6 then + wert:=0 + else if wert > 0.4 then + wert:=(1+cos((wert-0.4)/0.2*pi))/2 + else + wert:=1; + + werte[i+j*params.xsteps]:= + werte[i+j*params.xsteps]*wert; + if (i>0) and (i<im) then + werte[params.xsteps-i+j*params.xsteps]:= + werte[params.xsteps-i+j*params.xsteps]*wert; + if (j>0) and (j<jm) then + werte[i+(params.tsiz-j)*params.xsteps]:= + werte[i+(params.tsiz-j)*params.xsteps]*wert; + if (i>0) and (i<im) and (j>0) and (j<jm) then + werte[params.xsteps-i+(params.tsiz-j)*params.xsteps]:= + werte[params.xsteps-i+(params.tsiz-j)*params.xsteps]*wert; + end; +end; + // tWavelet ******************************************************************** function tWavelet.setzeTyp(s: string): boolean; @@ -1711,7 +1900,7 @@ begin end; gibAus(inttostr(werte.params.xsteps)+' '+inttostr(werte.params.tsiz)+' '+inttostr(length(werte.werte)),1); tmpFFTAlgo:=createFFTAlgorithmus(2*hlen,doRes,doResIms); - werte.fft(true,false,tmpFFTAlgo,nil,pvFehler); + werte.fft(true,false,tmpFFTAlgo,nil,0,pvFehler); tmpFFTAlgo.free; end; wtFrequenzfenster: begin |