unit werteunit;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, typenunit, math, process;

type
// tLLWerte ********************************************************************
  pTLLWerteSingle = ^tLLWerteSingle;
  pTLLWerteDouble = ^tLLWerteDouble;
  generic tLLWerte<wgen> = class(tObject)
{
    Diese Klasse stellt nur die Berechnungsroutinen und ähnliches bereit,
    nicht jedoch ein (wie auch immer geartetes) Userinterface.
    Auch die korrekte Parallelisierung obliegt dem übergeordneten Programmteil.
}
  private
  public
    werte:                                       array of wgen;
    params:                                      pTExtrainfos;
    constructor create(ps: pTExtrainfos); overload;
    constructor create(original: pTLLWerteSingle; ps: pTExtrainfos; xmin,xmax: longint); overload;
    constructor create(original: pTLLWerteDouble; ps: pTExtrainfos; xmin,xmax: longint); overload;
    procedure kopiereVon(st: boolean; original: pTLLWerteSingle); overload;
    procedure kopiereVon(st: boolean; original: pTLLWerteDouble); overload;
    procedure kopiereVon(st: boolean; original: pTLLWerteSingle; xmin,xmax: longint); overload;
    procedure kopiereVon(st: boolean; original: pTLLWerteDouble; xmin,xmax: longint); overload;
    procedure kopiereVon(st: boolean; original: pTLLWerteSingle; xmin,xmax,tmin,tmax: longint); overload;
    procedure kopiereVon(st: boolean; original: pTLLWerteDouble; xmin,xmax,tmin,tmax: longint); overload;
    procedure kopiereVerzerrt(original: pTLLWerteSingle; ZPs: tIntPointArray; ZGs: tExtPointArray; ZAs: tExtendedArray; xmin,xmax,tmin,tmax: longint; vb,nb: tTransformationen); overload;
    procedure kopiereVerzerrt(original: pTLLWerteDouble; ZPs: tIntPointArray; ZGs: tExtPointArray; ZAs: tExtendedArray; xmin,xmax,tmin,tmax: longint; vb,nb: tTransformationen); overload;
    destructor destroy; override;
    function liesDateien(dateien: tGenerischeInputDateiInfoArray): boolean;
    function fft(senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fen: tFenster; out pvFehler: extended): boolean; overload;
    function fft(xmin,xmax,tmin,tmax: longint; senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fen: tFenster; out pvFehler: extended): boolean; overload;
    procedure spiegle; overload;
    procedure spiegle(tmin,tmax: longint); overload;
    procedure fft2dNachbearbeitungA(nb: tFFTDatenordnung);
    procedure fft2dNachbearbeitungB(xmin,xmax: longint; nb: tFFTDatenordnung);
    procedure schreibeWert(var f: textfile; x,y: longint);
    function exprtofloat(st: boolean; s: string; cbgv: tCallBackGetValue): extended;
    procedure holeRam; overload;
    procedure holeRam(ausgaben: byte); overload;
    procedure holeRam(ausgaben: byte; gemaeszTXMinMax: boolean); overload;
    procedure gibMinMaxDichten(out wMi,wMa: extended; xmin,xmax,tmin,tmax: longint);
    function zuPixelWerten(whoehe,wbreite,xpmi,xmi,tmi: longint; xz,yz: extended; pPWerte: pTExtendedArray; pPAnzahlen: pTLongintArray): boolean;
    function findeSchwellwerte(xmi,xma,tmi,tma: longint; Schw: extended): tExtPointArray;
    procedure integriereSingle(qu: pTLLWerteSingle; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung);
    procedure integriereExtended(qu: pTLLWerteDouble; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung);
  end;
  tLLWerteSingle = specialize tLLWerte<single>;
  tLLWerteDouble = specialize tLLWerte<double>;
// andere Typen ****************************************************************
  tWavelet = class(tObject)
    freq,tfwhm,pvFehler: extended;
    hlen:                longint;
    typ:                 tWaveletTyp;
    werte:               tLLWerteDouble;
    mitFFT:              boolean;
    constructor create;
    destructor destroy; override;
    function setzeTyp(s: string): boolean;
    function berechneWerte: boolean;
  end;

implementation

constructor tLLWerte.create(ps: pTExtrainfos);
begin
 inherited create;
 params:=ps;
 setlength(werte,0);
end;

constructor tLLWerte.create(original: pTLLWerteSingle; ps: pTExtrainfos; xmin,xmax: longint);
begin
  inherited create;
  params:=ps;
  kopiereVon(false,original,xmin,xmax);
end;

constructor tLLWerte.create(original: pTLLWerteDouble; ps: pTExtrainfos; xmin,xmax: longint);
begin
  inherited create;
  params:=ps;
  kopiereVon(false,original,xmin,xmax);
end;

procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteSingle);
begin
  kopiereVon(st,original,0,original^.params^.xsteps-1);
end;

procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteDouble);
begin
  kopiereVon(st,original,0,original^.params^.xsteps-1);
end;

procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteSingle; xmin,xmax: longint);
begin
  kopiereVon(st,original,xmin,xmax,0,original^.params^.tsiz-1);
end;

procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteDouble; xmin,xmax: longint);
begin
  kopiereVon(st,original,xmin,xmax,0,original^.params^.tsiz-1);
end;

procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteSingle; xmin,xmax,tmin,tmax: longint);
var i,j: longint;
begin
 inherited create;
 tmax:=min(tmax,original^.params^.tsiz-1);
 tmin:=max(tmin,0);
 params^.tsiz:=tmax+1-tmin;
 xmax:=min(xmax,original^.params^.xsteps-1);
 xmin:=max(xmin,0);
 params^.xsteps:=xmax+1-xmin;
 params^.xstart:=original^.params^.xstart+(original^.params^.xstop-original^.params^.xstart)/(original^.params^.xsteps-1)*xmin;
 params^.xstop:=original^.params^.xstart+(original^.params^.xstop-original^.params^.xstart)/(original^.params^.xsteps-1)*xmax;
 params^.tstart:=original^.params^.tstart+(original^.params^.tstop-original^.params^.tstart)/(original^.params^.tsiz-1)*tmin;
 params^.tstop:=original^.params^.tstart+(original^.params^.tstop-original^.params^.tstart)/(original^.params^.tsiz-1)*tmax;
 params^.maxW:=0;
 params^.minW:=0;
 params^.np:=original^.params^.np;
 params^.beta:=original^.params^.beta;
 if not st then begin
   holeRam(0);
   for i:=xmin to xmax do
     for j:=tmin to tmax do
       werte[i-xmin+(j-tmin)*params^.xsteps]:=original^.werte[i+j*original^.params^.xsteps];
 end;
end;

procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteDouble; xmin,xmax,tmin,tmax: longint);
var i,j: longint;
begin
 inherited create;
 tmax:=min(tmax,original^.params^.tsiz-1);
 tmin:=max(tmin,0);
 params^.tsiz:=tmax+1-tmin;
 xmax:=min(xmax,original^.params^.xsteps-1);
 xmin:=max(xmin,0);
 params^.xsteps:=xmax+1-xmin;
 params^.xstart:=original^.params^.xstart+(original^.params^.xstop-original^.params^.xstart)/(original^.params^.xsteps-1)*xmin;
 params^.xstop:=original^.params^.xstart+(original^.params^.xstop-original^.params^.xstart)/(original^.params^.xsteps-1)*xmax;
 params^.tstart:=original^.params^.tstart+(original^.params^.tstop-original^.params^.tstart)/(original^.params^.tsiz-1)*tmin;
 params^.tstop:=original^.params^.tstart+(original^.params^.tstop-original^.params^.tstart)/(original^.params^.tsiz-1)*tmax;
 params^.maxW:=0;
 params^.minW:=0;
 params^.np:=original^.params^.np;
 params^.beta:=original^.params^.beta;
 if not st then begin
   holeRam(0);
   for i:=xmin to xmax do
     for j:=tmin to tmax do
       werte[i-xmin+(j-tmin)*params^.xsteps]:=original^.werte[i+j*original^.params^.xsteps];
 end;
end;

procedure tLLWerte.kopiereVerzerrt(original: pTLLWerteSingle; ZPs: tIntPointArray; ZGs: tExtPointArray; ZAs: tExtendedArray; xmin,xmax,tmin,tmax: longint; vb,nb: tTransformationen);
var i,j,k,l: longint;
    tmp:     extended;
