tFenster ist jetzt Klasse; zeitfrequenzanalyse funktioniert wieder

1 files changed, 151 insertions, 104 deletions

diff --git a/epostunit.pas b/epostunit.pas
index 28e4b44..581cccc 100644
--- a/epostunit.pas
+++ b/epostunit.pas
@@ -2631,153 +2631,193 @@ end;
 
 function tWerte.berechneZeitfrequenzanalyse(st: boolean; f: tMyStringlist; threads: longint; quelle: tWerte; Warn: tWarnstufe): boolean;
 var
-  i,j,tmin,tmax,tOf,Schritt: longint;
-  Zeit,pvFehler,freqMax:     extended;
-  Fenster:                   tFenster;
-  s:                         string;
+  i,j,tmin,tmax,qlen:  longint;
+  Zeit,pvFehler,total: extended;
+  wavelet:             tGauszFenster;
+  fenster:             tSin2Fenster;
+  s:                   string;
+  senkrecht,rundeAuf2: boolean;
 begin
   result:=false;
-  if (not st) and (quelle._xsteps<>1) and (quelle._tsiz<>1) then begin
+  if not ((quelle._xsteps<>1) xor (quelle._tsiz<>1)) then begin
     gibAus('Eine Zeitfrequenzanalyse geht nur auf eindimensionalen Daten! ('+inttostr(quelle._xsteps)+'x'+inttostr(quelle._tsiz)+')',3);
     exit;
   end;
   warteaufBeendigungDesLeseThreads;
-  Zeit:=now;
-  Fenster.Breite:=0;
-  Fenster.aktiv:=false;
-  Fenster.Rand:=0;
-  Schritt:=round(sqrt(quelle._tsiz));
+  senkrecht:=quelle._tsiz<>1;
   tmin:=0;
-  tmax:=quelle._tsiz-1;
-  if (quelle._xsteps=1) then
-    freqMax:=quelle._tsiz/(quelle.Transformationen.tstop-quelle.Transformationen.tstart)
+  rundeAuf2:=false;
+  if senkrecht then
+    tmax:=quelle._tsiz-1
   else
-    freqMax:=quelle._xsteps/(quelle.Transformationen.xstop-quelle.Transformationen.xstart);
+    tmax:=quelle._xsteps-1;
+  Zeit:=now;
+  wavelet:=tGauszFenster.create;
+  fenster:=tSin2Fenster.create;
   Genauigkeit:=gExtended;
   repeat
     if not f.metaReadln(s,true) then begin
       gibAus('Unerwartetes Dateiende!',3);
+      wavelet.free;
+      fenster.free;
       exit;
     end;
     if s='Ende' then break;
-    if startetMit('Threadanzahl:',s) then begin
-      threads:=strtoint(s);
+    if startetMit('Breite:',s) then begin
+      wavelet.breite:=quelle.kont2diskFak(senkrecht,quelle.exprtofloat(st,s));
+      wavelet.aktiv:=true;
       continue;
     end;
-    if startetMit('Fenster:',s) then begin
-      Fenster.Breite:=round(quelle.kont2diskFak('t',quelle.exprtofloat(st,erstesArgument(s,';'))));
-      Fenster.Rand:=round(quelle.kont2diskFak('t',quelle.exprtofloat(st,s)));
-      Fenster.aktiv:=true;
+    if startetMit('Rand:',s) then begin
+      fenster.rand:=round(quelle.kont2diskFak(senkrecht,quelle.exprtofloat(st,s)));
+      fenster.aktiv:=true;
       continue;
     end;
-    if startetMit('Schritt:',s) then begin
-      Schritt:=round(quelle.kont2diskFak('t',quelle.exprtofloat(st,s)));
+    if s='nicht auf Zweierpotenz runden' then begin
+      rundeAuf2:=false;
       continue;
     end;
