From 6dc078214273bf20c778728d612a1a120fa9a708 Mon Sep 17 00:00:00 2001 From: Erich Eckner Date: Wed, 25 Sep 2019 15:51:55 +0200 Subject: berechneAutokorrelation2d sollte nicht Betragsquadrate bilden, sonder korrekt quadrieren! --- epost.lps | 119 ++++++++++++++++++++++---------------------- epostunit.pas | 155 ++++++++++++++++++++++++++++++++++++++++++++-------------- werteunit.pas | 114 ++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 289 insertions(+), 99 deletions(-) diff --git a/epost.lps b/epost.lps index 7a89c3e..5f577eb 100644 --- a/epost.lps +++ b/epost.lps @@ -24,24 +24,24 @@ - - + + - - - + + + - + @@ -50,10 +50,10 @@ - - - + + + @@ -66,7 +66,7 @@ - + @@ -82,18 +82,18 @@ - + - + - + - + @@ -114,12 +114,12 @@ - + - + @@ -210,122 +210,119 @@ - - + + - - + - - + + - - + + - - + - + - + - + - + - + + - + - - + + - - + + - - + + - - + - + - + - - + + - - + + - - + + - - + + - - + + - - + - - + + - - + + - + - - + + diff --git a/epostunit.pas b/epostunit.pas index 0cc4cbf..37476d1 100644 --- a/epostunit.pas +++ b/epostunit.pas @@ -193,6 +193,7 @@ type procedure ermittlePhasenWinkel(threads: longint); procedure entspringe(threads: longint; entspringen: tEntspringModus); procedure fft2dNachbearbeitung(threads: longint; nB: tFFTDatenordnung; znt: boolean); + procedure fft2dQuadrieren(threads: longint; hcc,vcc: boolean); procedure schreibeWert(var f: textFile; p: tExtPoint; var letzterWert: extended; entspringen,verschiebung: extended; skalierung: string; linienIntegral: tLinienIntegral; tmpValues: tKnownValues); function exprToFloat(sT: boolean; s: string): extended; function findeZweitdominantestenPunkt(sT: boolean; f: tMyStringList): boolean; @@ -313,6 +314,13 @@ type constructor create(xMi,xMa: longint; pWerte: tWerte; endordnung: tFFTDatenordnung); procedure stExecute; override; end; + tFFT2dQThread = class(tLogThread) + xMin,xMax: longint; + pW: tWerte; + hcc,vcc: boolean; + constructor create(xMi,xMa: longint; pWerte: tWerte; hComplConj,vComplConj: boolean); + procedure stExecute; override; + end; tTauschThread = class(tLogThread) tMin,tMax: longint; pW: tWerte; @@ -5859,6 +5867,9 @@ var fensters: array[boolean] of tSin2Fenster; s: string; b: boolean; + c: char; + fftRichtungen: array[boolean] of char; + gespiegelt: array[boolean] of boolean; bekannteBefehle: tMyStringList; begin result:=false; @@ -5868,8 +5879,10 @@ begin gibAus('Eine 2d-Autokorrelation braucht (momentan) reelle Werte!',3); exit; end; - for b:=false to true do + for b:=false to true do begin fensters[b]:=tSin2Fenster.create; + gespiegelt[b]:=false; + end; bekannteBefehle:=tMyStringList.create; repeat if not f.metaReadln(s,true) then begin @@ -5881,14 +5894,18 @@ begin end; bekannteBefehle.clear; if istDasBefehl('Ende',s,bekannteBefehle,false) then break; - bekannteBefehle.add('''x-Fenster: ...'''); - bekannteBefehle.add('''t-Fenster: ...'''); - if (pos('-Fenster:',s)=2) and (s[1] in ['x','t']) then begin - b:=s[1]='t'; - erstesArgument(s,':'); - if b then fensters[b].rand:=round(kont2diskFak('t',exprToFloat(sT,s))) - else fensters[b].rand:=round(kont2diskFak('x',exprToFloat(sT,s))); - fensters[b].aktiv:=true; + for c in ['x','y','t'] do + if istDasBefehl(c+'-Fenster:',s,bekannteBefehle,true) then begin + fensters[c<>'x'].rand:=round(kont2diskFak(c,exprToFloat(sT,s))); + fensters[c<>'x'].aktiv:=true; + continue; + end; + if istDasBefehl('horizontal gespiegelt',s,bekannteBefehle,false) then begin + gespiegelt[false]:=true; + continue; + end; + if istDasBefehl('vertikal gespiegelt',s,bekannteBefehle,false) then begin + gespiegelt[true]:=true; continue; end; bekannteBefehle.sort; @@ -5922,40 +5939,39 @@ begin exit; end; - transformationen:=tFFTTransformation.create(transformationen,true,true,false); - gibAus('berechne t-FFT ...',3); - if not fft(threads,true,false,doRes,doResSmi,fensters[true],nil,pvFehler,warn) then begin - gibAus('Es traten Fehler auf!',3); - exit; + // sonst transformieren wir Zeilen, die nur Nullen beinhalten, in der 1. inversen FFT + if gespiegelt[true] and not gespiegelt[false] then begin + fftRichtungen[false]:='x'; + fftRichtungen[true]:='t'; + end + else begin + fftRichtungen[false]:='t'; + fftRichtungen[true]:='x'; end; - gibAus(' (Parseval-Fehler = '+floatToStr(pvFehler)+')',3); - gibAus('... fertig! '+timetostr(now-Zeit),3); - gibAus('berechne x-FFT ...',3); - if not fft(threads,false,false,doRes,doResSmi,fensters[false],nil,pvFehler,warn) then begin - gibAus('Es traten Fehler auf!',3); - exit; + + for b:=false to true do begin + gibAus('berechne '+fftRichtungen[b]+'-FFT ...',3); + if not fft(threads,fftRichtungen[b]='t',false,doRes,doResSmi,fensters[fftRichtungen[b]='t'],nil,pvFehler,warn) then begin + gibAus('Es traten Fehler auf!',3); + exit; + end; + gibAus(' (Parseval-Fehler = '+floatToStr(pvFehler)+')',3); + gibAus('... fertig! '+timetostr(now-Zeit),3); end; - gibAus(' (Parseval-Fehler = '+floatToStr(pvFehler)+')',3); - gibAus('... fertig! '+timetostr(now-Zeit),3); + for b:=false to true do + fensters[b].free; gibAus('Quadrieren ...',3); - fft2dNachbearbeitung(threads,doBetrQdr,false); - gibAus('... fertig! '+timetostr(now-Zeit),3); - gibAus('berechne inverse x-FFT ...',3); - if not fft(threads,false,true,doRes,doRes,nil,nil,pvFehler,warn) then begin - gibAus('Es traten Fehler auf!',3); - exit; - end; - gibAus(' (Parseval-Fehler = '+floatToStr(pvFehler)+')',3); + fft2dQuadrieren(threads,not gespiegelt[false],not gespiegelt[true]); gibAus('... fertig! '+timetostr(now-Zeit),3); - gibAus('berechne inverse t-FFT ...',3); - if not fft(threads,true,true,doRes,doRes,nil,nil,pvFehler,warn) then begin - gibAus('Es traten Fehler auf!',3); - exit; + for b:=true downto false do begin + gibAus('berechne inverse '+fftRichtungen[b]+'-FFT ...',3); + if not fft(threads,fftRichtungen[b]='t',true,doResSmi,doRes,nil,nil,pvFehler,warn) then begin + gibAus('Es traten Fehler auf!',3); + exit; + end; + gibAus(' (Parseval-Fehler = '+floatToStr(pvFehler)+')',3); + gibAus('... fertig! '+timetostr(now-Zeit),3); end; - gibAus(' (Parseval-Fehler = '+floatToStr(pvFehler)+')',3); - gibAus('... fertig! '+timetostr(now-Zeit),3); - for b:=false to true do - fensters[b].free; result:=true; end; @@ -7190,6 +7206,42 @@ begin end; end; +procedure tWerte.fft2dQuadrieren(threads: longint; hcc,vcc: boolean); +var + i: longint; + FQTs: array of tFFT2dQThread; + fertig: boolean; +begin + if istKomplex then + fehler('fft2Quadrieren kann (noch) keine voll-komplexen Werte quadrieren!'); + + // der "Rand" + case genauigkeit of + gSingle: sWerte.fft2dQuadrierenA(hcc,vcc); + gDouble: dWerte.fft2dQuadrierenA(hcc,vcc); + gExtended: eWerte.fft2dQuadrierenA(hcc,vcc); + end{of case}; + + // der Hauptteil (alles außer erster und mittlerer Zeile/Spalte) + setLength(FQTs,threads); + for i:=0 to length(FQTs)-1 do + FQTs[i]:=tFFT2dQThread.create( + round(i*(_xSteps div 2 -1)/length(FQTs))+1, + round((i+1)*(_xSteps div 2 -1)/length(FQTs)), + self, + hcc, + vcc); + repeat + sleep(10); + fertig:=true; + for i:=0 to length(FQTs)-1 do + fertig:=fertig and FQTs[i].fertig; + until fertig; + for i:=0 to length(FQTs)-1 do + FQTs[i].free; + gibAus('Alle FFT2dQThreads fertig!',1); +end; + function tWerte.exprToFloat(sT: boolean; s: string): extended; begin result:=matheunit.exprToFloat(sT,s,knownValues,@callBackGetValue); @@ -8121,6 +8173,33 @@ begin gibAus('... und fertig!',1); end; +// tFFT2dQThread *************************************************************** + +constructor tFFT2dQThread.create(xMi,xMa: longint; pWerte: tWerte; hComplConj,vComplConj: boolean); +begin + inherited create; + xMin:=xMi; + xMax:=xMa; + pW:=pWerte; + hcc:=hComplConj; + vcc:=vComplConj; + gibAus('FFT2d-Quadrierthread kreiert: '+intToStr(xMin)+'-'+intToStr(xMax),1); + suspended:=false; +end; + +procedure tFFT2dQThread.stExecute; +begin + gibAus('FFT2d-Quadrierthread gestartet: '+intToStr(xMin)+'-'+intToStr(xMax)+' ...',1); + if pW.istKomplex then + fehler('tFFT2dNBThread kann (noch) nicht komplex!'); + case pW.genauigkeit of + gSingle: pW.sWerte.fft2dQuadrierenB(xMin,xMax,hcc,vcc); + gDouble: pW.dWerte.fft2dQuadrierenB(xMin,xMax,hcc,vcc); + gExtended: pW.eWerte.fft2dQuadrierenB(xMin,xMax,hcc,vcc); + end{of case}; + gibAus('... und fertig!',1); +end; + // tTauschThread *************************************************************** constructor tTauschThread.create(tMi,tMa: longint; pWerte: tWerte); diff --git a/werteunit.pas b/werteunit.pas index eaa400b..02acddf 100644 --- a/werteunit.pas +++ b/werteunit.pas @@ -63,6 +63,8 @@ type procedure fft2dNachbearbeitungB(xMin,xMax: longint; nB: tFFTDatenordnung); procedure fft2dNachbearbeitungVerdoppeln(nB: tFFTDatenordnung); procedure fft2dNachbearbeitungKomplex(xMin,xMax: longint; nB: tFFTDatenordnung); + procedure fft2dQuadrierenA(hcc,vcc: boolean); + procedure fft2dQuadrierenB(xMin,xMax: longint; hcc,vcc: boolean); procedure tausche(tMin,tMax: longint); procedure schreibeWert(var f: textFile; x,y: extended; beschriftung: tExtPoint; var letzterWert: extended; entspringen, verschiebung: extended; skalierung: string; kvs: tKnownValues; linienIntegral: tLinienIntegral; cbgv: tCallBackGetValue); overload; inline; procedure schreibeWert(var f: textFile; xBeschr,yBeschr,wert: extended; skalierung: string; kvs: tKnownValues; cbgv: tCallBackGetValue); overload; inline; @@ -1274,6 +1276,118 @@ begin end{of case}; end; +procedure tLLWerte.fft2dQuadrierenA(hcc,vcc: boolean); +var + i,j,xS2,tS2XS: longint; + reRe,imRe,reIm: extended; +const + minusPlus: array[boolean] of extended = (-1,1); + zweiNull: array[boolean] of extended = (2,0); +begin + xS2:=params.xSteps div 2; + tS2XS:=(params.tSiz div 2) * params.xSteps; + // die vier "Ecken" sind rein reell + for i:=0 to 1 do + for j:=0 to 1 do + werte[j*tS2XS + i*xS2]:= + sqr(extended(werte[j*tS2XS + i*xS2])); + + // unterste und mittlere Zeile sind reell in t + for j:=0 to 1 do + for i:=1 to xS2-1 do begin + reRe:=werte[i+j*tS2XS]; + imRe:=werte[params.xSteps-i+j*tS2XS]; + + werte[i+j*tS2XS]:= // reRe + minusPlus[hcc]*sqr(imRe)+sqr(reRe); + werte[params.xSteps-i+j*tS2XS]:= // imRe + zweiNull[hcc]*imRe*reRe; + end; + + // linke und mittlere Spalte sind reell in x + for i:=0 to 1 do + for j:=1 to params.tSiz div 2 -1 do begin + reRe:=werte[i*xS2+j*params.xSteps]; + reIm:=werte[i*xS2+(params.tSiz-j)*params.xSteps]; + + werte[i*xS2+j*params.xSteps]:= // reRe + minusPlus[vcc]*sqr(reIm)+sqr(reRe); + werte[i*xS2+(params.tSiz-j)*params.xSteps]:= // reIm + zweiNull[vcc]*reIm*reRe; + end; +end; + +procedure tLLWerte.fft2dQuadrierenB(xMin,xMax: longint; hcc,vcc: boolean); +var + i,j: longint; + reRe,reIm,imRe,imIm: extended; +begin // bearbeitet nur den Hauptteil (außer erster und mittlerer Zeile/Spalte) nach! + case byte(hcc) + 2*byte(vcc) of + 0: + for i:=xMin to xMax do + for j:=1 to params.tSiz div 2 -1 do begin + reRe:=werte[i+j*params.xSteps]; + imRe:=werte[params.xSteps-i+j*params.xSteps]; + reIm:=werte[i+(params.tSiz-j)*params.xSteps]; + imIm:=werte[params.xSteps-i+(params.tSiz-j)*params.xSteps]; + + werte[i+j*params.xSteps]:= // reRe + sqr(imIm)-sqr(imRe)-sqr(reIm)+sqr(reRe); + werte[params.xSteps-i+j*params.xSteps]:= // imRe + -2*imIm*reIm+2*imRe*reRe; + werte[i+(params.tSiz-j)*params.xSteps]:= // reIm + -2*imIm*imRe+2*reIm*reRe; + werte[params.xSteps-i+(params.tSiz-j)*params.xSteps]:= // imIm + 2*imRe*reIm+2*imIm*reRe; + end; + 1: // horizontal komplex konjugiert + for i:=xMin to xMax do + for j:=1 to params.tSiz div 2 -1 do begin + reRe:=werte[i+j*params.xSteps]; + imRe:=werte[params.xSteps-i+j*params.xSteps]; + reIm:=werte[i+(params.tSiz-j)*params.xSteps]; + imIm:=werte[params.xSteps-i+(params.tSiz-j)*params.xSteps]; + + werte[i+j*params.xSteps]:= // reRe + -sqr(imIm)+sqr(imRe)-sqr(reIm)+sqr(reRe); + werte[params.xSteps-i+j*params.xSteps]:=0; // imRe + werte[i+(params.tSiz-j)*params.xSteps]:= // reIm + 2*imIm*imRe+2*reIm*reRe; + werte[params.xSteps-i+(params.tSiz-j)*params.xSteps]:=0; // imIm + end; + 2: // vertikal komplex konjugiert + for i:=xMin to xMax do + for j:=1 to params.tSiz div 2 -1 do begin + reRe:=werte[i+j*params.xSteps]; + imRe:=werte[params.xSteps-i+j*params.xSteps]; + reIm:=werte[i+(params.tSiz-j)*params.xSteps]; + imIm:=werte[params.xSteps-i+(params.tSiz-j)*params.xSteps]; + + werte[i+j*params.xSteps]:= // reRe + -sqr(imIm)-sqr(imRe)+sqr(reIm)+sqr(reRe); + werte[params.xSteps-i+j*params.xSteps]:= // imRe + 2*imIm*reIm+2*imRe*reRe; + werte[i+(params.tSiz-j)*params.xSteps]:=0; // reIm + werte[params.xSteps-i+(params.tSiz-j)*params.xSteps]:=0; // imIm + end; + 3: // horizontal & vertikal komplex konjugiert + for i:=xMin to xMax do + for j:=1 to params.tSiz div 2 -1 do begin + reRe:=werte[i+j*params.xSteps]; + imRe:=werte[params.xSteps-i+j*params.xSteps]; + reIm:=werte[i+(params.tSiz-j)*params.xSteps]; + imIm:=werte[params.xSteps-i+(params.tSiz-j)*params.xSteps]; + + werte[i+j*params.xSteps]:= // reRe + sqr(imIm)+sqr(imRe)+sqr(reIm)+sqr(reRe); + werte[params.xSteps-i+j*params.xSteps]:=0; // imRe + werte[i+(params.tSiz-j)*params.xSteps]:=0; // reIm + werte[params.xSteps-i+(params.tSiz-j)*params.xSteps]:= // imIm + -2*imRe*reIm+2*imIm*reRe; + end; + end{of case}; +end; + procedure tLLWerte.tausche(tMin,tMax: longint); var i,j,hLen: longint; -- cgit v1.2.3-70-g09d2