begin
  for i:=tmin to tmax do
    for j:=xmin to xmax do
      werte[j+i*params^.xsteps]:=0;
  for i:=0 to length(ZPs)-1 do
    for j:=0 to 1 do
      for k:=0 to 1 do
        if (ZPs[i]['x']+j>=xmin) and (ZPs[i]['x']+j<=xmax) and
           (ZPs[i]['y']+k>=tmin) and (ZPs[i]['y']+k<=tmax) then begin
          tmp:=original^.werte[i];
          if (vb.count>0) or (nb.count>0) then
            tmp:=(tmp-original^.params^.minW)/(original^.params^.maxW-original^.params^.minW);
          for l:=0 to vb.count-1 do
            vb[l].transformiereWert(tmp);
          tmp:=tmp * (ZGs[i]['x'] * (2*j-1) + 1-j) * (ZGs[i]['y'] * (2*k-1) + 1-k);
          werte[ZPs[i]['x']+j + (ZPs[i]['y']+k)*params^.xsteps]:=
            werte[ZPs[i]['x']+j + (ZPs[i]['y']+k)*params^.xsteps] +
            tmp;
        end;
  for i:=tmin to tmax do
    for j:=xmin to xmax do begin
      tmp:=werte[j + i*params^.xsteps] / ZAs[j + i*params^.xsteps];
      for k:=0 to nb.count-1 do
        nb[k].transformiereWert(tmp);
      werte[j + i*params^.xsteps]:=tmp;
    end;
end;

procedure tLLWerte.kopiereVerzerrt(original: pTLLWerteDouble; ZPs: tIntPointArray; ZGs: tExtPointArray; ZAs: tExtendedArray; xmin,xmax,tmin,tmax: longint; vb,nb: tTransformationen);
var i,j,k,l: longint;
    tmp:     extended;
begin
 for i:=tmin to tmax do
   for j:=xmin to xmax do
     werte[j+i*params^.xsteps]:=0;
 for i:=0 to length(ZPs)-1 do
   for j:=0 to 1 do
     for k:=0 to 1 do
       if (ZPs[i]['x']+j>=xmin) and (ZPs[i]['x']+j<=xmax) and
          (ZPs[i]['y']+k>=tmin) and (ZPs[i]['y']+k<=tmax) then begin
         tmp:=original^.werte[i];
         if (vb.count>0) or (nb.count>0) then
           tmp:=(tmp-original^.params^.minW)/(original^.params^.maxW-original^.params^.minW);
         for l:=0 to vb.count-1 do
           vb[l].transformiereWert(tmp);
         tmp:=tmp * (ZGs[i]['x'] * (2*j-1) + 1-j) * (ZGs[i]['y'] * (2*k-1) + 1-k);
         werte[ZPs[i]['x']+j + (ZPs[i]['y']+k)*params^.xsteps]:=
           werte[ZPs[i]['x']+j + (ZPs[i]['y']+k)*params^.xsteps] +
           tmp;
       end;
 for i:=tmin to tmax do
   for j:=xmin to xmax do begin
     tmp:=werte[j + i*params^.xsteps] / ZAs[j + i*params^.xsteps];
     for k:=0 to nb.count-1 do
       nb[k].transformiereWert(tmp);
     werte[j + i*params^.xsteps]:=tmp;
   end;
end;

destructor tLLWerte.destroy;
begin
 setlength(werte,0);
 inherited destroy;
end;

function tLLWerte.liesDateien(dateien: tGenerischeInputDateiInfoArray): boolean;
var i,j,k,l,tmpi,etsiz,spAnz,br: longint;
    f:                           file;
    tmps:                        single;
    tmpe,Zeit:                   extended;
    sa:                          tSingleArray;
    ea:                          tExtendedArray;
    ipp:                         tProcess;
    buf:                         tByteArray;
    etwasGelesen:                boolean;