-    if (quelle._xsteps=1) then begin
-      if startetMit('xmin:',s) then begin
-        tmin:=quelle.kont2disk('x',quelle.exprtofloat(st,s));
-        continue;
-      end;
-      if startetMit('xmax:',s) then begin
-        tmax:=quelle.kont2disk('x',quelle.exprtofloat(st,s));
+    if s='auf Zweierpotenz runden' then begin
+      rundeAuf2:=true;
+      continue;
+    end;
+    if senkrecht then begin
+      if startetMit('tmin:',s) or startetMit('ymin:',s) then begin
+        tmin:=round(quelle.kont2disk('t',quelle.exprtofloat(st,s)));
         continue;
       end;
-      if startetMit('wmax:',s) or startetMit('omegamax:',s) or startetMit('ωmax:',s) then begin
-        freqmax:=quelle.exprtofloat(st,s);
+      if startetMit('tmax:',s) or startetMit('ymax:',s)then begin
+        tmax:=round(quelle.kont2disk('t',quelle.exprtofloat(st,s)));
         continue;
       end;
     end;
-    if (quelle._tsiz=1) then begin
-      if startetMit('tmin:',s) then begin
-        tmin:=quelle.kont2disk('t',quelle.exprtofloat(st,s));
-        continue;
-      end;
-      if startetMit('tmax:',s) then begin
-        tmax:=quelle.kont2disk('t',quelle.exprtofloat(st,s));
+    if senkrecht then begin
+      if startetMit('xmin:',s) then begin
+        tmin:=round(quelle.kont2disk('x',quelle.exprtofloat(st,s)));
         continue;
       end;