begin
 result:=false;
 gibAus('... Dateien einlesen ...',1);
 zeit:=now;
 tmpi:=0;
 tmps:=0;
 etsiz:=0;
 spAnz:=-1;
 for i:=0 to length(dateien)-1 do begin
   gibAus(' '+dateien[i].Name,1);
   etwasGelesen:=false;
   if dateien[i] is tPipeInputDateiInfo then begin
     if ((dateien[i] as tPipeInputDateiInfo).bytesPerSample<>4) or
        ((dateien[i] as tPipeInputDateiInfo).Kodierung<>k32BitSignedInteger) then begin
       gibAus('Ich kann nur vier Bytes mit einem mal als Integer interpretiert aus einer Pipe einlesen!',3);
       exit;
     end;
     tmpe:=power(2,-31);
     ipp:=tProcess.create(nil);
     ipp.Options:=ipp.Options + [poUsePipes];
     ipp.Executable:=(dateien[i] as tPipeInputDateiInfo).Executable;
     ipp.Parameters.Text:=(dateien[i] as tPipeInputDateiInfo).ParametersText;
     ipp.execute;
     setlength(buf,0);
     br:=0;
     while ipp.running or (ipp.Output.NumBytesAvailable>0) do begin
       if ipp.Output.NumBytesAvailable > 0 then begin
         if br+ipp.Output.NumBytesAvailable>=length(buf) then
           setlength(buf,br+ipp.Output.NumBytesAvailable+65536*1024);
         tmpi:=ipp.Output.Read(buf[br],min(ipp.Output.NumBytesAvailable,length(buf)-br));
         if ((br+tmpi) shr 24) > (br shr 24) then gibAus(inttostr(br+tmpi)+' von '+inttostr(dateien[i].tsiz*dateien[i].xsteps*4)+' Bytes bisher gelesen ('+floattostrtrunc((br+tmpi)/(dateien[i].tsiz*dateien[i].xsteps*4)*100,2,true)+'%).',1);
         br:=br+tmpi;
       end
       else
         sleep(10);
     end;
     ipp.free;
     gibAus('insgesamt '+inttostr(br div 1024 div 1024)+' MB gelesen.',1);
     setlength(buf,br);
     if dateien[i].tsiz*dateien[i].xsteps*4>length(buf) then begin
       gibAus('Ich habe '+inttostr(length(buf))+' Bytes aus der Pipe gelesen, anstelle von wenigstens '+inttostr(dateien[i].tsiz*dateien[i].xsteps*4)+', wie erwartet!',3);
       setlength(buf,0);
       exit;
     end;
     tmpi:=length(buf)-4*dateien[i].tsiz*dateien[i].xsteps;  // der Offset des ersten Daten-Wortes
     for j:=dateien[i].tmin to dateien[i].tmax do
       for k:=dateien[i].xmin to dateien[i].xmax do
         werte[dateien[i].t0abs+ k-dateien[i].xmin + (j-dateien[i].tmin)*(dateien[i].xmax-dateien[i].xmin+1)]:=
           int32((((((buf[tmpi+3+4*(k+j*dateien[i].xsteps)] shl 8) or
                     buf[tmpi+2+4*(k+j*dateien[i].xsteps)]) shl 8) or
                    buf[tmpi+1+4*(k+j*dateien[i].xsteps)]) shl 8) or
                   buf[tmpi+4*(k+j*dateien[i].xsteps)]) * tmpe;
     setlength(buf,0);
     if etwasGelesen then begin
       gibAus('Ich habe diese Runde schon Daten gelesen!',3);
       exit;
     end;
     etwasGelesen:=true;
   end;
   if (dateien[i] is tSpaceTimeInputDateiInfo) or
      (dateien[i] is tTraceInputDateiInfo) then begin
     assign(f,dateien[i].Name);
     reset(f,1);
     blockread(f,tmpi,sizeof(integer));
     dec(tmpi);
     if tmpi-round(params^.tstart/dateien[i].groeszenFaktor)<>i then begin
       gibAus('Datei '''+dateien[i].Name+''' kommt nicht an '+inttostr(i)+'-ter Stelle, wie sie sollte, sondern an '+inttostr(tmpi-round(params^.tstart/dateien[i].groeszenFaktor))+'-ter.',3);
       writeln(tmpi);
       close(f);
       exit;
     end;
     if dateien[i] is tTraceInputDateiInfo then begin
       blockread(f,tmpi,sizeof(integer));  // #Traces
       spAnz:=tmpi;
     end;
     blockread(f,etsiz,sizeof(integer));
     if dateien[i] is tSpaceTimeInputDateiInfo then begin
       for j:=0 to etsiz-1 do begin
         case Dateien[i].Genauigkeit of
           gSingle: begin
             blockread(f,tmps,sizeof(single)); // xstart
             tmpe:=tmps;
           end;
           gExtended:
             blockread(f,tmpe,sizeof(extended)); // xstart
         end{of Case};
         tmpe:=tmpe*dateien[i].groeszenFaktor;
         if j=0 then params^.xstart:=tmpe;
         if tmpe<>params^.xstart then begin
           gibAus('Falscher linker Rand in '''+dateien[i].Name+''' im Schritt '+inttostr(j)+', nämlich '+myfloattostr(tmpe)+' statt '+myfloattostr(params^.xstart)+'!',3);
           close(f);
           exit;
         end;
         case Dateien[i].Genauigkeit of
           gSingle: begin
             blockread(f,tmps,sizeof(single)); // xstop
             tmpe:=tmps;
           end;
           gExtended:
             blockread(f,tmpe,sizeof(extended)); // xstop
         end{of Case};
         tmpe:=tmpe*dateien[i].groeszenFaktor;
         if j=0 then params^.xstop:=tmpe;
         if tmpe<>params^.xstop then begin
           gibAus('Falscher rechter Rand in '''+dateien[i].Name+''' im Schritt '+inttostr(j)+', nämlich '+myfloattostr(tmpe)+' statt '+myfloattostr(params^.xstop)+'!',3);
           close(f);
           exit;
         end;
         blockread(f,tmpi,sizeof(integer)); // xsteps
         if tmpi<>params^.xsteps then begin
           gibAus('Falsche Anzahl an x-Schritten in '''+dateien[i].Name+''' im Schritt '+inttostr(j)+', nämlich '+inttostr(tmpi)+' statt '+myfloattostr(params^.xsteps)+'!',3);
           close(f);
           exit;
         end;
         if ((sizeof(wgen) = sizeof(single)) and (Dateien[i].Genauigkeit=gSingle)) or
            ((sizeof(wgen) = sizeof(double)) and (Dateien[i].Genauigkeit=gExtended)) then
           blockread(f,werte[(j+Dateien[i].t0abs)*params^.xsteps],params^.xsteps*sizeof(wgen))
         else begin
           setlength(sa,params^.xsteps);
           blockread(f,sa[0],params^.xsteps*sizeof(single));
           for k:=0 to params^.xsteps-1 do
             werte[(j+Dateien[i].t0abs)*params^.xsteps+k]:=sa[k];
         end;
         if dateien[i].faktor/sqr(dateien[i].groeszenFaktor)<>1 then
           for k:=0 to params^.xsteps-1 do
             werte[(j+Dateien[i].t0abs)*params^.xsteps+k]:=werte[(j+Dateien[i].t0abs)*params^.xsteps+k]*dateien[i].faktor/sqr(dateien[i].groeszenFaktor);
         if etwasGelesen then begin
           gibAus('Ich habe diese Runde schon Daten gelesen!',3);
           exit;
         end;
         etwasGelesen:=true;
       end;
     end;
     if dateien[i] is tTraceInputDateiInfo then begin
       case Dateien[i].Genauigkeit of
         gSingle: begin
           setlength(sa,etsiz);
           setlength(ea,0);
         end;
         gExtended: begin
           setlength(sa,0);
           setlength(ea,etsiz);
         end;
       end{of case};
       for j:=0 to spAnz-1 do
         case Dateien[i].Genauigkeit of
           gSingle: begin
             blockread(f,tmps,sizeof(single));   // x
             if j=(Dateien[i] as tTraceInputDateiInfo).Spurnummer then begin
               params^.xstop:=tmps;
               params^.xstart:=params^.xstop;
             end;
             for k:=0 to length(FeldgroeszenNamen)-1 do begin
               blockread(f,sa[0],sizeof(single)*length(sa));
               if (j=(Dateien[i] as tTraceInputDateiInfo).Spurnummer) and
                  (k=(Dateien[i] as tTraceInputDateiInfo).Feldnummer) then begin
                 for l:=0 to length(sa)-1 do
                   werte[l+Dateien[i].t0abs]:=sa[l]*dateien[i].faktor/sqr(dateien[i].groeszenFaktor);
                 if etwasGelesen then begin
                   gibAus('Ich habe diese Runde schon Daten gelesen!',3);
                   exit;
                 end;
                 etwasGelesen:=true;
               end;
             end;
           end;
           gExtended: begin
             blockread(f,tmpe,sizeof(extended)); // x
             if j=(Dateien[i] as tTraceInputDateiInfo).Spurnummer then begin
               params^.xstop:=tmpe;
               params^.xstart:=params^.xstop;
             end;
             for k:=0 to length(FeldgroeszenNamen)-1 do begin
               blockread(f,ea[0],sizeof(extended)*length(ea));
               if (j=(Dateien[i] as tTraceInputDateiInfo).Spurnummer) and
                  (k=(Dateien[i] as tTraceInputDateiInfo).Feldnummer) then begin
                 for l:=0 to length(ea)-1 do
                   werte[l+dateien[i].t0abs]:=ea[l]*dateien[i].faktor/sqr(dateien[i].groeszenFaktor);
                 if etwasGelesen then begin
                   gibAus('Ich habe diese Runde schon Daten gelesen!',3);
                   exit;
                 end;
                 etwasGelesen:=true;
               end;
             end;
           end;
         end{of Case};
     end;
     if not eof(f) then begin
       gibAus('Zu viele Daten in '''+dateien[i].Name+'''!',3);
       close(f);
       exit;
     end;
     close(f);
   end;
   if not etwasGelesen then begin
     gibAus('Ich habe diese Runde keine Daten gelesen!',3);
     exit;
   end;
 end;
 gibAus('... fertig '+ZeitDarstellen(now-Zeit),1);
 result:=true;
end;

procedure tLLWerte.gibMinMaxDichten(out wMi,wMa: extended; xmin,xmax,tmin,tmax: longint);
var i,j:    longint;
begin
  wMi:=werte[xmin+tmin*params^.xsteps];
  wMa:=Werte[xmin+tmin*params^.xsteps];
  for i:=xmin to xmax do
    for j:=tmin to tmax do begin
      wMi:=min(wMi,werte[i+j*params^.xsteps]);
      wMa:=max(wMa,werte[i+j*params^.xsteps]);
    end;
end;

function tLLWerte.fft(senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fen: tFenster; out pvFehler: extended): boolean;
var len: longint;
begin
  len:=1;
  if senkrecht then begin
    while 2*len<=params^.tsiz do
      len:=len*2;
    result:=fft(0,params^.xsteps-1,(params^.tsiz-len) div 2,((params^.tsiz+len) div 2) - 1,true,invers,vor,nach,fen,pvFehler);
  end
  else begin
    while 2*len<=params^.xsteps do
      len:=len*2;
    result:=fft((params^.xsteps-len) div 2,((params^.xsteps+len) div 2) - 1,0,params^.tsiz-1,false,invers,vor,nach,fen,pvFehler);
  end;
end;

function tLLWerte.fft(xmin,xmax,tmin,tmax: longint; senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fen: tFenster; out pvFehler: extended): boolean;
var i,j,k,n,dist,absch,wnum,wstep,haL,
      pmax,pmin,smax,smin:             longint;
    in0,out0:                          boolean;
    ims,wRe,wIm:                       tExtendedArray;
    t1,t2,vorher,nachher,fenavg:       extended;
    umsortierung:                      tLongintArray;
const faktoren: array[tFFTDatenordnung,0..2] of longint = // (doResIms,doResSmi,doRes,doBetr,doBetrQdr);
  ((2,1,1),(2,1,1),(1,1,1),(1,1,1),(1,1,1));
begin
  result:=false;
  if senkrecht then begin
    pmax:=tmax;
    pmin:=tmin;
    smax:=xmax;
    smin:=xmin;
  end
  else begin
    pmax:=xmax;
    pmin:=xmin;
    smax:=tmax;
    smin:=tmin;
  end;
  if length(fen.Werte)<>pmax+1-pmin then
    fen.berechneWerte(pmax+1-pmin);
  n:=round(ln(pmax+1-pmin)/ln(2));
  if round(power(2,n))<>pmax+1-pmin then begin
    gibAus(inttostr(pmax+1-pmin)+' ist keine Zweierpotenz, ich bin zu dumm um Werte dieser Anzahl fourierzutransformieren!',3);
    exit;
    faktoren[doRes,0]:=faktoren[doRes,1];  // um so eine dämliche Warnung loszuwerden ...
  end;

  haL:=(pmax+1-pmin) div 2;
  if fen.aktiv then begin
    if senkrecht then begin
      for i:=xmin to xmax do  // Werte fenstern
        for j:=0 to tmax-tmin do
          werte[i+(j+tmin)*params^.xsteps]:=werte[i+(j+tmin)*params^.xsteps]*fen.Werte[j];
    end
    else begin
      for i:=0 to xmax-xmin do
        for j:=tmin to tmax do  // Werte fenstern
          werte[i+xmin+j*params^.xsteps]:=werte[i+xmin+j*params^.xsteps]*fen.Werte[i];
    end;
    fenavg:=0;
    for i:=0 to pmax-pmin do
      fenavg:=fenavg+fen.Werte[i];
    fenavg:=fenavg/(pmax-pmin+1);
    if fenavg<0.5 then begin
      gibAus('Sehr geringer Fensterdurchschnitt: '+floattostr(fenavg)+' ('+inttostr(fen.Breite)+' '+inttostr(length(fen.werte))+' '+inttostr(pmax-pmin+1)+')!',3);
    end;
  end;

  vorher:=0;
  if vor=doBetrQdr then begin
    for i:=xmin to xmax do
      for j:=tmin to tmax do
        vorher:=vorher + werte[i+j*params^.xsteps];
  end
  else begin
    for i:=xmin to xmax do
      for j:=tmin to tmax do
        vorher:=vorher + werte[i+j*params^.xsteps]*werte[i+j*params^.xsteps]
          *faktoren[vor,
                  Byte(((not senkrecht) and (i=xmin)) or (senkrecht and (j=tmin)))
               +2*Byte(((not senkrecht) and (2*i=xmin+xmax+1)) or (senkrecht and (2*j=tmin+tmax+1)))];
  end;

  setlength(umsortierung,pmax+1-pmin);
  for j:=0 to pmax-pmin do begin
    k:=0;
    for i:=0 to n-1 do             // Bits umkehren
      k:=2*k + byte(odd(j shr i));
    umsortierung[j]:=k;
  end;
  setlength(ims,pmax+1-pmin);
  setlength(wRe,haL);
  setlength(wIm,haL);
  for i:=0 to haL-1 do begin
    wRe[i]:=cos(pi*2*i/(pmax+1-pmin));
    wIm[i]:=sin(pi*2*i/(pmax+1-pmin));
  end;

  in0:=true;
  out0:=true;
  for i:=smin to smax do begin
    k:=1-2*byte(invers); // Eingang komplex konjugieren
    case vor of
      doResIms: begin     // Re_0 Re_1=Re_n-1, ... Re_n/2-1=Re_n/2+1 Re_n/2 Im_1=-Im_n-1, ... Im_n/2-1=-Im_n/2+1
        if pmin<>0 then begin // pmin != 0 wird bei dieser Eingangscodierung nicht geduldet (warum auch?)
          gibAus('Bei getrennten Imaginär- und Realteilen wird eine Verschiebung der Eingangsdaten in Transformationsrichtung nicht geduldet - wie kommt soetwas denn zustande?',3);
          exit;
        end;
        if senkrecht then begin
          for j:=1 to haL-1 do begin
            Ims[j]:=k*werte[i+(haL+j)*params^.xsteps];
            Ims[2*haL-j]:=-k*werte[i+(haL+j)*params^.xsteps];
            werte[i+(haL+j)*params^.xsteps]:=werte[i+(haL-j)*params^.xsteps];
          end;
        end
        else begin
          for j:=1 to haL-1 do begin
            Ims[j]:=k*werte[haL+j+i*params^.xsteps];
            Ims[2*haL-j]:=-k*werte[haL+j+i*params^.xsteps];
            werte[haL+j+i*params^.xsteps]:=werte[haL-j+i*params^.xsteps];
          end;
        end;
        ims[0]:=0;
        ims[haL]:=0;
      end;
      doResSmi: begin     // Re_0 Re_1=Re_n-1, ... Re_n/2-1=Re_n/2+1 Re_n/2 -Im_n/2+1=Im_n/2-1, ... -Im_n-1=Im_1
        if pmin<>0 then begin // pmin != 0 wird bei dieser Eingangscodierung nicht geduldet (warum auch?)
          gibAus('Bei getrennten Imaginär- und Realteilen wird eine Verschiebung der Eingangsdaten in Transformationsrichtung nicht geduldet - wie kommt soetwas denn zustande?',3);
          exit;
        end;
        if senkrecht then begin
          for j:=1 to haL-1 do begin
            Ims[j]:=k*werte[i+(2*haL-j)*params^.xsteps];
            Ims[2*haL-j]:=-Ims[j];
            werte[i+(2*haL-j)*params^.xsteps]:=werte[i+(j)*params^.xsteps];
          end;
        end
        else begin
          for j:=1 to haL-1 do begin
            Ims[j]:=k*werte[2*haL-j+i*params^.xsteps];
            Ims[2*haL-j]:=-Ims[j];
            werte[2*haL-j+i*params^.xsteps]:=werte[j+i*params^.xsteps];
          end;
        end;
        ims[0]:=0;
        ims[haL]:=0;
      end;
      doRes:             // im = 0
        for j:=0 to pmax-pmin do
          ims[j]:=0;
      doBetr,doBetrQdr: begin
        gibAus('Ich brauche mehr als Beträge oder Betragsquadrate um eine FFT auzurechnen!',3);
        exit;
      end;
    end{of case};
    if senkrecht then begin
      for j:=0 to tmax-tmin do begin
        if umsortierung[j]>j then begin
          t1:=werte[i+(j+tmin)*params^.xsteps];
          werte[i+(j+tmin)*params^.xsteps]:=werte[i+(umsortierung[j]+tmin)*params^.xsteps];
          werte[i+(umsortierung[j]+tmin)*params^.xsteps]:=t1;
          t1:=ims[j];
          ims[j]:=ims[umsortierung[j]];
          ims[umsortierung[j]]:=t1;
        end;
        in0:=in0 and (werte[i+(j+tmin)*params^.xsteps]=0) and (ims[j]=0);
      end;
    end
    else begin
      for j:=0 to xmax-xmin do begin
        if umsortierung[j]>j then begin
          t1:=werte[j+xmin+i*params^.xsteps];
          werte[j+xmin+i*params^.xsteps]:=werte[umsortierung[j]+xmin+i*params^.xsteps];
          werte[umsortierung[j]+xmin+i*params^.xsteps]:=t1;
          t1:=ims[j];
          ims[j]:=ims[umsortierung[j]];
          ims[umsortierung[j]]:=t1;
        end;
        in0:=in0 and (werte[j+xmin+i*params^.xsteps]=0) and (ims[j]=0);
      end;
    end;
    dist:=1;
    wstep:=haL;
    if senkrecht then begin
      while dist<tmax+1-tmin do begin
        absch:=0;
        while absch<tmax+1-tmin do begin
          wnum:=0;
          for j:=0 to dist-1 do begin
            // x_links:  [absch+j]
            // x_rechts: [absch+j+dist]
            t1:=wRe[wnum]*werte[i+(absch+j+tmin+dist)*params^.xsteps] - wIm[wnum]*ims[absch+j+dist];
            t2:=wRe[wnum]*ims[absch+j+dist] + wIm[wnum]*werte[i+(absch+j+tmin+dist)*params^.xsteps];

            werte[i+(absch+j+tmin+dist)*params^.xsteps]:=werte[i+(absch+j+tmin)*params^.xsteps]-t1;
            ims[absch+j+dist]:=ims[absch+j]-t2;
            werte[i+(absch+j+tmin)*params^.xsteps]:=werte[i+(absch+j+tmin)*params^.xsteps]+t1;
            ims[absch+j]:=ims[absch+j]+t2;
            wnum:=wnum+wstep;
          end;
          absch:=absch+2*dist;
        end;
        wstep:=wstep div 2;
        dist:=dist*2;
      end;
    end
    else begin
      while dist<xmax+1-xmin do begin
        absch:=0;
        while absch<xmax+1-xmin do begin
          wnum:=0;
          for j:=0 to dist-1 do begin
            // x_links:  [absch+j]
            // x_rechts: [absch+j+dist]
            t1:=wRe[wnum]*werte[absch+j+xmin+dist+i*params^.xsteps] - wIm[wnum]*ims[absch+j+dist];
            t2:=wRe[wnum]*ims[absch+j+dist] + wIm[wnum]*werte[absch+j+xmin+dist+i*params^.xsteps];

            werte[absch+j+xmin+dist+i*params^.xsteps]:=werte[absch+j+xmin+i*params^.xsteps]-t1;
            ims[absch+j+dist]:=ims[absch+j]-t2;
            werte[absch+j+xmin+i*params^.xsteps]:=werte[absch+j+xmin+i*params^.xsteps]+t1;
            ims[absch+j]:=ims[absch+j]+t2;
            wnum:=wnum+wstep;
          end;
          absch:=absch+2*dist;
        end;
        wstep:=wstep div 2;
        dist:=dist*2;
      end;
    end;
    if senkrecht then begin
      case nach of   // werte[ i + {tmin..tmax}*xsteps] & ims[{0..tmax-tmin}] -> werte[ i + {0..tmax-tmin}*xsteps ]
        doResIms: begin     // Re_0 Re_1 Re_2 ... Re_n/2 Im_1 Im_2 ... Im_n/2-1
          for j:=0 to haL do     // Realteile zusammensuchen
            werte[i+j*params^.xsteps]:=werte[i+(j+tmin)*params^.xsteps]/sqrt(tmax+1-tmin);
          for j:=1 to haL - 1 do // Imaginärteile zusammensuchen
            werte[i+(j+haL)*params^.xsteps]:=ims[j]/sqrt(tmax+1-tmin);
        end;
        doResSmi: begin     // Re_0 Re_1 Re_2 ... Re_n/2 Im_n/2-1 Im_n/2-2 ... Im_1
          for j:=0 to haL do     // Realteile zusammensuchen
            werte[i+j*params^.xsteps]:=werte[i+(j+tmin)*params^.xsteps]/sqrt(tmax+1-tmin);
          for j:=1 to haL - 1 do // Imaginärteile zusammensuchen
            werte[i+(2*haL-j)*params^.xsteps]:=ims[j]/sqrt(tmax+1-tmin);
        end;
        doBetr:                                  // Abs(Re + i*Im) = Sqrt(Re^2 + Im^2)
          for j:=0 to tmax-tmin do
            werte[i+j*params^.xsteps]:=sqrt((werte[i+(j+tmin)*params^.xsteps]*werte[i+(j+tmin)*params^.xsteps]+sqr(ims[j]))/(tmax+1-tmin));
        doBetrQdr:                               // Re^2 + Im^2
          for j:=0 to tmax-tmin do
            werte[i+j*params^.xsteps]:=(werte[i+(j+tmin)*params^.xsteps]*werte[i+(j+tmin)*params^.xsteps]+sqr(ims[j]))/(tmax+1-tmin);
        doRes:                                   // Re
          for j:=0 to tmax-tmin do
            werte[i+j*params^.xsteps]:=werte[i+(j+tmin)*params^.xsteps]/sqrt(tmax+1-tmin);
        else begin
          gibAus('Flieht, ihr Narren! (diese Meldung sollte eigentlich nicht auftreten können)',3);
          exit;
        end;
      end;
      for j:=0 to tmax-tmin do
        out0:=out0 and (werte[i+j*params^.xsteps]=0);
    end
    else begin
      case nach of   // werte[ {xtmin..xmax}+i*xsteps] & ims[{0..xmax-xmin}] -> werte[ {0..xmax-xmin}+i*xsteps ]
        doResIms: begin     // Re_0 Re_1 Re_2 ... Re_n/2 Im_1 Im_2 ... Im_n/2-1
          for j:=0 to haL do     // Realteile zusammensuchen
            werte[j+i*params^.xsteps]:=werte[j+xmin+i*params^.xsteps]/sqrt(xmax+1-xmin);
          for j:=1 to haL - 1 do // Imaginärteile zusammensuchen
            werte[j+haL+i*params^.xsteps]:=ims[j]/sqrt(xmax+1-xmin);
        end;
        doResSmi: begin     // Re_0 Re_1 Re_2 ... Re_n/2 Im_n/2-1 Im_n/2-2 ... Im_1
          for j:=0 to haL do     // Realteile zusammensuchen
            werte[j+i*params^.xsteps]:=werte[j+xmin+i*params^.xsteps]/sqrt(xmax+1-xmin);
          for j:=1 to haL - 1 do // Imaginärteile zusammensuchen
            werte[2*haL-j+i*params^.xsteps]:=ims[j]/sqrt(xmax+1-xmin);
        end;
        doBetr:                                  // Abs(Re + i*Im) = Sqrt(Re^2 + Im^2)
          for j:=0 to xmax-xmin do
            werte[j+i*params^.xsteps]:=sqrt((werte[j+xmin+i*params^.xsteps]*werte[j+xmin+i*params^.xsteps]+sqr(ims[j]))/(xmax+1-xmin));
        doBetrQdr:                               // Re^2 + Im^2
          for j:=0 to xmax-xmin do
            werte[j+i*params^.xsteps]:=(werte[j+xmin+i*params^.xsteps]*werte[j+xmin+i*params^.xsteps]+sqr(ims[j]))/(xmax+1-xmin);
        doRes:                                   // Re
          for j:=0 to xmax-xmin do
            werte[j+i*params^.xsteps]:=werte[j+xmin+i*params^.xsteps]/sqrt(xmax+1-xmin);
        else begin
          gibAus('Flieht, ihr Narren! (diese Meldung sollte eigentlich nicht auftreten können)',3);
          exit;
        end;
      end;
      for j:=0 to xmax-xmin do
        out0:=out0 and (werte[j+i*params^.xsteps]=0);
    end;
  end;

  gibAus(inttostr(xmin)+'-'+inttostr(xmax)+'x'+inttostr(tmin)+'-'+inttostr(tmax),1);
  if senkrecht then begin
    tmax:=tmax-tmin;
    tmin:=0;
  end
  else begin
    xmax:=xmax-xmin;
    xmin:=0;
  end;
  nachher:=0;
  if nach=doBetrQdr then begin
    for i:=xmin to xmax do
      for j:=tmin to tmax do
        nachher:=nachher + werte[i+j*params^.xsteps];
  end
  else begin
    for i:=xmin to xmax do
      for j:=tmin to tmax do
        nachher:=nachher + werte[i+j*params^.xsteps]*werte[i+j*params^.xsteps]
          *faktoren[nach,
                  Byte(((not senkrecht) and (i=xmin)) or (senkrecht and (j=tmin)))
               +2*Byte(((not senkrecht) and (2*i=xmin+xmax+1)) or (senkrecht and (2*j=tmin+tmax+1)))];
  end;
  if (nachher=0) and (vorher=0) then pvFehler:=0
  else pvFehler:=abs(nachher-vorher)/(nachher+vorher);

  gibAus(inttostr(byte(senkrecht))+' '+inttostr(byte(vor))+' '+inttostr(byte(nach))+' '+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);
  result:=not (in0 or out0);
end;

procedure tLLWerte.schreibeWert(var f: textfile; x,y: longint);
var xa,ta: extended;
begin
 if params^.xstop=params^.xstart then
   xa:=params^.xstop
 else
   xa:=params^.xstart+x/(params^.xsteps-1)*(params^.xstop-params^.xstart);
 if params^.tstop=params^.tstart then
   ta:=params^.tstop
 else
   ta:=params^.tstart+y/(params^.tsiz-1)*(params^.tstop-params^.tstart);
 writeln(f,floattostr(xa)+' '+floattostr(ta)+' '+floattostr(extended(werte[x+y*params^.xsteps])));
end;

procedure tLLWerte.spiegle;
begin
  spiegle(0,params^.tsiz-1);
end;

procedure tLLWerte.spiegle(tmin,tmax: longint);
var i,j:    longint;
    tmp:    wgen;
begin
  for i:=tmin to tmax do
    for j:=0 to params^.xsteps div 2 -1 do begin
      tmp:=werte[j+i*params^.xsteps];
      werte[j+i*params^.xsteps]:=werte[params^.xsteps-j-1+i*params^.xsteps];
      werte[params^.xsteps-j-1+i*params^.xsteps]:=tmp;
    end;
end;

procedure tLLWerte.fft2dNachbearbeitungA(nb: tFFTDatenordnung);
var i: longint;
begin
  case NB of
    doResIms,doResSmi: ;
    doBetr: begin
      werte[0]:=abs(extended(werte[0]));                                       // rein reell
      for i:=1 to params^.xsteps div 2 -1 do begin                             // unterste Zeile ist reell in t
        werte[i]:=sqrt(sqr(extended(werte[i]))+sqr(extended(werte[params^.xsteps-i])));
        werte[params^.xsteps-i]:=werte[i];
      end;
      for i:=1 to params^.xsteps div 2 -1 do begin                             // mittlere Zeile ist reell in t
        werte[i+(params^.tsiz div 2)*params^.xsteps]:=sqrt(sqr(extended(werte[i+(params^.tsiz div 2)*params^.xsteps]))+sqr(extended(werte[params^.xsteps-i+(params^.tsiz div 2)*params^.xsteps])));
        werte[params^.xsteps-i+(params^.tsiz div 2)*params^.xsteps]:=werte[i+(params^.tsiz div 2)*params^.xsteps];
      end;
      werte[params^.xsteps div 2]:=abs(extended(werte[params^.xsteps div 2])); // rein reell
      werte[params^.tsiz div 2]:=abs(extended(werte[params^.tsiz div 2]));     // rein reell
      for i:=1 to params^.tsiz div 2 -1 do begin                               // linkeste Spalte ist reell in x
        werte[i*params^.xsteps]:=sqrt(sqr(extended(werte[i*params^.xsteps]))+sqr(extended(werte[(params^.tsiz-i)*params^.xsteps])));
        werte[(params^.tsiz-i)*params^.xsteps]:=werte[i*params^.xsteps];
      end;
      for i:=1 to params^.tsiz div 2 -1 do begin                               // mittlere Spalte ist reell in x
        werte[params^.xsteps div 2 + i*params^.xsteps]:=sqrt(sqr(extended(werte[params^.xsteps div 2 + i*params^.xsteps]))+sqr(extended(werte[params^.xsteps div 2 + (params^.tsiz-i)*params^.xsteps])));
        werte[params^.xsteps div 2 + (params^.tsiz-i)*params^.xsteps]:=werte[params^.xsteps div 2 + i*params^.xsteps];
      end;
    end;
    doBetrQdr: begin
      werte[0]:=sqr(extended(werte[0]));                                       // rein reell
      for i:=1 to params^.xsteps div 2 -1 do begin                             // unterste Zeile ist reell in t
        werte[i]:=sqr(extended(werte[i]))+sqr(extended(werte[params^.xsteps-i]));
        werte[params^.xsteps-i]:=werte[i];
      end;
      for i:=1 to params^.xsteps div 2 -1 do begin                             // mittlere Zeile ist reell in t
        werte[i+(params^.tsiz div 2)*params^.xsteps]:=sqr(extended(werte[i+(params^.tsiz div 2)*params^.xsteps]))+sqr(extended(werte[params^.xsteps-i+(params^.tsiz div 2)*params^.xsteps]));
        werte[params^.xsteps-i+(params^.tsiz div 2)*params^.xsteps]:=werte[i+(params^.tsiz div 2)*params^.xsteps];
      end;
      werte[params^.xsteps div 2]:=sqr(extended(werte[params^.xsteps div 2])); // rein reell
      werte[params^.tsiz div 2]:=sqr(extended(werte[params^.tsiz div 2]));     // rein reell
      for i:=1 to params^.tsiz div 2 -1 do begin                               // linkeste Spalte ist reell in x
        werte[i*params^.xsteps]:=sqr(extended(werte[i*params^.xsteps]))+sqr(extended(werte[(params^.tsiz-i)*params^.xsteps]));
        werte[(params^.tsiz-i)*params^.xsteps]:=werte[i*params^.xsteps];
      end;
      for i:=1 to params^.tsiz div 2 -1 do begin                               // mittlere Spalte ist reell in x
        werte[params^.xsteps div 2 + i*params^.xsteps]:=sqr(extended(werte[params^.xsteps div 2 + i*params^.xsteps]))+sqr(extended(werte[params^.xsteps div 2 + (params^.tsiz-i)*params^.xsteps]));
        werte[params^.xsteps div 2 + (params^.tsiz-i)*params^.xsteps]:=werte[params^.xsteps div 2 + i*params^.xsteps];
      end;
    end;
  end{of case};
end;

procedure tLLWerte.fft2dNachbearbeitungB(xmin,xmax: longint; nb: tFFTDatenordnung);
var i,j:    longint;
begin // bearbeitet nur den Hauptteil (außer erster und mittlerer Zeile/Spalte) nach!
  case nb of
    doBetr: begin
      for i:=xmin to xmax do
        for j:=1 to params^.tsiz div 2 -1 do begin
          werte[i+j*params^.xsteps]:=
            sqrt(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
          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];
        end;
    end;
    doBetrQdr: begin
     for i:=xmin to xmax do
        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
          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];
        end;
    end;
  end{of case}
end;

function tLLWerte.exprtofloat(st: boolean; s: string; cbgv: tCallBackGetValue): extended;
var i,j,k,l,m:   longint;
    inv,neg,cbv: boolean;
    val1,val2:   extended;
const
    fkt1: array[0..5] of string = ('exp','sin','cos','tan','sqr','sqrt');
    fkt2: array[0..1] of string = ('min','max');
begin
  s:=trimAll(s);

  for i:=0 to length(fkt1)-1 do
    while fktpos(fkt1[i],s)>0 do begin
      j:=fktpos(fkt1[i],s)+length(fkt1[i]);
      while (j<=length(s)) and (s[j]<>'(') do
        inc(j);
      m:=j+1;                // erstes Zeichen innerhalb der Klammern
      k:=1;
      while (j<length(s)) and (k>0) do begin
        inc(j);
        case s[j] of
          '(': inc(k);
          ')': dec(k);
        end;
      end;
      k:=fktpos(fkt1[i],s); // erstes Zeichen des Funktionsnamens
      val1:=exprtofloat(st,copy(s,m,j-m),cbgv);
      case i of
        0: val1:=exp(val1);
        1: val1:=sin(val1);
        2: val1:=cos(val1);
        3: val1:=tan(val1);
        4: val1:=sqr(val1);
        5: val1:=sqrt(val1);
      end{of case};
      s:=copy(s,1,k-1) + floattostr(val1) + copy(s,j+1,length(s));
    end;

  for i:=0 to length(fkt2)-1 do
    while fktpos(fkt2[i],s)>0 do begin
      j:=fktpos(fkt2[i],s)+length(fkt2[i]);
      while (j<=length(s)) and (s[j]<>'(') do
        inc(j);
      m:=j+1;                // erstes Zeichen innerhalb der Klammern
      k:=1;
      l:=-1;
      while (j<length(s)) and (k>0) do begin
        if (k=1) and (s[j] in [',',';']) then
          l:=j;              // das Komma/Semikolon
        inc(j);
        case s[j] of
          '(': inc(k);
          ')': dec(k);
        end;
      end;
      k:=fktpos(fkt1[i],s); // erstes Zeichen des Funktionsnamens
      val1:=exprtofloat(st,copy(s,m,l-m),cbgv);
      val2:=exprtofloat(st,copy(s,l+1,j-l-1),cbgv);
      case i of
        0: val1:=min(val1,val2);
        1: val1:=max(val1,val2);
      end{of case};
      s:=copy(s,1,k-1) + floattostr(val1) + copy(s,j+1,length(s));
    end;

  while pos('(',s)>0 do begin
    i:=pos('(',s);
    j:=1;
    while j>0 do begin
      inc(i);
      case s[i] of
        '(': inc(j);
        ')': dec(j);
      end;
    end;
    s:=copy(s,1,pos('(',s)-1)+
       floattostr(exprtofloat(st,copy(s,pos('(',s)+1,i-pos('(',s)-1),cbgv))+
       copy(s,i+1,length(s)-i);
  end;
  if (binOpPos('+',s)>0) or (binOpPos('-',s)>0) then begin
    result:=0;
    inv:=false;
    repeat
      i:=min(binOpPos('+',s),binOpPos('-',s));
      if i=0 then i:=max(binOpPos('+',s),binOpPos('-',s));
      if i=0 then i:=length(s)+1;
      if inv then result:=result-exprtofloat(st,copy(s,1,i-1),cbgv)
             else result:=result+exprtofloat(st,copy(s,1,i-1),cbgv);
      inv:=s[i-byte(i>length(s))]='-';
      delete(s,1,i);
    until s='';
    exit;
  end;
  if (binOpPos('*',s)>0) or (binOpPos('/',s)>0) then begin
    result:=1;
    inv:=false;
    repeat
      i:=min(binOpPos('*',s),binOpPos('/',s));
      if i=0 then i:=max(binOpPos('*',s),binOpPos('/',s));
      if i=0 then i:=length(s)+1;
      if inv then result:=result/exprtofloat(st,copy(s,1,i-1),cbgv)
             else result:=result*exprtofloat(st,copy(s,1,i-1),cbgv);
      inv:=s[i-byte(i>length(s))]='/';
      delete(s,1,i);
    until s='';
    exit;
  end;
  if binOpPos('^',s)>0 then begin
    i:=binOpPos('^',s);
    result:=power(exprtofloat(st,copy(s,1,i-1),cbgv),
                  exprtofloat(st,copy(s,i+1,length(s)-i),cbgv));
    exit
  end;
  neg:=startetMit('-',s);
  cbv:=false;
  for i:=1 to length(s) do
    cbv:=cbv or not (s[i] in ['0'..'9','.',',','e','E']);
  if not cbv then result:=strtofloat(s)
  else if s='np' then result:=params^.np
  else if s='maxw' then result:=params^.maxW
  else if s='minw' then result:=params^.minW
  else if s='beta' then result:=params^.beta
  else if s='xstart' then result:=params^.xstart
  else if s='xstop' then result:=params^.xstop
  else if s='tstart' then result:=params^.tstart
  else if s='tstop' then result:=params^.tstop
  else if st then result:=1
  else if assigned(cbgv) then result:=cbgv(s)
  else result:=nan;

  result:=result*(2*byte(not neg)-1);
end;

procedure tLLWerte.holeRam;
begin
  holeRam(1);
end;

procedure tLLWerte.holeRam(ausgaben: byte);
begin
  holeRam(ausgaben,false);
end;

procedure tLLWerte.holeRam(ausgaben: byte; gemaeszTXMinMax: boolean);
var Zeit:  extended;
    br,ho: longint;
begin
  Zeit:=now;
  if gemaeszTXMinMax then begin
    br:=params^.xsteps_;
    ho:=params^.tsiz_;
  end
  else begin
    br:=params^.xsteps;
    ho:=params^.tsiz;
  end;
  gibAus('Fordere '+inttostr(floor(ho*br*sizeof(wgen)/1024/1024))+' MB RAM an ('+inttostr(br)+' x-Schritte mal '+inttostr(ho)+' t-Schritte; bisher '+inttostr(belegterSpeicher div 1024)+' MB belegt). ...',ausgaben);
  setlength(werte,br*ho);
  gibAus('... fertig '+ZeitDarstellen(now-Zeit),ausgaben);
end;

function tLLWerte.zuPixelWerten(whoehe,wbreite,xpmi,xmi,tmi: longint; xz,yz: extended; pPWerte: pTExtendedArray; pPAnzahlen: pTLongintArray): boolean;
var i,j,k,l,
      xv,xb,tv,tb: longint;
    b:             boolean;
begin
  result:=false;
  for i:=0 to length(pPWerte^)-1 do begin
    pPWerte^[i]:=0;
    pPAnzahlen^[i]:=0;
  end;
  b:=false;
  for j:=0 to whoehe-1 do
    for i:=0 to wbreite-1 do begin
      xv:=min(params^.xsteps-1,max(0,ceil((i+max(0,xpmi)-1/2)/xz+xmi)));
      xb:=min(params^.xsteps-1,max(0,ceil((i+max(0,xpmi)+1/2)/xz+xmi-1)));
      tv:=min(params^.tsiz-1,max(0,ceil((j-1/2)/yz+tmi)));
      tb:=min(params^.tsiz-1,max(0,ceil((j+1/2)/yz+tmi-1)));
      if xv>xb then begin
        xv:=xv-byte((i>0) or (xpmi>0));
        xb:=xb+byte((xpmi<0) or ((i=0) and (xpmi=0)));
      end;
      if tv>tb then begin
        tv:=tv-byte(j>0);
        tb:=tb+byte(j=0);
      end;
      if (xv>xb) or (tv>tb) then begin
        gibAus('Keine Inputwerte für Position '+inttostr(i)+':'+inttostr(j)+'!',1);
        exit;
      end;
      pPanzahlen^[j*wbreite+i]:=(xb-xv+1)*(tb-tv+1);
      for k:=xv to xb do
        for l:=tv to tb do begin
          b:=b or (werte[k+l*params^.xsteps]<>0);
          pPWerte^[j*wbreite+i]:=pPWerte^[j*wbreite+i]+min(params^.maxW,max(params^.minW,werte[k+l*params^.xsteps]));
        end;
    end;
  if not b then
    gibAus('Habe nur Nullen im Input!',1);
  result:=true;
end;

function tLLWerte.findeSchwellwerte(xmi,xma,tmi,tma: longint; Schw: extended): tExtPointArray;
var i,j,k,l,m,vz: longint;
    dx,dy,x0,y0:  extended;
begin
  setlength(result,0);
  gibAus('Schwellwerte finden ('+inttostr(xmi)+'-'+inttostr(xma)+') '+floattostr(Schw)+' ...',1);
  i:=0;
  dx:=(params^.xstop-params^.xstart)/(params^.xsteps-1);
  x0:=params^.xstart;
  dy:=(params^.tstop-params^.tstart)/(params^.tsiz-1);
  y0:=params^.tstart;
  for j:=xmi to xma do begin
    for k:=tmi to tma do begin
      vz:=0;
      for l:=0 to 1 do
        for m:=0 to 1 do
          vz:=
             vz or (byte(werte[j-l+(k-m)*params^.xsteps]>=Schw) shl 1) or
                    byte(werte[j-l+(k-m)*params^.xsteps]<=Schw);
      if vz=3 then begin
        if i>=length(result) then
          setlength(result,length(result)+Speicherhappen);
        result[i]['x']:=j*dx+x0;
        result[i]['y']:=k*dy+y0;
        inc(i);
      end;
    end;
    if (xma-j) and $ff = 0 then
      gibAus('x = '+inttostr(j)+' ('+inttostr(xmi)+'-'+inttostr(xma)+')',1);
  end;
  gibAus('... fertig!',1);
  setlength(result,i);
end;

procedure tLLWerte.integriereSingle(qu: pTLLWerteSingle; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung);
var i,j:        longint;
    int,faktor: extended;
begin
  case richtung of
    irHorizontal: begin
      faktor:=(qu^.params^.xstop-qu^.params^.xstart)/(qu^.params^.xsteps-1);
      for i:=tmi to tma do begin
        int:=0;
        for j:=0 to xmi-1 do
          int:=int+qu^.werte[j + i*qu^.params^.xsteps];
        for j:=xmi to xma do begin
          int:=int+qu^.werte[j + i*qu^.params^.xsteps];
          werte[j-xof + (i-tof)*params^.xsteps]:=int*faktor;
        end;
      end;
    end;
    irEinfall: begin
      faktor:=
        sqrt(
          sqr((qu^.params^.xstop-qu^.params^.xstart)/(qu^.params^.xsteps-1)) +
          sqr((qu^.params^.tstop-qu^.params^.tstart)/(qu^.params^.tsiz-1)));
      gibAus('dx = '+floattostr((qu^.params^.xstop-qu^.params^.xstart)/(qu^.params^.xsteps-1)),1);
      gibAus('dt = '+floattostr((qu^.params^.tstop-qu^.params^.tstart)/(qu^.params^.tsiz-1)),1);
      for i:=tmi to tma do begin  // von links eintretendes (inkl. Ecke links unten)
        int:=0;
        for j:=1 to min(xmi,i) do
          int:=int+qu^.werte[xmi-j + (i-j)*qu^.params^.xsteps];
        for j:=0 to min(tma-i,xma-xmi) do begin
          int:=int+qu^.werte[xmi+j + (i+j)*qu^.params^.xsteps];
          werte[j+xmi-xof + (i+j-tof)*params^.xsteps]:=int*faktor;
        end;
      end;
      for i:=xmi+1 to xma do begin  // von unten eintretendes (exkl. Ecke links unten)
        int:=0;
        for j:=1 to min(tmi,i) do
          int:=int+qu^.werte[i-j + (tmi-j)*qu^.params^.xsteps];
        for j:=0 to min(tma-tmi,xma-i) do begin
          int:=int+qu^.werte[i+j + (tmi+j)*qu^.params^.xsteps];
          werte[i+j-xof + (tmi+j-tof)*params^.xsteps]:=int*faktor;
        end;
      end;
    end;
    irAusfall: begin
      faktor:=
        sqrt(
          sqr((qu^.params^.xstop-qu^.params^.xstart)/(qu^.params^.xsteps-1)) +
          sqr((qu^.params^.tstop-qu^.params^.tstart)/(qu^.params^.tsiz-1)));
      gibAus('dx = '+floattostr((qu^.params^.xstop-qu^.params^.xstart)/(qu^.params^.xsteps-1)),1);
      gibAus('dt = '+floattostr((qu^.params^.tstop-qu^.params^.tstart)/(qu^.params^.tsiz-1)),1);
      for i:=tmi to tma do begin  // nach links austretendes (inkl. Ecke links oben)
        int:=0;
        for j:=1 to min(xmi,qu^.params^.tsiz-1-i) do
          int:=int+qu^.werte[xmi-j + (i+j)*qu^.params^.xsteps];
        for j:=0 to min(i-tmi,xma-xmi) do begin
          int:=int+qu^.werte[xmi+j + (i-j)*qu^.params^.xsteps];
          werte[j+xmi-xof + (i-j-tof)*params^.xsteps]:=int*faktor;
        end;
      end;
      for i:=xmi+1 to xma do begin  // nach oben austretendes (exkl. Ecke links oben)
        int:=0;
        for j:=1 to min(qu^.params^.tsiz-1-tma,i) do
          int:=int+qu^.werte[i-j + (tma+j)*qu^.params^.xsteps];
        for j:=0 to min(tma-tmi,xma-i) do begin
          int:=int+qu^.werte[i+j + (tma-j)*qu^.params^.xsteps];
          werte[i+j-xof + (tma-j-tof)*params^.xsteps]:=int*faktor;
        end;
      end;
    end;
  end{of case};
end;

procedure tLLWerte.integriereExtended(qu: pTLLWerteDouble; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung);
var i,j:        longint;
    int,faktor: extended;
begin
  case richtung of
    irHorizontal: begin
      faktor:=(qu^.params^.xstop-qu^.params^.xstart)/(qu^.params^.xsteps-1);
      for i:=tmi to tma do begin
        int:=0;
        for j:=0 to xmi-1 do
          int:=int+qu^.werte[j + i*qu^.params^.xsteps];
        for j:=xmi to xma do begin
          int:=int+qu^.werte[j + i*qu^.params^.xsteps];
          werte[j-xof + (i-tof)*params^.xsteps]:=int*faktor;
        end;
      end;
    end;
    irEinfall: begin
      faktor:=
        sqrt(
          sqr((qu^.params^.xstop-qu^.params^.xstart)/(qu^.params^.xsteps-1)) +
          sqr((qu^.params^.tstop-qu^.params^.tstart)/(qu^.params^.tsiz-1)));
      gibAus('dx = '+floattostr((qu^.params^.xstop-qu^.params^.xstart)/(qu^.params^.xsteps-1)),1);
      gibAus('dt = '+floattostr((qu^.params^.tstop-qu^.params^.tstart)/(qu^.params^.tsiz-1)),1);
      for i:=tmi to tma do begin  // von links eintretendes (inkl. Ecke links unten)
        int:=0;
        for j:=1 to min(xmi,i) do
          int:=int+qu^.werte[xmi-j + (i-j)*qu^.params^.xsteps];
        for j:=0 to min(tma-i,xma-xmi) do begin
          int:=int+qu^.werte[xmi+j + (i+j)*qu^.params^.xsteps];
          werte[j+xmi-xof + (i+j-tof)*params^.xsteps]:=int*faktor;
        end;
      end;
      for i:=xmi+1 to xma do begin  // von unten eintretendes (exkl. Ecke links unten)
        int:=0;
        for j:=1 to min(tmi,i) do
          int:=int+qu^.werte[i-j + (tmi-j)*qu^.params^.xsteps];
        for j:=0 to min(tma-tmi,xma-i) do begin
          int:=int+qu^.werte[i+j + (tmi+j)*qu^.params^.xsteps];
          werte[i+j-xof + (tmi+j-tof)*params^.xsteps]:=int*faktor;
        end;
      end;
    end;
    irAusfall: begin
      faktor:=
        sqrt(
          sqr((qu^.params^.xstop-qu^.params^.xstart)/(qu^.params^.xsteps-1)) +
          sqr((qu^.params^.tstop-qu^.params^.tstart)/(qu^.params^.tsiz-1)));
      gibAus('dx = '+floattostr((qu^.params^.xstop-qu^.params^.xstart)/(qu^.params^.xsteps-1)),1);
      gibAus('dt = '+floattostr((qu^.params^.tstop-qu^.params^.tstart)/(qu^.params^.tsiz-1)),1);
      for i:=tmi to tma do begin  // nach links austretendes (inkl. Ecke links oben)
        int:=0;
        for j:=1 to min(xmi,qu^.params^.tsiz-1-i) do
          int:=int+qu^.werte[xmi-j + (i+j)*qu^.params^.xsteps];
        for j:=0 to min(i-tmi,xma-xmi) do begin
          int:=int+qu^.werte[xmi+j + (i-j)*qu^.params^.xsteps];
          werte[j+xmi-xof + (i-j-tof)*params^.xsteps]:=int*faktor;
        end;
      end;
      for i:=xmi+1 to xma do begin  // nach oben austretendes (exkl. Ecke links oben)
        int:=0;
        for j:=1 to min(qu^.params^.tsiz-1-tma,i) do
          int:=int+qu^.werte[i-j + (tma+j)*qu^.params^.xsteps];
        for j:=0 to min(tma-tmi,xma-i) do begin
          int:=int+qu^.werte[i+j + (tma-j)*qu^.params^.xsteps];
          werte[i+j-xof + (tma-j-tof)*params^.xsteps]:=int*faktor;
        end;
      end;
    end;
  end{of case};
end;

// tWavelet ********************************************************************

function tWavelet.setzeTyp(s: string): boolean;
begin
  result:=true;
  if s='Sin2' then begin
    typ:=wtSin2;
    exit;
  end;
  if s='Frequenzfenster' then begin
    typ:=wtFrequenzfenster;
    exit;
  end;
  gibAus('Kenne Wavelettyp '''+s+''' nicht!',3);
  result:=false;
end;

function tWavelet.berechneWerte: boolean;
var i:       longint;
    tmp:     extended;
    fenster: tFenster;
    nur0:    boolean;
begin
  result:=false;
  werte.params^.xsteps:=2;
  werte.params^.xstart:=0;
  werte.params^.xstop:=1;
  werte.params^.tstart:=0;
  werte.params^.tstop:=1;

  if mitFFT then begin
    if round(power(2,round(ln(werte.params^.tsiz)/ln(2))))<>werte.params^.tsiz then begin
      gibAus('Waveletlänge muss eine Zweierpotenz sein bei Faltung mittels FFT, '+inttostr(werte.params^.tsiz)+' ist das nicht!',3);
      exit;
    end;
    werte.holeRam(1);
    case typ of
      wtSin2: begin
        hlen:=round(tfwhm);
        for i:=0 to hlen do begin
          tmp:=sqr(cos(i*pi/2/tfwhm));
          if i>0 then begin
            werte.werte[(werte.params^.tsiz-i)*2]:=  tmp*cos(-i*freq*2*pi);
            werte.werte[(werte.params^.tsiz-i)*2+1]:=tmp*sin(-i*freq*2*pi);
          end;
          werte.werte[2*i]:=  tmp*cos(i*freq*2*pi);
          werte.werte[2*i+1]:=tmp*sin(i*freq*2*pi);
        end;
        for i:=hlen+1 to werte.params^.tsiz-1-hlen do begin
          werte.werte[2*i]:=0;
          werte.werte[2*i+1]:=0;
        end;
        fenster.aktiv:=false;
    //    fenster.Breite:=0;
    //    fenster.berechneWerte(werte.tsiz);
        gibAus(inttostr(werte.params^.xsteps)+' '+inttostr(werte.params^.tsiz)+' '+inttostr(length(werte.werte)),1);
        werte.fft(0,1,0,werte.params^.tsiz-1,true,false,doRes,doResIms,fenster,pvFehler);
      end;
      wtFrequenzfenster: begin
        hlen:=werte.params^.tsiz div 2;
        for i:=0 to hlen do begin
          werte.werte[2*i]:=Byte(tfwhm*Abs(i/werte.params^.tsiz-freq)<=1); // Re=1 <=> |f-f_Mitte| < Laenge/T_fwhm
          werte.werte[2*i+1]:=0; // Dummy-Wavelet
        end;
        for i:=hlen+1 to 2*hlen-1 do begin
          werte.werte[2*i]:=0;   // Ims = 0
          werte.werte[2*i+1]:=0; // Dummy-Wavelet
        end;
(*        gibAus('dump:',3);
        for i:=0 to werte.tsiz-1 do
          if (i=0) or (i=werte.tsiz-1) or
             (werte.werte[2*i+1]<>0) or
             (werte.werte[2*i]<>werte.werte[2*(i-1)]) or
             (werte.werte[2*i]<>werte.werte[2*(i+1)]) then
          gibAus(inttostr(i)+' '+floattostr(werte.werte[2*i])+' '+floattostr(werte.werte[2*i+1]),3);
        gibAus('ende',3);
        exit; *)
      end;
    end{of case};
  end
  else begin
    if typ <> wtSin2 then begin
      gibAus('Ich kann nur das ''Sin2''-Wavelet im Zeitbereich definieren!',3);
      exit;
    end;
    hlen:=round(tfwhm);
    werte.params^.tsiz:=2*hlen+1;
    setlength(werte.werte,werte.params^.xsteps*werte.params^.tsiz);
    for i:=-hlen to hlen do begin
      tmp:=sqr(cos(i*pi/2/tfwhm));
      werte.werte[2*(i+hlen)]:=  tmp*cos(i*freq*2*pi);
      werte.werte[2*(i+hlen)+1]:=tmp*sin(i*freq*2*pi);
    end;
  end;
  nur0:=true;
  for i:=0 to length(werte.werte)-1 do
    nur0:=nur0 and (werte.werte[i]=0);
  if nur0 then gibAus('Das Wavelet hat nur Nullen!',1);
  result:=not nur0;
end;

constructor tWavelet.create;
var ps: pTExtrainfos;
begin
  inherited create;
  getmem(ps,sizeof(tExtrainfos));
  ps^:=tExtrainfos.create;
  werte:=tLLWerteDouble.create(ps);
  pvFehler:=0;
end;

destructor tWavelet.destroy;
begin
  werte.params^.Free;
  freemem(werte.params,sizeof(tExtrainfos));
  werte.free;
  inherited destroy;
end;

end.