-      if startetMit('kmax:',s) then begin
-        freqmax:=quelle.exprtofloat(st,s);
+      if startetMit('xmax:',s) then begin
+        tmax:=round(quelle.kont2disk('x',quelle.exprtofloat(st,s)));
         continue;
       end;
     end;
     gibAus('Verstehe Option '''+s+''' nicht bei Erstellung einer Zeitfrequenzanalyse!',3);
+    wavelet.free;
+    fenster.free;
     exit;
   until false;
-  tmin:=max(0,tmin);
-  tmax:=min(quelle._tsiz-1,tmax);
-  Schritt:=max(1,Schritt);
-
-  _tsiz:=round(power(2,ceil(ln(Fenster.Breite)/ln(2))));
-  _xsteps:=1 + ((tmax-tmin+1-(Fenster.Breite-1)) div Schritt);
-  if _xsteps<=0 then begin
-    gibAus('Das angegebene Fenster passt nicht zwischen Anfangs- und Endzeit! ('+
-      floattostrtrunc(disk2kontFak('t',Fenster.Breite),2,true)+'>='+
-      floattostrtrunc(disk2kontFak('t',tmax-tmin+1),2,true)+')',3);
-    exit;
+
+  if senkrecht then
+    qlen:=quelle._tsiz
+  else
+    qlen:=quelle._xsteps;
+
+  tmin:=max(tmin,0);
+  tmax:=min(tmax,qlen-1);
+  _xsteps:=tmax+1-tmin;
+
+  if rundeAuf2 then begin
+    _xsteps:=round(power(2,round(ln(_xsteps)/ln(2))));
+    if _xsteps>qlen then
+      _xsteps:=_xsteps div 2;
+    tmax:=(tmin+tmax+_xsteps) div 2;
+    tmin:=tmax+1-_xsteps;
+    if tmin<0 then begin
+      tmin:=0;
+      tmax:=_xsteps-1;
+    end;
+    if tmax>=qlen then begin
+      tmax:=qlen-1;
+      tmin:=tmax+1-_xsteps;
+    end;
+    if tmin<0 then begin
+      gibAus('Das Fenster passt nicht in die Werte - das sollte nicht passieren können! ('+
+        inttostr(tmin)+'..'+
+        inttostr(tmax)+': '+
+        inttostr(qlen)+' ('+
+        inttostr(quelle._xsteps)+'x'+
+        inttostr(quelle._tsiz)+'))',3);
+      wavelet.free;
+      fenster.free;
+      exit;
+    end;
   end;
 
-  Transformationen:=tAgglomeration.create;
+  _tsiz:=_xsteps;
+  if wavelet.breite>_xsteps then begin
+    gibAus('Die angegebene Breite ist größer als die Anzahl der Werte! ('+
+      floattostrtrunc(wavelet.breite,2,true)+'>='+
+      inttostr(_xsteps)+')',3);
+    wavelet.free;
+    fenster.free;
+    exit;
+  end;
 
-  if quelle._tsiz<>1 then begin
-    for i:=0 to _xsteps-1 do
-      (Transformationen as tAgglomeration).addKomponente(quelle.Transformationen);
-    (Transformationen as tAgglomeration).schritt:=(quelle.disk2kont('t',tmax)-quelle.disk2kont('t',tmin))*(1+1/(tmax-tmin));
-  end
+  if fenster.aktiv then
+    fenster.breite:=qlen-fenster.rand
   else begin
-    for i:=0 to _tsiz-1 do
-      (Transformationen as tAgglomeration).addKomponente(quelle.Transformationen);
-    (Transformationen as tAgglomeration).schritt:=(quelle.disk2kont('x',tmax)-quelle.disk2kont('x',tmin))*(1+1/(tmax-tmin));
+    fenster.breite:=qlen+1;
+    fenster.rand:=0;
   end;
+  fenster.berechneWerte(qlen,true);
+
+  if senkrecht then
+    Transformationen:=tKoordinatenAusschnitt.create(quelle.Transformationen,0,0,tmin,tmax)
+  else
+    Transformationen:=tKoordinatenAusschnitt.create(quelle.Transformationen,tmin,tmax,0,0);
+  Transformationen:=tDiagonaleAgglomeration.create(Transformationen);
 
   if not st then begin
-    tOf:=(_tsiz-Fenster.Breite) div 2;
+    total:=0;
+
+    tmin:=tmin - _xsteps div 2;
+    while tmin<0 do
+      tmin:=tmin+qlen;
+
     eWerte.holeRam(3);
     gibAus('kopiere Inhalt ...',3);
-    for j:=0 to tOf-1 do
-      for i:=0 to _xsteps-1 do
-        eWerte.werte[i+j*_xsteps]:=0;
-    for j:=tOf+Fenster.Breite to _tsiz-1 do
-      for i:=0 to _xsteps-1 do
-        eWerte.werte[i+j*_xsteps]:=0;
     case quelle.Genauigkeit of
       gSingle:
         for i:=0 to _xsteps-1 do
-          for j:=0 to Fenster.Breite-1 do  begin
-            eWerte.werte[i + (j+tOf)*_xsteps]:=
-              quelle.sWerte.werte[i*Schritt+j+tmin];
-          end;
+          for j:=0 to _tsiz-1 do
+            eWerte.werte[j + i*_xsteps]:=
+              quelle.sWerte.werte[(i+j+tmin) mod qlen] * fenster.werte[(i+j+tmin) mod qlen];
       gDouble:
         for i:=0 to _xsteps-1 do
-          for j:=0 to Fenster.Breite-1 do  begin
-            eWerte.werte[i + (j+tOf)*_xsteps]:=
-              quelle.dWerte.werte[i*Schritt+j+tmin];
-          end;
+          for j:=0 to _tsiz-1 do
+            eWerte.werte[j + i*_xsteps]:=
+              quelle.dWerte.werte[(i+j+tmin) mod qlen] * fenster.werte[(i+j+tmin) mod qlen];
       gExtended:
         for i:=0 to _xsteps-1 do
-          for j:=0 to Fenster.Breite-1 do
-            eWerte.werte[i + (j+tOf)*_xsteps]:=
-              quelle.eWerte.werte[i*Schritt+j+tmin];
+          for j:=0 to _tsiz-1 do
+            eWerte.werte[j + i*_xsteps]:=
+              quelle.eWerte.werte[(i+j+tmin) mod qlen] * fenster.werte[(i+j+tmin) mod qlen];
     end{of case};
-    gibAus('... fertig, berechne Fouriertransformation ...',3);
-    if not fft(threads,true,false,doRes,doBetrQdr,Fenster,pvFehler,Warn) then begin
+    for i:=0 to length(eWerte.werte)-1 do
+      total:=total+sqr(eWerte.werte[i]);
+    gibAus('... fertig ('+floatToStr(total)+'), berechne Fouriertransformation ...',3);
+    if not fft(threads,senkrecht,false,doRes,doBetrQdr,wavelet,pvFehler,Warn) then begin
       gibAus('Es traten Fehler auf!',3);
+      wavelet.free;
+      fenster.free;
       exit;
     end;
-    gibAus(' (Parseval-Fehler = '+floattostr(pvFehler)+')',3);
-  end;
-  Transformationen:=tFFTTransformation.create(Transformationen,false,true);
-  if (Transformationen.tstop<=freqmax) or (freqmax<=0) then
-    _tsiz:=_tsiz div 2
-  else begin
-    freqmax:=Transformationen.tstop * round((_tsiz div 2)/Transformationen.tstop*freqmax) / (_tsiz div 2);
-    _tsiz:=round((_tsiz div 2)/Transformationen.tstop*freqmax);
+    total:=0;
+    for i:=0 to length(eWerte.werte)-1 do
+      total:=total+sqr(eWerte.werte[i]);
+    gibAus(' (Parseval-Fehler = '+floattostr(pvFehler)+') -> '+floattostr(total),3);
   end;
+  Transformationen:=tFFTTransformation.create(Transformationen,not senkrecht,senkrecht);
+  if senkrecht then              // die zweite Hälfte der Werte ist redundant
+    _tsiz:=_tsiz div 2 + 1
+  else
+    _xsteps:=_xsteps div 2 + 1;
   Transformationen:=tKoordinatenAusschnitt.create(Transformationen,0,_xsteps-1,0,_tsiz-1);
   if not st then
     eWerte.holeRAM(0);
   gibAus('... fertig '+timetostr(now-Zeit),3);
+  wavelet.free;
+  fenster.free;
   result:=true;
 end;
 
@@ -3029,7 +3069,7 @@ function tWerte.berechneFFT(st: boolean; f: tMyStringlist; threads: longint; War
 var
   Zeit,pvFehler: extended;
   NB:            tFFTDatenordnung;
-  Fenster:       tFenster;
+  Fenster:       tSin2Fenster;
   senkrecht:     boolean;
   s:             string;
 begin
@@ -3037,13 +3077,13 @@ begin
   warteaufBeendigungDesLeseThreads;
   Zeit:=now;
   NB:=doBetrQdr;
-  Fenster.Breite:=0;
-  Fenster.Rand:=0;
-  Fenster.aktiv:=false;
+  Fenster:=tSin2Fenster.create;
   senkrecht:=true;
   repeat
     if not f.metaReadln(s,true) then begin
       gibAus('Unerwartetes Dateiende!',3);
+      Fenster.free;
+      exit;
     end;
     if s='Ende' then break;
     if startetMit('Nachbereitung:',s) then begin
@@ -3067,6 +3107,7 @@ begin
       continue;
     end;
     gibAus('Verstehe Option '''+s+''' nicht bei Erstellung einer FFT!',3);
+    Fenster.free;
     exit;
   until false;
   if senkrecht then
@@ -3077,6 +3118,7 @@ begin
     gibAus('berechne FFT ...',3);
     if not fft(threads,senkrecht,false,doRes,NB,Fenster,pvFehler,Warn) then begin
       gibAus('Es traten Fehler auf!',3);
+      Fenster.free;
       exit;
     end;
     gibAus(' (Parseval-Fehler = '+floattostr(pvFehler)+')',3);
@@ -3086,6 +3128,7 @@ begin
     eWerte.holeRam(0);
     gibAus('... fertig! '+timetostr(now-Zeit),3);
   end;
+  Fenster.free;
   result:=true;
 end;
 
@@ -3093,7 +3136,7 @@ function tWerte.berechneFFT2d(st: boolean; f: tMyStringlist; threads: longint; W
 var
   Zeit,pvFehler:    extended;
   NB,preOrd:        tFFTDatenordnung;
-  Fensters:         array[boolean] of tFenster;
+  Fensters:         array[boolean] of tSin2Fenster;
   s:                string;
   b,spiegeln:       boolean;
 begin
@@ -3101,15 +3144,15 @@ begin
   warteaufBeendigungDesLeseThreads;
   Zeit:=now;
   NB:=doBetrQdr;
-  for b:=false to true do begin
-    Fensters[b].Breite:=0;
-    Fensters[b].Rand:=0;
-    Fensters[b].aktiv:=false;
-  end;
+  for b:=false to true do
+    Fensters[b]:=tSin2Fenster.create;
   spiegeln:=false;
   repeat
     if not f.metaReadln(s,true) then begin
       gibAus('Unerwartetes Dateiende!',3);
+      for b:=false to true do
+        Fensters[b].free;
+      exit;
     end;
     if s='Ende' then break;
     if startetMit('Nachbereitung:',s) then begin
@@ -3130,6 +3173,8 @@ begin
       continue;
     end;
     gibAus('Verstehe Option '''+s+''' nicht bei Erstellung einer zweidimensionalen FFT!',3);
+    for b:=false to true do
+      Fensters[b].free;
     exit;
   until false;
   Fensters[true].Breite:=_tsiz-Fensters[true].Breite;
@@ -3139,6 +3184,8 @@ begin
 
   if st then begin
     result:=true;
+    for b:=false to true do
+      Fensters[b].free;
     exit;
   end;
 
@@ -3174,6 +3221,8 @@ begin
       fft2dNachbearbeitung(threads,nb);                // die Hauptarbeit
   end{of case};
   gibAus('... fertig! '+timetostr(now-Zeit),3);
+  for b:=false to true do
+    Fensters[b].free;
   result:=true;
 end;
 
@@ -4835,7 +4884,6 @@ var
   i,j,k,hl:        longint;
   sus,suc,tmp,pvF: extended;
   in0,out0:        boolean;
-  fenster:         tFenster;
   tmpW:            tWerte;
   tmpFFTAlgo:      tFFTAlgorithmus;
 begin
@@ -4847,10 +4895,9 @@ begin
   if wl.mitFFT then begin
     for i:=xmi to xma do
       for j:=0 to qu._tsiz-1 do in0:=in0 and (zi.eWerte.werte[i+j*zi._xsteps]=0);
-    fenster.aktiv:=false;
     gibAus('FFT berechnen ...',1);
     tmpFFTAlgo:=createFFTAlgorithmus(zi._tsiz,doRes,doResIms);
-    zi.eWerte.fft(true,false,tmpFFTAlgo,fenster,pvF);
+    zi.eWerte.fft(true,false,tmpFFTAlgo,nil,pvF);
     tmpFFTAlgo.free;
     pvFehler:=pvF+pvFehler;
     if wl.typ=wtSin2 then // Das Sin²-Wavelet besteht eigntlich aus zwei!
@@ -4883,11 +4930,11 @@ begin
     end;
     gibAus('... fertig, iFFT berechnen ...',1);
     tmpFFTAlgo:=createFFTAlgorithmus(zi._tsiz,doResIms,doBetrQdr);
-    zi.eWerte.fft(xmi,xma,true,true,tmpFFTAlgo,fenster,pvF);
+    zi.eWerte.fft(xmi,xma,true,true,tmpFFTAlgo,nil,pvF);
     pvFehler:=pvF+pvFehler;
     case wl.typ of
       wtSin2: begin // Das Sin²-Wavelet besteht eigntlich aus zwei!
-        tmpW.eWerte.fft(true,true,tmpFFTAlgo,fenster,pvF);
+        tmpW.eWerte.fft(true,true,tmpFFTAlgo,nil,pvF);
         pvFehler:=(pvF+pvFehler)/3;
         for i:=xmi to xma do
           for j:=0 to zi._tsiz-1 do begin
@@ -5432,7 +5479,7 @@ begin
     for i:=von+1 to bis do
       if Orte[i]['y'] <> Orte[von]['y'] then begin
         gibAus(' interner Quicksort-Fehler: "komisch, die Orte sind doch unterschiedlich ..."',1);
-        halt;
+        exit;
       end;
     result:=true;
     exit;