unit werteunit;

{$mode objfpc}{$H+}

interface

uses
  classes, sysutils, typenunit, math, process, lowlevelunit, matheunit, fftunit,
  randomunit;

type
// tLLWerte ********************************************************************
  pTLLWerteSingle = ^tLLWerteSingle;
  pTLLWerteDouble = ^tLLWerteDouble;
  pTLLWerteExtended = ^tLLWerteExtended;
  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
    procedure sortiereNachWert(var positionen: tInt64PointArray; logarithmischesPivot: boolean);
    function dominanzQuadrat(von,nach: tInt64Point; xFak,yFak: extended): extended; inline;
    function abstandsQuadratFuerDominanz(dist: tExtPoint; spitze: tInt64Point; xFak,yFak: extended): extended; inline; overload;
    function abstandsQuadratFuerDominanz(dist,spitze: tInt64Point; xFak,yFak: extended): extended; inline; overload;
    function torusAbstandsQuadrat(von,nach: tInt64Point; metrik: tExtPoint; istVollKomplex: longint): extended; inline;
  public
    werte:                                       array of wGen;
    params:                                      tExtraInfos;
    constructor create(ps: tExtraInfos); overload;
    constructor create(original: pTLLWerteSingle; ps: tExtraInfos; xMin,xMax: longint); overload;
    constructor create(original: pTLLWerteDouble; ps: tExtraInfos; xMin,xMax: longint); overload;
    constructor create(original: pTLLWerteExtended; ps: tExtraInfos; xMin,xMax: longint); overload;
    procedure kopiereVon(sT: boolean; original: pTLLWerteSingle); overload;
    procedure kopiereVon(sT: boolean; original: pTLLWerteDouble); overload;
    procedure kopiereVon(sT: boolean; original: pTLLWerteExtended); 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: pTLLWerteExtended; 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 kopiereVon(sT: boolean; original: pTLLWerteExtended; xMin,xMax,tMin,tMax: longint); overload;
    procedure kopiereVonNach(original: pTLLWerteSingle; qxmin,qxmax,qtmin,qtmax,zxmin,ztmin: longint); overload;
    procedure kopiereVonNach(original: pTLLWerteDouble; qxmin,qxmax,qtmin,qtmax,zxmin,ztmin: longint); overload;
    procedure kopiereVonNach(original: pTLLWerteExtended; qxmin,qxmax,qtmin,qtmax,zxmin,ztmin: longint); overload;
    procedure kopiereVerzerrt(original: pTLLWerteSingle; zPs: tIntPointArray; zGs: tExtPointArray; zAs: tExtendedArray; xMin,xMax,tMin,tMax: longint; vB,nB: tTransformation; vA,nA: longint); overload;
    procedure kopiereVerzerrt(original: pTLLWerteDouble; zPs: tIntPointArray; zGs: tExtPointArray; zAs: tExtendedArray; xMin,xMax,tMin,tMax: longint; vB,nB: tTransformation; vA,nA: longint); overload;
    procedure kopiereVerzerrt(original: pTLLWerteExtended; zPs: tIntPointArray; zGs: tExtPointArray; zAs: tExtendedArray; xMin,xMax,tMin,tMax: longint; vB,nB: tTransformation; vA,nA: longint); overload;
    procedure kopiereLOVerzerrt(original: pTLLWerteSingle; xMin,xMax,tMin,tMax: longint; verhHo,verhVe: extended); overload;
    procedure kopiereLOVerzerrt(original: pTLLWerteDouble; xMin,xMax,tMin,tMax: longint; verhHo,verhVe: extended); overload;
    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; hg: tExtendedArray; out pvFehler: extended): boolean; overload; inline;
    function fft(sMin,sMax: longint; senkrecht,invers: boolean; algo: tFFTAlgorithmus; fen: tFenster; hg: tExtendedArray; out pvFehler: extended): boolean; overload;
    procedure radonTransformation(xMin,xMax: longint; xStep,yStep: extended; qu: pTLLWerteSingle); overload;
    procedure radonTransformation(xMin,xMax: longint; xStep,yStep: extended; qu: pTLLWerteDouble); overload;
    procedure radonTransformation(xMin,xMax: longint; xStep,yStep: extended; qu: pTLLWerteExtended); overload;
    procedure spiegle; overload;
    procedure spiegle(tMin,tMax: longint); overload;
    procedure fft2dNachbearbeitungA(nB: tFFTDatenordnung);
    procedure fft2dNachbearbeitungB(xMin,xMax: longint; nB: tFFTDatenordnung);
    procedure fft2dNachbearbeitungVerdoppeln(nB: tFFTDatenordnung);
    procedure fft2dNachbearbeitungKomplex(xMin,xMax: longint; nB: tFFTDatenordnung);
    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;
    procedure schreibeWertIntegriert(var f: textfile; i: longint; hor: boolean; var letzterWert: extended; entspringen: extended);
    procedure erzeugeBinning(senkrecht,linien: boolean; x0,dX: extended; s: string);
    procedure holeRAM; inline; overload;
    procedure holeRAM(ausgaben: byte); inline; overload;
    procedure holeRAM(ausgaben: byte; gemaeszTXMinMax: boolean); inline; overload;
    procedure gibMinMaxDichten(out wMi,wMa: extended; out pMi,pMa: tInt64Point; out mMi,mMa: boolean; xMin,xMax,tMin,tMax: longint);
    function zuPixelWerten(wHoehe,wBreite,xPMi,xMi,tMi: longint; xZ,yZ: extended; mo: boolean; pPWerte: pTExtendedArray; pPAnzahlen: pTLongintArray): boolean;
    function findeSchwellwerte(xMi,xMa,tMi,tMa: longint; Schw: extended): tExtPointArray;
    procedure integriere(qu: pTLLWerteSingle; xMi,xMa,tMi,tMa,xOf,tOf: longint; richtung: tIntegrationsRichtung); overload;
    procedure integriere(qu: pTLLWerteDouble; xMi,xMa,tMi,tMa,xOf,tOf: longint; richtung: tIntegrationsRichtung); overload;
    procedure integriere(qu: pTLLWerteExtended; xMi,xMa,tMi,tMa,xOf,tOf: longint; richtung: tIntegrationsRichtung); overload;
    procedure gauszFit(amplituden,breiten,positionen,ueberlappe,hintergruende: pTLLWerteExtended; von,bis: longint; senkrecht: boolean; fensterBreite,maxBreite,maxVerschiebung: extended; positionsMitten: tExtendedArray);
    function ermittleRandDurchschnitt: extended;
    function ermittleRandMinimum: extended;
    function ermittleRandPerzentil(p: extended): extended;
    procedure integriereVertikal(xMi,xMa,tMi,tMa: int64; hg: pTExtendedArray);
    procedure integriereVertikalMitRand(xMi,xMa,tMi,tMa,tRa: int64; hg: pTExtendedArray);
    procedure kantenFilter(betraege: tLLWerte; xFak,yFak: extended; filterTyp: tKantenFilterTyp); overload;
    procedure kantenFilter(betraege: tLLWerte; xFak,yFak: extended; filterTyp: tKantenFilterTyp; einseitig: boolean; out maxPos: tInt64Point); overload;
    procedure fenstereWerte(xMi,xMa,tMi,tMa: int64; xFen,tFen: tFenster; hg: tExtendedArray);
    procedure verschiebe(richtung: tInt64Point; xV,xB,tV,tB: longint);
    procedure ermittlePhasenWinkel(xMi,xMa: longint);
    procedure macheKomplex(tMi,tMa: int64; kmm: tKomplexMachModus; mT: tMersenneTwister);
    procedure entspringe(mi,ma: int64; em: tEntspringModus);
    procedure entferneHeiszePixel(schwelle,minusSchwelle,plusSchwelle: extended);
    procedure quotioent(dend: pTLLWerteSingle; sor: pTLLWerteSingle; xMi,xMa,xOf,tMi,tMa,tOf: int64; eps: extended); overload;
    procedure quotioent(dend: pTLLWerteSingle; sor: pTLLWerteDouble; xMi,xMa,xOf,tMi,tMa,tOf: int64; eps: extended); overload;
    procedure quotioent(dend: pTLLWerteSingle; sor: pTLLWerteExtended; xMi,xMa,xOf,tMi,tMa,tOf: int64; eps: extended); overload;
    procedure quotioent(dend: pTLLWerteDouble; sor: pTLLWerteSingle; xMi,xMa,xOf,tMi,tMa,tOf: int64; eps: extended); overload;
    procedure quotioent(dend: pTLLWerteDouble; sor: pTLLWerteDouble; xMi,xMa,xOf,tMi,tMa,tOf: int64; eps: extended); overload;
    procedure quotioent(dend: pTLLWerteDouble; sor: pTLLWerteExtended; xMi,xMa,xOf,tMi,tMa,tOf: int64; eps: extended); overload;
    procedure quotioent(dend: pTLLWerteExtended; sor: pTLLWerteSingle; xMi,xMa,xOf,tMi,tMa,tOf: int64; eps: extended); overload;
    procedure quotioent(dend: pTLLWerteExtended; sor: pTLLWerteDouble; xMi,xMa,xOf,tMi,tMa,tOf: int64; eps: extended); overload;
    procedure quotioent(dend: pTLLWerteExtended; sor: pTLLWerteExtended; xMi,xMa,xOf,tMi,tMa,tOf: int64; eps: extended); overload;
    procedure produkt(f1: pTLLWerteSingle; f2: pTLLWerteSingle; xMi,xMa,xOf,tMi,tMa,tOf: int64; konj: boolean; daO: tFFTDatenordnung); overload;
    procedure produkt(f1: pTLLWerteSingle; f2: pTLLWerteDouble; xMi,xMa,xOf,tMi,tMa,tOf: int64; konj: boolean; daO: tFFTDatenordnung); overload;
    procedure produkt(f1: pTLLWerteSingle; f2: pTLLWerteExtended; xMi,xMa,xOf,tMi,tMa,tOf: int64; konj: boolean; daO: tFFTDatenordnung); overload;
    procedure produkt(f1: pTLLWerteDouble; f2: pTLLWerteSingle; xMi,xMa,xOf,tMi,tMa,tOf: int64; konj: boolean; daO: tFFTDatenordnung); overload;
    procedure produkt(f1: pTLLWerteDouble; f2: pTLLWerteDouble; xMi,xMa,xOf,tMi,tMa,tOf: int64; konj: boolean; daO: tFFTDatenordnung); overload;
    procedure produkt(f1: pTLLWerteDouble; f2: pTLLWerteExtended; xMi,xMa,xOf,tMi,tMa,tOf: int64; konj: boolean; daO: tFFTDatenordnung); overload;
    procedure produkt(f1: pTLLWerteExtended; f2: pTLLWerteSingle; xMi,xMa,xOf,tMi,tMa,tOf: int64; konj: boolean; daO: tFFTDatenordnung); overload;
    procedure produkt(f1: pTLLWerteExtended; f2: pTLLWerteDouble; xMi,xMa,xOf,tMi,tMa,tOf: int64; konj: boolean; daO: tFFTDatenordnung); overload;
    procedure produkt(f1: pTLLWerteExtended; f2: pTLLWerteExtended; xMi,xMa,xOf,tMi,tMa,tOf: int64; konj: boolean; daO: tFFTDatenordnung); overload;
    procedure wertAusUmgebungMitteln(x,y: longint); inline;
    procedure extrahiereKanten(xMi,xMa,tMi,tMa: longint; vert: boolean; expo: int64);
    procedure nullenEinfuegen(hor,vert: boolean);
    procedure skaliere(tMi,tMa: int64; skalierung: string; trafo: tTransformation; kvs: tKnownValues; cbgv: tCallBackGetValue);
  end;
  tLLWerteSingle = specialize tLLWerte<single>;
  tLLWerteDouble = specialize tLLWerte<double>;
  tLLWerteExtended = specialize tLLWerte<extended>;
// andere Typen ****************************************************************
  tWavelet = class(tObject)
    freq,tfwhm,pvFehler: extended;
    hLen:                longint;
    typ:                 tWaveletTyp;
    werte:               tLLWerteDouble;
    mitFFT:              boolean;
    constructor create(globaleWerte: tKnownValues);
    destructor destroy; override;
    function setzeTyp(s: string): boolean;
    function berechneWerte: boolean;
  end;

const
  anzSergeyFelder = 12;

implementation

uses systemunit;

// tLLWerte ********************************************************************

procedure tLLWerte.sortiereNachWert(var positionen: tInt64PointArray; logarithmischesPivot: boolean);
var
  mins,maxs:                tExtendedArray;
  pivot,wert,wertLi,wertRe: extended;
  vons,biss:                tLongintArray;
  tmp:                      tInt64Point;
  i,li,re,cnt:              int64;
begin
  setLength(vons,1);
  vons[0]:=0;
  setLength(biss,1);
  biss[0]:=length(positionen)-1;
  setLength(mins,1);
  setLength(maxs,1);
  mins[0]:=0;
  maxs[0]:=0;
  for i:=vons[0] to biss[0] do begin
    wert:=werte[positionen[i]['x']+positionen[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;

    if logarithmischesPivot then begin
      pivot:=sqrt(maxs[cnt-1]*mins[cnt-1]);
      if (pivot<=mins[cnt-1]) or (pivot>=maxs[cnt-1]) then
        pivot:=(maxs[cnt-1]+mins[cnt-1])/2;
    end
    else
      pivot:=(maxs[cnt-1]+mins[cnt-1])/2;
    if pivot>=maxs[cnt-1] then
      pivot:=mins[cnt-1];

    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[positionen[li]['x']+positionen[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[positionen[re]['x']+positionen[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']:=positionen[re]['x'];
      tmp['y']:=positionen[re]['y'];
      positionen[re]['x']:=positionen[li]['x'];
      positionen[re]['y']:=positionen[li]['y'];
      positionen[li]['x']:=tmp['x'];
      positionen[li]['y']:=tmp['y'];
      inc(li);
      dec(re);
    end;

    vons[cnt]:=li;
    biss[cnt-1]:=re;

    inc(cnt);
  end;

  for i:=1 to length(positionen)-1 do
    if werte[positionen[i-1]['x']+positionen[i-1]['y']*params.xSteps]<werte[positionen[i]['x']+positionen[i]['y']*params.xSteps]then
      fehler('Interner Fehler: Werte sind nicht sortiert, sollten sie aber sein!');
end;

function tLLWerte.dominanzQuadrat(von,nach: tInt64Point; xFak,yFak: extended): extended;
var
  ref:       extended;
  vo,nac,mi: tExtPoint;
  pt:        tInt64Point;
begin
  ref:=werte[von['x']+von['y']*params.xSteps];
  vo:=extPoint(von);
  nac:=extPoint(symmetrischModulo(nach-von,int64Point(params.xSteps,params.tSiz))+von);
  mi:=vo;
  while round(vo-nac)<>int64Point(0,0) do begin
    mi:=0.5*(vo+nac);
    pt:=modulo(round(mi),int64Point(params.xSteps,params.tSiz));
    if werte[pt['x']+pt['y']*params.xSteps]>=ref then
      nac:=mi
    else
      vo:=mi;
  end;
  vo:=symmetrischModulo(extPoint(von)-mi,int64Point(params.xSteps,params.tSiz));
  result:=abstandsQuadratFuerDominanz(vo,von,xFak,yFak);
end;

function tLLWerte.abstandsQuadratFuerDominanz(dist: tExtPoint; spitze: tInt64Point; xFak,yFak: extended): extended;
begin
  result:=(sqr(dist['x']*xFak)+sqr(dist['y']*yFak)) * (1+power(max(werte[spitze['x']+spitze['y']*params.xSteps],0),0.7));
end;

function tLLWerte.abstandsQuadratFuerDominanz(dist,spitze: tInt64Point; xFak,yFak: extended): extended;
begin
  result:=(sqr(dist['x']*xFak)+sqr(dist['y']*yFak)) * (1+power(max(werte[spitze['x']+spitze['y']*params.xSteps],0),0.7));
end;

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

{$DEFINE tLLWerte_create}
constructor tLLWerte.create(original: pTLLWerteSingle; ps: tExtraInfos; xMin,xMax: longint);
{$INCLUDE werteunit.inc}
constructor tLLWerte.create(original: pTLLWerteDouble; ps: tExtraInfos; xMin,xMax: longint);
{$INCLUDE werteunit.inc}
constructor tLLWerte.create(original: pTLLWerteExtended; ps: tExtraInfos; xMin,xMax: longint);
{$INCLUDE werteunit.inc}
{$UNDEF tLLWerte_create.create}

{$DEFINE tLLWerte_kopiereVon}
procedure tLLWerte.kopiereVon(sT: boolean; original: pTLLWerteSingle);
{$INCLUDE werteunit.inc}
procedure tLLWerte.kopiereVon(sT: boolean; original: pTLLWerteDouble);
{$INCLUDE werteunit.inc}
procedure tLLWerte.kopiereVon(sT: boolean; original: pTLLWerteExtended);
{$INCLUDE werteunit.inc}
{$UNDEF tLLWerte_kopiereVon}

{$DEFINE tLLWerte_kopiereVon_xMiMa}
procedure tLLWerte.kopiereVon(sT: boolean; original: pTLLWerteSingle; xMin,xMax: longint);
{$INCLUDE werteunit.inc}
procedure tLLWerte.kopiereVon(sT: boolean; original: pTLLWerteDouble; xMin,xMax: longint);
{$INCLUDE werteunit.inc}
procedure tLLWerte.kopiereVon(sT: boolean; original: pTLLWerteExtended; xMin,xMax: longint);
{$INCLUDE werteunit.inc}
{$UNDEF tLLWerte_kopiereVon_xMiMa}

{$DEFINE tLLWerte_kopiereVon_xTMiMa}
procedure tLLWerte.kopiereVon(sT: boolean; original: pTLLWerteSingle; xMin,xMax,tMin,tMax: longint);
{$INCLUDE werteunit.inc}
procedure tLLWerte.kopiereVon(sT: boolean; original: pTLLWerteDouble; xMin,xMax,tMin,tMax: longint);
{$INCLUDE werteunit.inc}
procedure tLLWerte.kopiereVon(sT: boolean; original: pTLLWerteExtended; xMin,xMax,tMin,tMax: longint);
{$INCLUDE werteunit.inc}
{$UNDEF tLLWerte_kopiereVon_xTMiMa}

{$DEFINE tLLWerte_kopiereVonNach}
procedure tLLWerte.kopiereVonNach(original: pTLLWerteSingle; qxmin,qxmax,qtmin,qtmax,zxmin,ztmin: longint);
{$INCLUDE werteunit.inc}
procedure tLLWerte.kopiereVonNach(original: pTLLWerteDouble; qxmin,qxmax,qtmin,qtmax,zxmin,ztmin: longint);
{$INCLUDE werteunit.inc}
procedure tLLWerte.kopiereVonNach(original: pTLLWerteExtended; qxmin,qxmax,qtmin,qtmax,zxmin,ztmin: longint);
{$INCLUDE werteunit.inc}
{$UNDEF tLLWerte_kopiereVonNach}

{$DEFINE tLLWerte_kopiereVerzerrt}
procedure tLLWerte.kopiereVerzerrt(original: pTLLWerteSingle; zPs: tIntPointArray; zGs: tExtPointArray; zAs: tExtendedArray; xMin,xMax,tMin,tMax: longint; vB,nB: tTransformation; vA,nA: longint);
{$INCLUDE werteunit.inc}
procedure tLLWerte.kopiereVerzerrt(original: pTLLWerteDouble; zPs: tIntPointArray; zGs: tExtPointArray; zAs: tExtendedArray; xMin,xMax,tMin,tMax: longint; vB,nB: tTransformation; vA,nA: longint);
{$INCLUDE werteunit.inc}
procedure tLLWerte.kopiereVerzerrt(original: pTLLWerteExtended; zPs: tIntPointArray; zGs: tExtPointArray; zAs: tExtendedArray; xMin,xMax,tMin,tMax: longint; vB,nB: tTransformation; vA,nA: longint);
{$INCLUDE werteunit.inc}
{$UNDEF tLLWerte_kopiereVerzerrt}

{$DEFINE tLLWerte_kopiereLOVerzerrt}
procedure tLLWerte.kopiereLOVerzerrt(original: pTLLWerteSingle; xMin,xMax,tMin,tMax: longint; verhHo,verhVe: extended);
{$INCLUDE werteunit.inc}
procedure tLLWerte.kopiereLOVerzerrt(original: pTLLWerteDouble; xMin,xMax,tMin,tMax: longint; verhHo,verhVe: extended);
{$INCLUDE werteunit.inc}
procedure tLLWerte.kopiereLOVerzerrt(original: pTLLWerteExtended; xMin,xMax,tMin,tMax: longint; verhHo,verhVe: extended);
{$INCLUDE werteunit.inc}
{$UNDEF tLLWerte_kopiereLOVerzerrt}

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;
  tmpD:                        double;
  tmpE,Zeit:                   extended;
  sA:                          tSingleArray;
  da:                          tDoubleArray;
  ea:                          tExtendedArray;
  iPP:                         tProcess;
  buf:                         tByteArray;
  etwasGelesen:                boolean;
  s:                           string;
begin
  result:=false;
  gibAus('... Dateien einlesen ...',1);
  Zeit:=now;
  tmpI:=0;
  tmpS:=0;
  tmpD:=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);
        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;
            gDouble: begin
              blockread(f,tmpD,sizeOf(double)); // xStart
              tmpE:=tmpD;
            end;
            gExtended:
              blockread(f,tmpE,sizeOf(extended)); // xStart
          end{of Case};
          tmpE:=tmpE*dateien[i].groeszenFaktor;
          if j=0 then params.transformationen.xStart:=tmpE;
          if tmpE<>params.transformationen.xStart then begin
            gibAus('Falscher linker Rand in '''+dateien[i].name+''' im Schritt '+intToStr(j)+', nämlich '+myFloatToStr(tmpE)+' statt '+myFloatToStr(params.transformationen.xStart)+'!',3);
            close(f);
            exit;
          end;
          case dateien[i].genauigkeit of
            gSingle: begin
              blockread(f,tmpS,sizeOf(single)); // xStop
              tmpE:=tmpS;
            end;
            gDouble: begin
              blockread(f,tmpD,sizeOf(double)); // xStop
              tmpE:=tmpD;
            end;
            gExtended:
              blockread(f,tmpE,sizeOf(extended)); // xStop
          end{of Case};
          tmpE:=tmpE*dateien[i].groeszenFaktor;
          if j=0 then params.transformationen.xStop:=tmpE;
          if tmpE<>params.transformationen.xStop then begin
            gibAus('Falscher rechter Rand in '''+dateien[i].name+''' im Schritt '+intToStr(j)+', nämlich '+myFloatToStr(tmpE)+' statt '+myFloatToStr(params.transformationen.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) = wertGroesze(dateien[i].genauigkeit) 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 power(dateien[i].gamma,3)/sqr(dateien[i].groeszenFaktor)<>1 then    // gamma^3 als Skalierungsfaktor für Dichten ?
            for k:=0 to params.xSteps-1 do
              werte[(j+dateien[i].t0Abs)*params.xSteps+k]:=werte[(j+dateien[i].t0Abs)*params.xSteps+k]*power(dateien[i].gamma,3)/sqr(dateien[i].groeszenFaktor);
        end;
        if etwasGelesen then begin
          gibAus('Ich habe diese Runde schon Daten gelesen!',3);
          exit;
        end;
        etwasGelesen:=true;
      end;
      if dateien[i] is tTraceInputDateiInfo then begin
        case dateien[i].genauigkeit of
          gSingle: begin
            setLength(sA,etsiz);
            setLength(da,0);
            setLength(ea,0);
          end;
          gDouble: begin
            setLength(sA,0);
            setLength(da,etsiz);
            setLength(ea,0);
          end;
          gExtended: begin
            setLength(sA,0);
            setLength(da,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.transformationen.xStop:=tmpS;
                params.transformationen.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]*sqr(dateien[i].gamma)/sqr(dateien[i].groeszenFaktor);
                  if etwasGelesen then begin
                    gibAus('Ich habe diese Runde schon Daten gelesen!',3);
                    exit;
                  end;
                  etwasGelesen:=true;
                end;
              end;
            end;
            gDouble: begin
              blockread(f,tmpD,sizeOf(double)); // x
              if j=(dateien[i] as tTraceInputDateiInfo).spurNummer then begin
                params.transformationen.xStop:=tmpD;
                params.transformationen.xStart:=params.xStop;
              end;
              for k:=0 to length(feldGroeszenNamen)-1 do begin
                blockread(f,da[0],sizeOf(double)*length(da));
                if (j=(dateien[i] as tTraceInputDateiInfo).spurNummer) and
                  (k=(dateien[i] as tTraceInputDateiInfo).feldNummer) then begin
                  for l:=0 to length(da)-1 do
                    werte[l+dateien[i].t0Abs]:=da[l]*sqr(dateien[i].gamma)/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.transformationen.xStop:=tmpE;
                params.transformationen.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]*sqr(dateien[i].gamma)/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 dateien[i] is tPhaseSpaceInputDateiInfo then begin
      if i<>0 then begin
        gibAus('Ich kann Phasenraumdateien nicht kaskadieren!',3);
        close(f);
        exit;
      end;
      assign(f,dateien[i].name);
      reset(f,1);
      seek(f,filePos(f)
        + 4*wertGroesze(dateien[i].genauigkeit) // xStart,xStop,tStart,tStop
        + 2*sizeOf(longint));                   // xSteps,tSiz
      blockread(f,werte[0],params.xSteps*params.tSiz*wertGroesze(dateien[i].genauigkeit));
      close(f);
      etwasGelesen:=true;
    end;
    if dateien[i] is tSergeyInputDateiInfo then begin
      etsiz:=dateien[i].tSiz;

      setLength(buf,sizeOf(extended)*anzSergeyFelder);
      tmpI:=wertGroesze(dateien[i].genauigkeit);

      assign(f,dateien[i].name+'traces/traces.dat');
      reset(f,1);

      setLength(buf,tmpI*anzSergeyFelder);

      for j:=0 to etsiz-1 do begin
        if (dateien[i] as tSergeyInputDateiInfo).feldNummer>0 then
          blockread(f,buf[0],(dateien[i] as tSergeyInputDateiInfo).feldNummer*tmpI);
        case dateien[i].genauigkeit of
          gSingle: begin
            blockread(f,tmpS,tmpI);
            werte[j]:=tmpS;
          end;
          gDouble: begin
            blockread(f,tmpD,tmpI);
            werte[j]:=tmpD;
          end;
          gExtended: begin
            blockread(f,tmpE,tmpI);
            werte[j]:=tmpE;
          end;
        end{of Case};
        if (dateien[i] as tSergeyInputDateiInfo).feldNummer<anzSergeyFelder-1 then
          blockread(f,buf[0],(anzSergeyFelder-1-(dateien[i] as tSergeyInputDateiInfo).feldNummer)*tmpI);
      end;

      setLength(buf,0);
      if not eof(f) then begin
        gibAus('Zu viele Daten in '''+dateien[i].name+'traces/traces.dat''!',3);
        close(f);
        exit;
      end;
      close(f);
      etwasGelesen:=true;
    end;
    if dateien[i] is tAndorInputDateiInfo then begin
      etsiz:=dateien[i].tSiz;

      assign(f,dateien[i].name);
      reset(f,1);
      if (dateien[i] as tAndorInputDateiInfo).istHintergrund then
        seek(f,(dateien[i] as tAndorInputDateiInfo).hintergrundStart)
      else
        seek(f,(dateien[i] as tAndorInputDateiInfo).datenStart);
      tmpI:=wertGroesze(dateien[i].genauigkeit);
      if tmpI = sizeOf(wGen) then
        blockread(f,werte[0],dateien[i].xSteps*dateien[i].tSiz*tmpI)
      else case dateien[i].genauigkeit of
        gSingle:
          for j:=0 to dateien[i].xSteps*dateien[i].tSiz-1 do begin
            blockread(f,tmpS,tmpI);
            werte[j]:=tmpS;
          end;
        gDouble:
          for j:=0 to dateien[i].xSteps*dateien[i].tSiz-1 do begin
            blockread(f,tmpD,tmpI);
            werte[j]:=tmpD;
          end;
        gExtended:
          for j:=0 to dateien[i].xSteps*dateien[i].tSiz-1 do begin
            blockread(f,tmpE,tmpI);
            werte[j]:=tmpE;
          end;
      end{of case};
      try
        readAnAndorString(f,s,false);
      except
        s:='falsch';
      end;
      if s<>'' then begin
        gibAus('Syntax-Fehler in '''+dateien[i].name+''': vmtl. zu viele/wenige Daten.',3);
        closeFile(f);
        exit;
      end;
      close(f);
      etwasGelesen:=true;
    end;
    if not etwasGelesen then begin
      gibAus('Ich habe diese Runde keine Daten gelesen!',3);
      exit;
    end;
  end;
  params.refreshKnownValues;
(*  for j:=0 to params.tSiz-1 do
    for i:=0 to params.xSteps-1 do begin
      werte[i+j*params.xSteps]:=100*sqr(sin(i/4+j/40));
      if i<50 then
        werte[i+j*params.xSteps]:=werte[i+j*params.xSteps]*(1-cos(i/50*pi))/2;
      if (params.xSteps-i)<50 then
        werte[i+j*params.xSteps]:=werte[i+j*params.xSteps]*(1-cos((params.xSteps-i)/50*pi))/2;
      if j<50 then
        werte[i+j*params.xSteps]:=werte[i+j*params.xSteps]*(1-cos(j/50*pi))/2;
      if (params.tSiz-j)<50 then
        werte[i+j*params.xSteps]:=werte[i+j*params.xSteps]*(1-cos((params.tSiz-j)/50*pi))/2;
      werte[i+j*params.xSteps]:=werte[i+j*params.xSteps] + 50 + 10*random;
    end; *)
  gibAus('... fertig '+timetostr(now-Zeit),1);
  result:=true;
end;

procedure tLLWerte.gibMinMaxDichten(out wMi,wMa: extended; out pMi,pMa: tInt64Point; out mMi,mMa: boolean; xMin,xMax,tMin,tMax: longint);
var
  i,j:    longint;
begin
  wMi:=werte[xMin+tMin*params.xSteps];
  wMa:=wMi;
  pMi:=int64Point(xMin,tMin);
  pMa:=pMi;
  mMi:=false;
  mMa:=false;
  for i:=xMin to xMax do
    for j:=tMin+byte(i=xMin) to tMax do begin
      if wMi=werte[i+j*params.xSteps] then
        mMi:=true;
      if wMi > werte[i+j*params.xSteps] then begin
        wMi:=werte[i+j*params.xSteps];
        pMi:=int64Point(i,j);
        mMi:=false;
      end;
      if wMa=werte[i+j*params.xSteps] then
        mMa:=true;
      if wMa < werte[i+j*params.xSteps] then begin
        wMa:=werte[i+j*params.xSteps];
        pMa:=int64Point(i,j);
        mMa:=false;
      end;
    end;
end;

function tLLWerte.fft(senkrecht,invers: boolean; algo: tFFTAlgorithmus; fen: tFenster; hg: tExtendedArray; out pvFehler: extended): boolean;
begin
  if senkrecht then
    result:=fft(0,params.xSteps-1,senkrecht,invers,algo,fen,hg,pvFehler)
  else
    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; hg: tExtendedArray; out pvFehler: extended): boolean;
var
  i,j,pMax,pStep,sStep,imShift: longint;
  in0,out0,imPart,vollKomplex:  boolean;
  vorher,nachher,offset:        extended;
begin
  result:=false;
  vollKomplex:=algo.inOrdnung in [doAlleResIms,doAlleResSmi];
  imShift:=(params.tSiz div 2)*params.xSteps;
  if senkrecht then begin
    pMax:=params.tSiz div (1+byte(vollKomplex)) -1;
    pStep:=params.xSteps;
    sStep:=1;
  end
  else begin
    pMax:=params.xSteps-1;
    pStep:=1;
    sStep:=params.xSteps;
  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:=0;
    if not invers then begin
      if length(hg)=1 then
        offset:=hg[0]
      else if length(hg)<>0 then begin
        if length(hg)<>params.xSteps then begin
          gibAus('fft: Der abzuziehende Hintergrund hat die falsche Länge: '+intToStr(length(hg))+' vs. '+intToStr(params.xSteps),1);
          exit;
        end;
        if senkrecht then
          for j:=0 to pMax do
            for i:=sMin to sMax do
              werte[i+j*params.xSteps]:=werte[i+j*params.xSteps] - hg[i]
        else
          for j:=sMin to sMax do
            for i:=0 to pMax do
              werte[i+j*params.xSteps]:=werte[i+j*params.xSteps] - hg[i];
      end;
    end;
    if length(fen.werte)<>pMax+1 then begin
      gibAus('Die Breite des FFT-Fensters ('+intToStr(length(fen.werte))+') ist nicht gleich der Breite der Werte ('+intToStr(pMax+1)+')!',1);
      exit;
    end;
    for imPart:=false to vollKomplex do  // Werte fenstern
      for i:=sMin to sMax do
        for j:=0 to pMax do
          werte[i*sStep+j*pStep + byte(imPart)*imShift]:=
            (werte[i*sStep+j*pStep + byte(imPart)*imShift]-offset*byte(not imPart))*fen.werte[j];
  end;

  if not senkrecht and vollKomplex and (algo.inOrdnung <> algo.outOrdnung) then begin
    gibAus('Eine waagerechte fft von (senkrecht) '+fftDoToStr(algo.inOrdnung)+' zu '+fftDoToStr(algo.outOrdnung)+' ist nicht in situ möglich!',1);
    exit;
  end;

  vorher:=0;
  nachher:=0;
  in0:=true;
  out0:=true;

  case sizeOf(wGen) of
    sizeOf(single):
      for i:=sMin to sMax do begin
        if not senkrecht and (algo.inOrdnung = doAlleResIms) then
          algo.laden(invers,pSingle(@(werte[i*sStep])),pSingle(@(werte[imShift   +   i*sStep])),pStep)
        else if not senkrecht and (algo.inOrdnung = doAlleResSmi) then      // "Smi" bezieht sich auf die Dimension, in der nicht transformiert wird!
          algo.laden(invers,pSingle(@(werte[i*sStep])),pSingle(@(werte[2*imShift-(1+i)*sStep])),pStep)
        else
          algo.laden(invers,pSingle(@(werte[i*sStep])),pStep);
        algo.summen(true,vorher,in0);
        algo.ausfuehren;
        algo.summen(false,nachher,out0);
        if not senkrecht and (algo.outOrdnung = doAlleResIms) then
          algo.speichern(invers,pSingle(@(werte[i*sStep])),pSingle(@(werte[imShift   +   i*sStep])),pStep)
        else if not senkrecht and (algo.outOrdnung = doAlleResSmi) then     // "Smi" bezieht sich auf die Dimension, in der nicht transformiert wird!
          algo.speichern(invers,pSingle(@(werte[i*sStep])),pSingle(@(werte[2*imShift-(1+i)*sStep])),pStep)
        else
          algo.speichern(invers,pSingle(@(werte[i*sStep])),pStep);
      end;
    sizeOf(double):
      for i:=sMin to sMax do begin
        if not senkrecht and (algo.inOrdnung = doAlleResIms) then
          algo.laden(invers,pDouble(@(werte[i*sStep])),pDouble(@(werte[imShift   +   i*sStep])),pStep)
        else if not senkrecht and (algo.inOrdnung = doAlleResSmi) then      // "Smi" bezieht sich auf die Dimension, in der nicht transformiert wird!
          algo.laden(invers,pDouble(@(werte[i*sStep])),pDouble(@(werte[2*imShift-(1+i)*sStep])),pStep)
        else
          algo.laden(invers,pDouble(@(werte[i*sStep])),pStep);
        algo.summen(true,vorher,in0);
        algo.ausfuehren;
        algo.summen(false,nachher,out0);
        if not senkrecht and (algo.outOrdnung = doAlleResIms) then
          algo.speichern(invers,pDouble(@(werte[i*sStep])),pDouble(@(werte[imShift   +   i*sStep])),pStep)
        else if not senkrecht and (algo.outOrdnung = doAlleResSmi) then     // "Smi" bezieht sich auf die Dimension, in der nicht transformiert wird!
          algo.speichern(invers,pDouble(@(werte[i*sStep])),pDouble(@(werte[2*imShift-(1+i)*sStep])),pStep)
        else
          algo.speichern(invers,pDouble(@(werte[i*sStep])),pStep);
      end;
    sizeOf(extended):
      for i:=sMin to sMax do begin
        if not senkrecht and (algo.inOrdnung = doAlleResIms) then
          algo.laden(invers,pExtended(@(werte[i*sStep])),pExtended(@(werte[imShift   +   i*sStep])),pStep)
        else if not senkrecht and (algo.inOrdnung = doAlleResSmi) then      // "Smi" bezieht sich auf die Dimension, in der nicht transformiert wird!
          algo.laden(invers,pExtended(@(werte[i*sStep])),pExtended(@(werte[2*imShift-(1+i)*sStep])),pStep)
        else
          algo.laden(invers,pExtended(@(werte[i*sStep])),pStep);
        algo.summen(true,vorher,in0);
        algo.ausfuehren;
        algo.summen(false,nachher,out0);
        if not senkrecht and (algo.outOrdnung = doAlleResIms) then
          algo.speichern(invers,pExtended(@(werte[i*sStep])),pExtended(@(werte[imShift   +   i*sStep])),pStep)
        else if not senkrecht and (algo.outOrdnung = doAlleResSmi) then     // "Smi" bezieht sich auf die Dimension, in der nicht transformiert wird!
          algo.speichern(invers,pExtended(@(werte[i*sStep])),pExtended(@(werte[2*imShift-(1+i)*sStep])),pStep)
        else
          algo.speichern(invers,pExtended(@(werte[i*sStep])),pStep);
      end;
  end{of case};

  if (nachher=0) and (vorher=0) then pvFehler:=0
  else pvFehler:=abs(nachher-vorher)/(nachher+vorher);

  if invers and (length(hg)<>0) then begin // Hintergrund addieren
    if length(hg)=1 then
      for i:=sMin to sMax do
        for j:=0 to pMax do
          werte[i*sStep+j*pStep]:=
            werte[i*sStep+j*pStep]+hg[0]
    else if length(hg)<>0 then begin
      if length(hg)<>params.xSteps then begin
        gibAus('fft: Der zu addierende Hintergrund hat die falsche Länge: '+intToStr(length(hg))+' vs. '+intToStr(params.xSteps),1);
        exit;
      end;
      if senkrecht then
        for j:=0 to pMax do
          for i:=sMin to sMax do
            werte[i+j*params.xSteps]:=werte[i+j*params.xSteps] + hg[i]
      else
        for j:=sMin to sMax do
          for i:=0 to pMax do
            werte[i+j*params.xSteps]:=werte[i+j*params.xSteps] + hg[i];
    end;
  end;

  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);
  result:=not (in0 or out0);
end;

{$DEFINE tLLWerte_radonTransformation}
procedure tLLWerte.radonTransformation(xMin,xMax: longint; xStep,yStep: extended; qu: pTLLWerteSingle);
{$INCLUDE werteunit.inc}
procedure tLLWerte.radonTransformation(xMin,xMax: longint; xStep,yStep: extended; qu: pTLLWerteDouble);
{$INCLUDE werteunit.inc}
procedure tLLWerte.radonTransformation(xMin,xMax: longint; xStep,yStep: extended; qu: pTLLWerteExtended);
{$INCLUDE werteunit.inc}
{$UNDEF tLLWerte_radonTransformation}

procedure tLLWerte.schreibeWert(var f: textfile; x,y: extended; beschriftung: tExtPoint; var letzterWert: extended; entspringen,verschiebung: extended; skalierung: string; kvs: tKnownValues; linienIntegral: tLinienIntegral; cbgv: tCallBackGetValue);
var
  tmp: extended;
  i:   longint;
begin
  if linienIntegral.schritte=0 then
    tmp:=werte[round(x)+round(y)*params.xSteps]+verschiebung
  else begin
    tmp:=verschiebung;
    for i:=0 to linienIntegral.schritte-1 do
      tmp:=
        tmp+
        werte[
          round(x+linienIntegral.von['x']+i*linienIntegral.schritt['x'])+
          round(y+linienIntegral.von['y']+i*linienIntegral.schritt['y'])*params.xSteps];
  end;
  if isNaN(letzterWert) or (entspringen<0) then
    letzterWert:=tmp;
  if entspringen>0 then begin
    while tmp>=letzterWert+entspringen/2 do
      tmp:=tmp-entspringen;
    while tmp<letzterWert-entspringen/2 do
      tmp:=tmp+entspringen;
  end;
  letzterWert:=tmp;
  schreibeWert(f,beschriftung['x'],beschriftung['y'],tmp,skalierung,kvs,cbgv);
end;

procedure tLLWerte.schreibeWert(var f: textfile; xBeschr,yBeschr,wert: extended; skalierung: string; kvs: tKnownValues; cbgv: tCallBackGetValue);
begin
  if skalierung<>'1' then
    wert:=wert*exprToFloat(false,skalierung,kvs,cbgv);
  writeln(f,myFloatToStr(xBeschr)+' '+myFloatToStr(yBeschr)+' '+myFloatToStr(wert));
end;

procedure tLLWerte.schreibeWertIntegriert(var f: textfile; i: longint; hor: boolean; var letzterWert: extended; entspringen: extended);
var
  j:   longint;
  tmp: extended;
begin
  tmp:=0;
  if hor then begin
    for j:=0 to params.xSteps-1 do
      tmp:=tmp+werte[j+i*params.xSteps];
  end
  else begin
    for j:=0 to params.tSiz-1 do
      tmp:=tmp+werte[i+j*params.xSteps];
  end;
  if isNaN(letzterWert) or (entspringen<0) then
    letzterWert:=tmp;
  if entspringen>0 then begin
    while tmp>=letzterWert + entspringen/2 do
      tmp:=tmp-entspringen;
    while tmp<letzterWert - entspringen/2 do
      tmp:=tmp+entspringen;
  end;
  letzterWert:=tmp;
  if hor then
    schreibeWert(f,(params.xSteps-1)/2,i,tmp,'1',nil,nil)
  else
    schreibeWert(f,i,(params.tSiz-1)/2,tmp,'1',nil,nil);
end;

procedure tLLWerte.erzeugeBinning(senkrecht,linien: boolean; x0,dX: extended; s: string);
var
  f:     textfile;
  i:     longint;
  sum,x: extended;
begin
  assignFile(f,s);
  rewrite(f);
  sum:=0;
  while x0<0 do
    x0:=x0+dX;
  for i:=0 to params.xSteps*params.tSiz-1 do
    if i+1>x0 then begin
      sum:=sum+werte[i]*(x0-i);
      if senkrecht then
        x:=x0/(params.tSiz-1)*(params.tStop-params.tStart)+params.tStart
      else
        x:=x0/(params.xSteps-1)*(params.xStop-params.xStart)+params.xStart;
      writeln(f,floatToStr(x)+' '+floatToStr(sum/dX));
      if linien then begin
        writeln(f,floatToStr(x)+' 0');
        writeln(f)
      end;
      sum:=werte[i]*(i+1-x0);
      x0:=x0+dX;
    end
    else
      sum:=sum+werte[i];
  closeFile(f);
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,j: longint;
begin
  case nB of
    doResIms,doResSmi: ;
    doBetr: begin
      // die vier "Ecken" sind rein reell
      for i:=0 to 1 do
        for j:=0 to 1 do
          werte[j*(params.tSiz div 2)*params.xSteps + i*params.xSteps div 2]:=
            abs(extended(werte[j*(params.tSiz div 2)*params.xSteps + i*params.xSteps div 2]));

      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;
      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
      // die vier "Ecken" sind rein reell
      for i:=0 to 1 do
        for j:=0 to 1 do
          werte[j*(params.tSiz div 2)*params.xSteps + i*params.xSteps div 2]:=
            sqr(extended(werte[j*(params.tSiz div 2)*params.xSteps + i*params.xSteps div 2]));

      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;
      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;
    else
      fehler('Nachbearbeitungsmethode '+fftDoToStr(nB)+' nicht implementiert!');
  end{of case};
end;

procedure tLLWerte.fft2dNachbearbeitungB(xMin,xMax: longint; nB: tFFTDatenordnung);
var
  i,j:                 longint;
  reRe,reIm,imRe,imIm: extended;
begin // bearbeitet nur den Hauptteil (außer erster und mittlerer Zeile/Spalte) nach!
  case nB of
    doResIms,doResSmi: ;
    doBetr: begin
      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]:=
            sqrt(
              sqr(reRe-imIm)  // Re^2
             +sqr(reIm+imRe)  // Im^2
            );
          werte[params.xSteps-i+(params.tSiz-j)*params.xSteps]:=
            werte[i+j*params.xSteps];
          werte[params.xSteps-i+j*params.xSteps]:=
            sqrt(
              sqr(reRe+imIm)  // Re^2
             +sqr(reIm-imRe)  // Im^2
            );
          werte[i+(params.tSiz-j)*params.xSteps]:=
            werte[params.xSteps-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
          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]:=
              sqr(reRe-imIm)  // Re^2
             +sqr(reIm+imRe); // Im^2
          werte[params.xSteps-i+(params.tSiz-j)*params.xSteps]:=
            werte[i+j*params.xSteps];
          werte[params.xSteps-i+j*params.xSteps]:=
              sqr(reRe+imIm)  // Re^2
             +sqr(reIm-imRe); // Im^2
          werte[i+(params.tSiz-j)*params.xSteps]:=
            werte[params.xSteps-i+j*params.xSteps];
        end;
    end;
    else
      fehler('Nachbearbeitungsmethode '+fftDoToStr(nB)+' nicht implementiert!');
  end{of case}
end;

procedure tLLWerte.fft2dNachbearbeitungVerdoppeln(nB: tFFTDatenordnung);
var
  i,j,xS2,tS4:         int64;
  reRe,reIm,imRe,imIm: extended;
begin
  // Vorsicht:
  //   Spalten sind zuerst tSiz/2 lang und dann tSiz.
  //   Zeilen behalten ihre Länge von xSteps.

  // Layout vorher: (doResSmi in x und t)
  //
  // reIm | imIm
  // -----+-----
  // reRe | imRe
  //

  xS2:=params.xSteps div 2;
  tS4:=params.tSiz div 4;
  case nB of
    doAlleResSmi:
      for j:=tS4 downto 0 do
        for i:=xS2 downto 0 do begin
          reRe:=werte[i + j*2*xS2];
          if (i<>0) and (i<>xS2) then
            imRe:=werte[-i + (j+1)*2*xS2]
          else
            imRe:=0;            // erste und mittlere Spalte ist x-reell
          if (j<>0) and (j<>tS4) then
            reIm:=werte[i + (2*tS4-j)*2*xS2]
          else
            reIm:=0;            // erste und mittlere Zeile ist t-reell
          if (i<>0) and (i<>xS2) and (j<>0) and (j<>tS4) then
            imIm:=werte[-i + (2*tS4-j+1)*2*xS2]
          else
            imIm:=0;            // erste und mittlere Zeile und Spalte sind x-t-reell

          werte[i     +          j *2*xS2]:=reRe-imIm;
          werte[i     + (4*tS4-1-j)*2*xS2]:=reIm+imRe;
          if (i<>0) and (i<>xS2) then begin
            werte[2*xS2-i +          j *2*xS2]:= reRe+imIm;
            werte[2*xS2-i + (4*tS4-1-j)*2*xS2]:=-imRe+reIm;
          end;
          if (j<>0) and (j<>tS4) then begin
            werte[i     + (2*tS4-1-j)*2*xS2]:=reRe+imIm;
            werte[i     +   (2*tS4+j)*2*xS2]:=imRe-reIm;
            if (i<>0) and (i<>xS2) then begin
              werte[2*xS2-i + (2*tS4-1-j)*2*xS2]:= reRe-imIm;
              werte[2*xS2-i +   (2*tS4+j)*2*xS2]:=-imRe-reIm;
            end;
          end;
        end;
    doAlleResIms:
      for j:=tS4 downto 0 do
        for i:=xS2 downto 0 do begin
          reRe:=werte[i + j*2*xS2];
          if (i<>0) and (i<>xS2) then
            imRe:=werte[-i + (j+1)*2*xS2]
          else
            imRe:=0;            // erste und mittlere Spalte ist x-reell
          if (j<>0) and (j<>tS4) then
            reIm:=werte[i + (2*tS4-j)*2*xS2]
          else
            reIm:=0;            // erste und mittlere Zeile ist t-reell
          if (i<>0) and (i<>xS2) and (j<>0) and (j<>tS4) then
            imIm:=werte[-i + (2*tS4-j+1)*2*xS2]
          else
            imIm:=0;            // erste und mittlere Zeile und Spalte sind x-t-reell

          werte[i     +        j *2*xS2]:=reRe-imIm;
          werte[i     + (2*tS4+j)*2*xS2]:=reIm+imRe;
          if (i<>0) and (i<>xS2) then begin
            werte[2*xS2-i +        j *2*xS2]:= reRe+imIm;
            werte[2*xS2-i + (2*tS4+j)*2*xS2]:=-imRe+reIm;
          end;
          if (j<>0) and (j<>tS4) then begin
            werte[i     + (2*tS4-j)*2*xS2]:=reRe+imIm;
            werte[i     + (4*tS4-j)*2*xS2]:=imRe-reIm;
            if (i<>0) and (i<>xS2) then begin
              werte[2*xS2-i + (2*tS4-j)*2*xS2]:= reRe-imIm;
              werte[2*xS2-i + (4*tS4-j)*2*xS2]:=-imRe-reIm;
            end;
          end;
        end;
    else
      fehler('Nachbearbeitungsmethode '+fftDoToStr(nB)+' verdoppelt die Anzahl der Daten nicht!');
  end{of case};
end;

procedure tLLWerte.fft2dNachbearbeitungKomplex(xMin,xMax: longint; nB: tFFTDatenordnung);
var
  i,j,hLen: longint;
begin
  hLen:=params.tSiz div 2;
  if not params.istKomplex then
    fehler('Komplexe FFT2dNachbearbeitung geht nur mit vollkomplexen Werten!');
  case nB of
    doResIms,doResSmi,doRes:
      fehler('Vollkomplexe Werte können nicht auf die Hälfte komprimiert werden!');
    doBetr:
      for j:=0 to hLen-1 do
        for i:=xMin to xMax do
          werte[i+j*params.xSteps]:=
            sqrt(
              sqr(extended(werte[i+j*params.xSteps]))+
              sqr(extended(werte[i+(j+hLen)*params.xSteps]))
            );
    doBetrQdr:
      for j:=0 to hLen-1 do
        for i:=xMin to xMax do
          werte[i+j*params.xSteps]:=
            sqr(extended(werte[i+j*params.xSteps]))+
            sqr(extended(werte[i+(j+hLen)*params.xSteps]));
    doAlleResIms: ;
    doAlleResSmi:
      fehler('Nachbearbeitungsmethode '+fftDoToStr(nB)+' ist aus Faulheit noch nicht implementiert -- mach''s doch einfach!');
    doGetrennt:
      fehler('fft2dNachbearbeitungKomplex kann Werte nicht getrennt speichern!');
    else
      fehler('fft2dNachbearbeitungKomplex kennt Nachbearbeitungsmethode '+fftDoToStr(nB)+' noch nicht!');
  end{of case};
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;
  if (ausgaben and __ausgabenMaske) <> 0 then
    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);
  if (ausgaben and __ausgabenMaske) <> 0 then
    gibAus('... fertig '+timetostr(now-Zeit),ausgaben);
end;

function tLLWerte.zuPixelWerten(wHoehe,wBreite,xPMi,xMi,tMi: longint; xZ,yZ: extended; mo: boolean; pPWerte: pTExtendedArray; pPAnzahlen: pTLongintArray): boolean;
var
  i,j,k,l,
    xV,xB,tV,tB: longint;
  hLen:          int64;
  b,imPart:      boolean;
  wert:          extended;
begin
  result:=false;
  for i:=0 to length(pPWerte^)-1 do
    pPWerte^[i]:=0;
  for i:=0 to length(pPAnzahlen^)-1 do
    pPAnzahlen^[i]:=0;
  b:=false;
  hLen:=params.tSiz div (1+byte(params.istKomplex));
  for imPart:=false to params.istKomplex do
    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(hLen-1,max(0,ceil((j-1/2)/yZ+tMi)))+byte(imPart)*hLen;
        tB:=min(hLen-1,max(0,ceil((j+1/2)/yZ+tMi-1)))+byte(imPart)*hLen;
        if xV>xB then begin
          if (i>0) or (xPMi>0) then
            dec(xV)
          else
            inc(xB);
        end;
        if tV>tB then begin
          if j>0 then
            dec(tV)
          else
            inc(tB);
        end;
        if (xV>xB) or (tV>tB) then begin
          gibAus('Keine Inputwerte für Position '+intToStr(i)+':'+intToStr(j)+'!',1);
          exit;
        end;
        if not imPart then
          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);
            wert:=werte[k+l*params.xSteps];
            if mo then begin
              while wert > params.maxW do
                wert:=wert - params.maxW + params.minW;
              while wert < params.minW do
                wert:=wert + params.maxW - params.minW;
            end
            else
              wert:=min(params.maxW,max(params.minW,wert));
            pPWerte^[(j+wHoehe*byte(imPart))*wBreite+i]:=pPWerte^[(j+wHoehe*byte(imPart))*wBreite+i]+wert;
          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)+' x '+intToStr(tMi)+'-'+intToStr(tMa)+') '+floatToStr(Schw)+' ...',1);
  i:=0;
  if params.xSteps>1 then
    dX:=(params.xStop-params.xStart)/(params.xSteps-1)
  else
    dX:=0;
  x0:=params.xStart;
  if params.tSiz>1 then
    dY:=(params.tStop-params.tStart)/(params.tSiz-1)
  else
    dY:=0;
  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;

{$DEFINE tLLWerte_integriere}
procedure tLLWerte.integriere(qu: pTLLWerteSingle; xMi,xMa,tMi,tMa,xOf,tOf: longint; richtung: tIntegrationsRichtung);
{$INCLUDE werteunit.inc}
procedure tLLWerte.integriere(qu: pTLLWerteDouble; xMi,xMa,tMi,tMa,xOf,tOf: longint; richtung: tIntegrationsRichtung);
{$INCLUDE werteunit.inc}
procedure tLLWerte.integriere(qu: pTLLWerteExtended; xMi,xMa,tMi,tMa,xOf,tOf: longint; richtung: tIntegrationsRichtung);
{$INCLUDE werteunit.inc}
{$UNDEF tLLWerte_integriere}

procedure tLLWerte.gauszFit(amplituden,breiten,positionen,ueberlappe,hintergruende: pTLLWerteExtended; von,bis: longint; senkrecht: boolean; fensterBreite,maxBreite,maxVerschiebung: extended; positionsMitten: tExtendedArray);
var
  i,j,ii,zaehl,                    // Zähler
    qpSchritt,qsSchritt,           // Schrittlängen in der Quelle
    zpSchritt,zdSchritt,           // Schrittlängen im Ziel
    pLen,offset,                   // Datenlänge und Offset in der Quelle
    wiMin,wiMax,                   // momentane Fensterränder (für den Zählindex) in der Quelle
    verbesserung,                  // <0: schlechter, 1: besser, 2: viel besser
    aParams,nParams:     longint;  // Index der alten und neuen Parameter
  parameter:             array[0..1,boolean,0..3] of extended;
    // dim0 nummeriert Parametersätze
    // dim1: Ableitung (true) oder werte (false)
    // dim2: position, 1/breite, amplitude, hintergrund
  fehlers:               array[0..1] of extended;
  t0,t1,t2,t3,                     // Zwischenergebnisse
    qdrSumm,                       // Quadratesumme im betrachteten Fenster
    pMi,pMa,                       // Achsenenden in Datenrichtung
    schrittFaktor:       extended;
  ignBr,ignVersch:       boolean;  // Breite/Verschiebung am Anschlag!
begin
  if senkrecht then begin
    qpSchritt:=params.xSteps; // Schritt in der Quelle parallel zur Fit-Richtung
    qsSchritt:=1;             // Schritt in der Quelle senkrecht zur Fit-Richtung
    zpSchritt:=params.xSteps; // Schritt im Ziel in positionsMitten-Richtung
    zdSchritt:=1;             // Schritt im Ziel in Daten-Richtung (= senkrecht zur Fit-Richtung)
    pMi:=params.tStart;
    pMa:=params.tStop;
    pLen:=params.tSiz;
  end
  else begin
    qsSchritt:=params.xSteps;
    qpSchritt:=1;
    zdSchritt:=length(positionsMitten);
    zpSchritt:=1;
    pMi:=params.xStart;
    pMa:=params.xStop;
    pLen:=params.xSteps;
  end;
  maxBreite:=1/maxBreite;

  for i:=von to bis do
    for j:=0 to length(positionsMitten)-1 do begin
      offset:=i*qsSchritt;
      aParams:=0;

      wiMin:=0;
      wiMax:=pLen-1;
      if fensterBreite>0 then begin
        wiMin:=max(wiMin,round(positionsMitten[j] - fensterBreite / 2));
        wiMax:=min(wiMax,round(positionsMitten[j] + fensterBreite / 2));
      end;

      qdrSumm:=0;
      t0:=0;
      t1:=0;
      t2:=0;
      schrittFaktor:=0;
      for ii:=wiMin to wiMax do begin
        t3:=werte[offset+qpSchritt*ii];
        t3:=abs(t3);
        qdrSumm:=qdrSumm + sqr(t3);
        t0:=t0 + t3;
        t1:=t1 + ii*t3;
        t2:=t2 + sqr(ii)*t3;
        if t3>schrittFaktor then schrittFaktor:=t3;
      end;

      t1:=t1/t0;
      t2:=t2/t0;

      parameter[aParams,false,0]:=t1;                     // Erwartungswert
      parameter[aParams,false,2]:=schrittFaktor;          // Maximalwert
      parameter[aParams,false,1]:=parameter[aParams,false,2]/t0*sqrt(pi);
//      parameter[aParams,false,1]:=1/sqrt(2*(t2-sqr(t1))); // beinhalted Standardabweichung  <-- schlechter!!
      parameter[aParams,false,3]:=0;

      if (maxBreite>=0) and (parameter[aParams,false,1]<2*maxBreite) then
        parameter[aParams,false,1]:=maxBreite*2;

      nParams:=aParams;
      {$DEFINE gauszFitBerechneWerte}
      {$I gauszFit.inc} // werte + Gradienten berechnen
      {$UNDEF}

      schrittFaktor:=1;
      if maxVerschiebung>=0 then
        schrittFaktor:=min(schrittFaktor,maxVerschiebung/max(abs(parameter[nParams,true,0]),1e-50));
      if maxBreite>=0 then
        schrittFaktor:=min(schrittFaktor,abs(parameter[nParams,false,2]-maxBreite)/max(abs(parameter[nParams,true,1]),1e-50));

      zaehl:=0;
      repeat
        nParams:=(aParams+1) mod length(parameter);
        for ii:=0 to length(parameter[aParams,false])-1 do
          parameter[nParams,false,ii]:=
            parameter[aParams,false,ii] - schrittFaktor * parameter[aParams,true,ii]; // * byte(ii<3);

        ignVersch:=false;
        if maxVerschiebung>0 then begin
          if parameter[nParams,false,0]-positionsMitten[j] > maxVerschiebung then begin
            parameter[nParams,false,0]:=positionsMitten[j] + maxVerschiebung;
            ignVersch:=true;
          end;
          if positionsMitten[j]-parameter[nParams,false,0] > maxVerschiebung then begin
            parameter[nParams,false,0]:=positionsMitten[j] - maxVerschiebung;
            ignVersch:=true;
          end;
        end;

        ignBr:=false;
        if maxBreite>0 then
          if parameter[nParams,false,1] < maxBreite then begin
            parameter[nParams,false,1]:=maxBreite;
            ignBr:=true;
          end;

        {$DEFINE gauszFitBerechneWerte}
        {$I gauszFit.inc} // werte + Gradienten berechnen
        {$UNDEF}

        if fehlers[aParams]>2*fehlers[nParams] then begin
          verbesserung:=2;
          schrittFaktor:=schrittFaktor * 2;
          aParams:=nParams;
        end
        else if fehlers[aParams]>fehlers[nParams] then begin
          verbesserung:=1;
          schrittFaktor:=schrittFaktor * 1.3;
          aParams:=nParams;
        end
        else begin
          if verbesserung>0 then
            verbesserung:=-1
          else
            dec(verbesserung);
          schrittFaktor:=schrittFaktor * 0.5;
        end;

        if schrittFaktor<1e-50 then
          fehler('Sehr kleiner Schrittfaktor ('+floatToStr(schrittFaktor)+') beim Fitten!');

        inc(zaehl);
        {$IFDEF gauszFitStatus}
        if (zaehl mod 10000)=0 then
          gibAus(
            floatToStr(fehlers[aParams])+' '+
            floatToStr(qdrSumm)+' '+
            intToStr(byte((verbesserung<-10) or (fehlers[aParams]*100<qdrSumm)))+' '+
            intToStr(byte(ignVersch or (abs(schrittFaktor*parameter[aParams,true,0])<1)))+' '+
            intToStr(byte(ignBr or (abs(schrittFaktor*parameter[aParams,true,1]/parameter[aParams,false,1])<0.01)))+' '+
            intToStr(byte(abs(schrittFaktor*parameter[aParams,true,2]/parameter[aParams,false,2])<0.01))+' '+
            intToStr(byte(abs(schrittFaktor*parameter[aParams,true,3]/max(abs(parameter[aParams,false,3]),1e-10))<0.01))+' '+
            intToStr(verbesserung),
            3);
        {$ENDIF}

      until ((zaehl>=100000) or (verbesserung<-10) or (fehlers[aParams]*10<qdrSumm)) and
            (ignVersch or (abs(schrittFaktor*parameter[aParams,true,0])<1)) and
            (ignBr or (abs(schrittFaktor*parameter[aParams,true,1]/parameter[aParams,false,1])<0.01)) and
            (abs(schrittFaktor*parameter[aParams,true,2]/parameter[aParams,false,2])<0.01) and
            (abs(schrittFaktor*parameter[aParams,true,3]/max(abs(parameter[aParams,false,3]),1e-10))<0.01);
      {$IFDEF gauszFitStatus}
      gibAus('',3);
      {$ENDIF}

      if assigned(positionen) then
        positionen^.werte[j*zpSchritt + i*zdSchritt]:=pMi+(pMa-pMi)/(pLen-1)*parameter[aParams,false,0];
      if assigned(breiten) then
        breiten^.werte[j*zpSchritt + i*zdSchritt]:=(pMa-pMi)/(pLen-1)/parameter[aParams,false,1];
      if assigned(amplituden) then
        amplituden^.werte[j*zpSchritt + i*zdSchritt]:=parameter[aParams,false,2];
      if assigned(hintergruende) then
        hintergruende^.werte[j*zpSchritt + i*zdSchritt]:=parameter[aParams,false,3];
      if assigned(ueberlappe) then
        ueberlappe^.werte[j*zpSchritt + i*zdSchritt]:=fehlers[aParams]/qdrSumm;
    end;
end;

function tLLWerte.ermittleRandDurchschnitt: 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;

function tLLWerte.ermittleRandMinimum: extended;
var
  i: int64;
begin
  result:=werte[0];
  for i:=0 to params.xSteps-1 do
    if werte[i]<result then
      result:=werte[i];
  for i:=1 to params.tSiz-2 do
    if werte[i*params.xSteps]<result then
      result:=werte[i*params.xSteps];
  for i:=0 to params.xSteps-1 do
    if werte[i + (params.tSiz-1)*params.xSteps]<result then
      result:=werte[i + (params.tSiz-1)*params.xSteps];
  for i:=1 to params.tSiz-2 do
    if werte[params.xSteps-1 + i*params.xSteps]<result then
      result:=werte[params.xSteps-1 + i*params.xSteps];
end;

function tLLWerte.ermittleRandPerzentil(p: extended): extended;
var
  i:    int64;
  posi: tInt64PointArray;
begin
  setLength(posi,2*(params.xSteps+params.tSiz-2));
  for i:=0 to length(posi)-1 do
    posi[i]['x']:=0;
  for i:=0 to params.xSteps-1 do begin
    posi[i]:=
      int64Point(i,0);
    posi[i+params.xSteps]:=
      int64Point(i,params.tSiz-1);
  end;
  for i:=1 to params.tSiz-2 do begin
    posi[i-1+2*params.xSteps]:=
      int64Point(0,i);
    posi[i-3+2*params.xSteps+params.tSiz]:=
      int64Point(params.xSteps-1,i);
  end;
  sortiereNachWert(posi,false);
  i:=round(p*(length(posi)-1));
  if i<0 then i:=0;
  if i>=length(posi) then i:=length(posi)-1;
  gibAus(
    floatToStr(extended(werte[posi[0]['x']+posi[0]['y']*params.xSteps])) +
    ' .. ' +
    floatToStr(extended(werte[posi[length(posi)-1]['x']+posi[length(posi)-1]['y']*params.xSteps])) +
    ' -(' +
    floatToStr(p) +
    ')-> ' +
    floatToStr(extended(werte[posi[i]['x']+posi[i]['y']*params.xSteps])),
    1
  );
  result:=werte[posi[i]['x']+posi[i]['y']*params.xSteps];
end;

procedure tLLWerte.integriereVertikal(xMi,xMa,tMi,tMa: int64; hg: pTExtendedArray);
var
  i,j: int64;
begin
  for i:=xMi to xMa do begin
    hg^[i]:=0;
    for j:=tMi to tMa do
      hg^[i]:=hg^[i]+werte[i+j*params.xSteps];
    hg^[i]:=hg^[i]/(tMa+1-tMi);
  end;
end;

procedure tLLWerte.integriereVertikalMitRand(xMi,xMa,tMi,tMa,tRa: int64; hg: pTExtendedArray);
var
  i,j:        int64;
  positionen: array of tInt64Point;
begin
  setLength(positionen,tMa-tMi+1);
  tRa:=min(tRa,(tMa-tMi) div 2);
  for i:=xMi to xMa do begin
    for j:=tMi to tMa do
      positionen[j-tMi]:=int64Point(i,j);
    sortiereNachWert(positionen,false);
    hg^[i]:=0;
    for j:=tRa to length(positionen)-tRa-1 do
      hg^[i]:=hg^[i]+werte[positionen[j]['x']+positionen[j]['y']*params.xSteps];
    hg^[i]:=hg^[i]/(tMa+1-tMi-2*tRa);
  end;
end;

procedure tLLWerte.kantenFilter(betraege: tLLWerte; xFak,yFak: extended; filterTyp: tKantenFilterTyp);
var
  dummy: tInt64Point;
begin
  kantenFilter(betraege,xFak,yFak,filterTyp,false,dummy);
end;

procedure tLLWerte.kantenFilter(betraege: tLLWerte; xFak,yFak: extended; filterTyp: tKantenFilterTyp; einseitig: boolean; out maxPos: tInt64Point);
var
  maxima:                      tInt64PointArray;
  i,iM,j,jM,jP,mCnt,di,dj:     int64;
  wert,minWert,maxWert,radius: extended;
  istVollKomplex:              byte; // 1=nein, 2=ja
begin
  istVollKomplex:=params.tSiz div betraege.params.tSiz;
  if (betraege.params.tSiz*istVollKomplex <> params.tSiz) or
    (betraege.params.xSteps <> params.xSteps) or
    (not (istVollKomplex in [1,2])) then
    fehler(
      'Die Dimension der Beträge ('+intToStr(betraege.params.xSteps)+' x '+intToStr(betraege.params.tSiz)+') stimmt nicht mit '+
      'der Dimension der Werte ('+intToStr(params.xSteps)+' x '+intToStr(params.tSiz)+') -- oder deren ersten t-Hälfte -- überein ('+intToStr(istVollKomplex)+')!');
  if einseitig then begin
    if istVollKomplex<>2 then
      fehler('Kann nur voll komplexe Werte einseitig per Kantenfilter filtern ('+intToStr(istVollKomplex)+')!');
    if filterTyp<>kfHochpass then
      fehler('Kann nur einen Hochpass als einseitigen Kantenfilter verwenden!');
  end;

  setLength(maxima,0);
  mCnt:=0;
  for j:=0 to betraege.params.tSiz-1 do begin
    jM:=j - 1 + betraege.params.tSiz*byte(j=0);
    jP:=j + 1 - betraege.params.tSiz*byte(j=betraege.params.tSiz-1);
    for i:=0 to betraege.params.xSteps div 2 + 1 do begin
      iM:=i - 1 + betraege.params.xSteps*byte(i=0);
      wert:=betraege.werte[i+j*betraege.params.xSteps];
      if
        (wert > betraege.werte[iM+j*betraege.params.xSteps]) and
        (wert > betraege.werte[i+jM*betraege.params.xSteps]) and
        (wert > betraege.werte[(i+1)+j*betraege.params.xSteps]) and
        (wert > betraege.werte[i+jP*betraege.params.xSteps]) and
        (2*i<>params.xSteps) and
        (2*j*istVollKomplex<>params.tSiz) 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);

  betraege.sortiereNachWert(maxima,true);
  maxPos:=maxima[0];
  maxWert:=0;
  for i:=0 to length(maxima)-1 do begin
    minWert:=betraege.abstandsQuadratFuerDominanz(symmetrischModulo(maxima[i],int64Point(betraege.params.xSteps,betraege.params.tSiz)),maxima[i],xFak,yFak);
    for j:=0 to i-1 do begin
      wert:=betraege.dominanzQuadrat(maxima[i],maxima[j],xFak,yFak);
      if wert<minWert then begin
        minWert:=wert;
        if minWert<=maxWert then
          break;
      end;
    end;
    if (i=0) or (minWert>maxWert) then begin
      maxWert:=minWert;
      maxPos:=maxima[i];
    end;
  end;

  if istVollKomplex=1 then begin
    iM:=params.xSteps div 2 + 1;
    jM:=params.tSiz div 2 + 1;
    radius:=sqr(min(maxPos['x'],iM-maxPos['x'])*xFak) + sqr(min(maxPos['y'],jM-maxPos['y'])*yFak);
    for j:=0 to jM do
      for i:=0 to iM do begin
        wert:=sqrt((sqr(i)*xFak+sqr(j)*yFak)/radius);
        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;
        if filterTyp=kfHochpass then
          wert:=1-wert;

        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
  else if einseitig then begin
    iM:=params.xSteps;
    jM:=params.tSiz div 2;
    for j:=0 to jM-1 do
      for i:=0 to iM-1 do begin
        maxWert:=
          sqrt(
            min(
              torusAbstandsQuadrat(
                int64Point(i,j),
                int64Point(0,0),
                extPoint(xFak,yFak),
                istVollKomplex
              ), // minimaler Abstand zu den Ecken
              torusAbstandsQuadrat(
                int64Point(i,j),
                int64Point(params.xSteps,params.tSiz div istVollKomplex) - maxPos,
                extPoint(xFak,yFak),
                istVollKomplex
              ) // minimaler Abstand zum punktgespiegelten Maximum (welches ja unterdrückt werden soll)
            )
          );
        minWert:=
          sqrt(
            torusAbstandsQuadrat(
              int64Point(i,j),
              maxPos,
              extPoint(xFak,yFak),
              istVollKomplex
            ) // minimaler Abstand zum Maximum
          );

        wert:=(minWert-maxWert)/(minWert+maxWert); // klein <=> eher nehmen; 0 <=> 50% nehmen; groß <=> eher nicht nehmen
        if wert <= -0.2 then
          wert:=1
        else if wert >= 0.2 then
          wert:=0
        else
          wert:=0.5*(1 - sin(wert/0.2*pi/2));

        werte[i+j*params.xSteps]:=            // Re
          werte[i+j*params.xSteps]*wert;
        werte[i+(jM+j)*params.xSteps]:=       // Im
          werte[i+(jM+j)*params.xSteps]*wert;
      end;
  end
  else begin
    iM:=params.xSteps;
    jM:=params.tSiz div 2;
    radius:=sqr(min(maxPos['x'],iM-maxPos['x']))*xFak + sqr(min(maxPos['y'],jM-maxPos['y']))*yFak;
    for j:=0 to jM-1 do begin
      dj:=j;
      if jM-dj<dj then
        dj:=jM-dj;
      for i:=0 to iM-1 do begin
        di:=i;
        if iM-di<di then
          di:=iM-di;
        wert:=sqrt((sqr(di)*xFak+sqr(dj)*yFak)/radius);
        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;
        if filterTyp=kfHochpass then
          wert:=1-wert;

        werte[i+j*params.xSteps]:=            // Re
          werte[i+j*params.xSteps]*wert;
        werte[i+(jM+j)*params.xSteps]:=       // Im
          werte[i+(jM+j)*params.xSteps]*wert;
      end;
    end;
  end;
end;

function tLLWerte.torusAbstandsQuadrat(von,nach: tInt64Point; metrik: tExtPoint; istVollKomplex: longint): extended;
var
  tmp: tInt64Point;
  c:   char;
begin
  tmp:=von-nach;
  for c:='x' to 'y' do
    if tmp[c]<0 then
      tmp[c]:=-tmp[c];
  if 2*tmp['x']>params.xSteps then
    tmp['x']:=params.xSteps-tmp['x'];
  if istVollKomplex*2*tmp['y']>params.tSiz then
    tmp['y']:=params.tSiz div istVollKomplex - tmp['y'];
  result:=metrik['x']*sqr(tmp['x'])+metrik['y']*sqr(tmp['y']);
end;

procedure tLLWerte.fenstereWerte(xMi,xMa,tMi,tMa: int64; xFen,tFen: tFenster; hg: tExtendedArray);
var
  i,j:    int64;
  offset: extended;
begin
  offset:=0;
  if length(hg)=1 then
    offset:=hg[0]
  else if length(hg)<>0 then begin
    if length(hg)<>params.xSteps then begin
      gibAus('Der abzuziehende Hintergrund hat die falsche Länge: '+intToStr(length(hg))+' vs. '+intToStr(params.xSteps),1);
      exit;
    end;
    for j:=tMi to tMa do
      for i:=xMi to xMa do
        werte[i+j*params.xSteps]:=werte[i+j*params.xSteps] - hg[i];
  end;

  {$DEFINE tLLWerte_fenstereWerte_fensterMultiplikation}
  case 4*byte(offset<>0) + 2*byte(tFen.aktiv) + byte(xFen.aktiv) of
    1:
      {$DEFINE hatXFenster}
      {$INCLUDE werteunit.inc}
    3:
      {$DEFINE hatTFenster}
      {$INCLUDE werteunit.inc}
    2:
      {$UNDEF hatXFenster}
      {$INCLUDE werteunit.inc}
    6:
      {$DEFINE hatOffset}
      {$INCLUDE werteunit.inc}
    7:
      {$UNDEF hatXFenster}
      {$INCLUDE werteunit.inc}
    5:
      {$UNDEF hatTFenster}
      {$INCLUDE werteunit.inc}
    4:
      {$UNDEF hatXFenster}
      {$INCLUDE werteunit.inc}
      {$UNDEF hatOffset}
  end{of case};
  {$UNDEF tLLWerte_fenstereWerte_fensterMultiplikation}
end;

procedure tLLWerte.verschiebe(richtung: tInt64Point; xV,xB,tV,tB: longint);
var
  xS,tS,x,t,xN,tN,xM,tM: longint;
  imPart:                boolean;
  wert:                  wGen;
begin
  xM:=params.xSteps;
  tM:=params.tSiz div (1+byte(params.istKomplex));
  if (xV<0) or (xB<0) or (tV<0) or (tB<0) or
    (xV>=xM) or (xB>=xM) or (tV>=tM) or (tB>=tM) then
    fehler('Fehler: Das Startrechteck ('+intToStr(xV)+'-'+intToStr(xB)+'x'+intToStr(tV)+'-'+intToStr(tB)+') liegt nicht vollsändig in den Daten ('+intToStr(xM)+'x'+intToStr(tM)+').');
  for imPart:=false to params.istKomplex do
    for xS:=xV to xB do
      for tS:=tV to tB do begin
        wert:=werte[xS+(tS+byte(imPart)*tM)*xM];
        xN:=xS;
        tN:=tS;
        repeat
          x:=xN;
          t:=tN;
          xN:=(x + richtung['x']) mod xM;
          tN:=(t + richtung['y']) mod tM;
          werte[x+(t+byte(imPart)*tM)*xM]:=werte[xN+(tN+byte(imPart)*tM)*xM];
        until (xN=xS) and (tN=tS);
        werte[x+(t+byte(imPart)*tM)*xM]:=wert;
      end;
end;

procedure tLLWerte.ermittlePhasenWinkel(xMi,xMa: longint);
var
  i,j,ts2: longint;
begin
  ts2:=params.tSiz div 2;
  for j:=0 to ts2-1 do
    for i:=xMi to xMa do
      if (werte[i+(j+ts2)*params.xSteps]=0) and (werte[i+j*params.xSteps]=0) then
        werte[i+j*params.xSteps]:=0
      else
        werte[i+j*params.xSteps]:=
          arctan2(
            werte[i+(j+ts2)*params.xSteps],
            werte[i+j*params.xSteps]
          );
end;

procedure tLLWerte.macheKomplex(tMi,tMa: int64; kmm: tKomplexMachModus; mT: tMersenneTwister);
var
  i,j,hLen: int64;
  arg:      extended;
begin
  hLen:=length(werte) div 2;
  case kmm of
    kmmImNull:
      fillChar(werte[hLen+tMi*params.xSteps],(tMa+1-tMi)*params.xSteps*sizeOf(wGen),0);
    kmmReNull: begin
      move(werte[tMi*params.xSteps],werte[hLen+tMi*params.xSteps],(tMa+1-tMi)*params.xSteps*sizeOf(wGen));
      fillChar(werte[tMi*params.xSteps],(tMa+1-tMi)*params.xSteps*sizeOf(wGen),0);
    end;
    kmmPhZuf:
      for j:=tMi to tMa do
        for i:=0 to params.xSteps-1 do begin
          arg:=2*pi*mT.random;
          werte[i+j*params.xSteps+hLen]:=
            werte[i+j*params.xSteps]*sin(arg);
          werte[i+j*params.xSteps]:=
            werte[i+j*params.xSteps]*cos(arg);
        end;
    else
      fehler('Komplexmachmodus nicht implementiert ('+intToStr(integer(kmm))+')!');
  end{of case};
end;

procedure tLLWerte.entspringe(mi,ma: int64; em: tEntspringModus);
var
  i,j,pStep,sStep,pMax,jStart,bwLen: int64;
  vglWert,bwFak:                     extended;
  bisherigeWerte:                    tExtendedArray;
begin
  case em.modus of
    emKein:
      exit;
    emHorizontal: begin
      pStep:=1;
      sStep:=params.xSteps;
      pMax:=params.xSteps-1;
    end;
    emVertikal: begin
      pStep:=params.xSteps;
      sStep:=1;
      pMax:=params.tSiz-1;
    end;
    else
      fehler('Entspringmodus '''+tEntspringModusToStr(em)+''' nicht implementiert!');
  end{of case};
  assert(em.modus in [emHorizontal,emVertikal],'Hoppla, der Entspringmodus '''+tEntspringModusToStr(em)+''' ist doch noch nicht implementiert!');
  bwLen:=max(1,round(em.parameter[0]));
  bwFak:=1/bwLen;
  setLength(bisherigeWerte,bwLen);
  jStart:=max(0,min(pMax,round(em.parameter[1])));
  for i:=mi to ma do begin // senkrecht
    for j:=0 to bwLen-1 do
      bisherigeWerte[j]:=0;
    vglWert:=0;
    for j:=jStart to pMax do begin // parallel, nach rechts
      while werte[i*sStep+j*pStep]-vglWert >= pi do
        werte[i*sStep+j*pStep]:=werte[i*sStep+j*pStep]-2*pi;
      while werte[i*sStep+j*pStep]-vglWert < -pi do
        werte[i*sStep+j*pStep]:=werte[i*sStep+j*pStep]+2*pi;
      vglWert:=vglWert + bwFak * (werte[i*sStep+j*pStep] - bisherigeWerte[j mod bwLen]);
      bisherigeWerte[j mod bwLen]:=werte[i*sStep+j*pStep];
    end;
    vglWert:=0;
    for j:=jStart to jStart+bwLen-1 do begin
      if j <= pMax then
        bisherigeWerte[j mod bwLen]:=werte[i*sStep+j*pStep]
      else
        bisherigeWerte[j mod bwLen]:=0;
      vglWert:=vglWert + bwFak * bisherigeWerte[j mod bwLen];
    end;
    for j:=jStart-1 downto 0 do begin // parallel, nach links
      while werte[i*sStep+j*pStep]-vglWert >= pi do
        werte[i*sStep+j*pStep]:=werte[i*sStep+j*pStep]-2*pi;
      while werte[i*sStep+j*pStep]-vglWert < -pi do
        werte[i*sStep+j*pStep]:=werte[i*sStep+j*pStep]+2*pi;
      vglWert:=vglWert + bwFak * (werte[i*sStep+j*pStep] - bisherigeWerte[j mod bwLen]);
      bisherigeWerte[j mod bwLen]:=werte[i*sStep+j*pStep];
    end;
  end;
end;

procedure tLLWerte.entferneHeiszePixel(schwelle,minusSchwelle,plusSchwelle: extended);
var
  i,j,iM,jM,iP,jP,mCnt:     int64;
  heiszePixel:              tIntPointArray;
  pixelHitze:               tExtendedArray;
  wert,weLi,weOb,weRe,weUn: extended;
begin
  if params.istKomplex then
    fehler('Heiße Pixel können (noch) nicht aus voll komplexen Werten entfernt werden!');
  if minusSchwelle>plusSchwelle then
    fehler('Unterer Schwellwert ('+floatToStr(minusSchwelle)+') liegt über dem oberen Schwellwert ('+floatToStr(plusSchwelle)+')!');
  schwelle:=(schwelle+1)/4;
  setLength(heiszePixel,0);
  setLength(pixelHitze,0);
  mCnt:=0;
  for j:=0 to params.tSiz-1 do
    for i:=0 to params.xSteps-1 do
      if werte[i+j*params.xSteps]>plusSchwelle then
        werte[i+j*params.xSteps]:=plusSchwelle
      else if werte[i+j*params.xSteps]<minusSchwelle then
        werte[i+j*params.xSteps]:=minusSchwelle;

  for j:=0 to params.tSiz-1 do begin
    jM:=j - 1 + params.tSiz*byte(j=0);
    jP:=j + 1 - params.tSiz*byte(j=params.tSiz-1);
    for i:=0 to params.xSteps-1 do begin
      iM:=i - 1 + params.xSteps*byte(i=0);
      iP:=i + 1 - params.xSteps*byte(i=params.xSteps-1);
      wert:=abs(extended(werte[i+j*params.xSteps]));
      weLi:=abs(extended(werte[iM+j*params.xSteps]));
      weUn:=abs(extended(werte[i+jP*params.xSteps]));
      weRe:=abs(extended(werte[iP+j*params.xSteps]));
      weOb:=abs(extended(werte[i+jM*params.xSteps]));
      if
        (wert > weLi) and
        (wert > weOb) and
        (wert > weRe) and
        (wert > weUn) and
        (wert > (weLi+weOb+weRe+weUn)*schwelle) then begin

        if length(heiszePixel)<=mCnt then begin
          setLength(heiszePixel,length(heiszePixel)+1024);
          setLength(pixelHitze,length(pixelHitze)+1024);
        end;
        heiszePixel[mCnt]['x']:=i;
        heiszePixel[mCnt]['y']:=j;
        pixelHitze[mCnt]:=
          (werte[iM+j*params.xSteps] +
           werte[i+jP*params.xSteps] +
           werte[iP+j*params.xSteps] +
           werte[i+jM*params.xSteps])/4;
        inc(mCnt);
      end;
    end;
  end;
  setLength(heiszePixel,mCnt);
  setLength(pixelHitze,mCnt);
  for i:=0 to length(heiszePixel)-1 do
    werte[heiszePixel[i]['x']+heiszePixel[i]['y']*params.xSteps]:=
      pixelHitze[i];
end;

{$DEFINE tLLWerte_quotient}
procedure tLLWerte.quotioent(dend: pTLLWerteSingle; sor: pTLLWerteSingle; xMi,xMa,xOf,tMi,tMa,tOf: int64; eps: extended);
{$INCLUDE werteunit.inc}
procedure tLLWerte.quotioent(dend: pTLLWerteSingle; sor: pTLLWerteDouble; xMi,xMa,xOf,tMi,tMa,tOf: int64; eps: extended);
{$INCLUDE werteunit.inc}
procedure tLLWerte.quotioent(dend: pTLLWerteSingle; sor: pTLLWerteExtended; xMi,xMa,xOf,tMi,tMa,tOf: int64; eps: extended);
{$INCLUDE werteunit.inc}
procedure tLLWerte.quotioent(dend: pTLLWerteDouble; sor: pTLLWerteSingle; xMi,xMa,xOf,tMi,tMa,tOf: int64; eps: extended);
{$INCLUDE werteunit.inc}
procedure tLLWerte.quotioent(dend: pTLLWerteDouble; sor: pTLLWerteDouble; xMi,xMa,xOf,tMi,tMa,tOf: int64; eps: extended);
{$INCLUDE werteunit.inc}
procedure tLLWerte.quotioent(dend: pTLLWerteDouble; sor: pTLLWerteExtended; xMi,xMa,xOf,tMi,tMa,tOf: int64; eps: extended);
{$INCLUDE werteunit.inc}
procedure tLLWerte.quotioent(dend: pTLLWerteExtended; sor: pTLLWerteSingle; xMi,xMa,xOf,tMi,tMa,tOf: int64; eps: extended);
{$INCLUDE werteunit.inc}
procedure tLLWerte.quotioent(dend: pTLLWerteExtended; sor: pTLLWerteDouble; xMi,xMa,xOf,tMi,tMa,tOf: int64; eps: extended);
{$INCLUDE werteunit.inc}
procedure tLLWerte.quotioent(dend: pTLLWerteExtended; sor: pTLLWerteExtended; xMi,xMa,xOf,tMi,tMa,tOf: int64; eps: extended);
{$INCLUDE werteunit.inc}
{$UNDEF tLLWerte_quotient}

{$DEFINE tLLWerte_produkt}
procedure tLLWerte.produkt(f1: pTLLWerteSingle; f2: pTLLWerteSingle; xMi,xMa,xOf,tMi,tMa,tOf: int64; konj: boolean; daO: tFFTDatenordnung);
{$INCLUDE werteunit.inc}
procedure tLLWerte.produkt(f1: pTLLWerteSingle; f2: pTLLWerteDouble; xMi,xMa,xOf,tMi,tMa,tOf: int64; konj: boolean; daO: tFFTDatenordnung);
{$INCLUDE werteunit.inc}
procedure tLLWerte.produkt(f1: pTLLWerteSingle; f2: pTLLWerteExtended; xMi,xMa,xOf,tMi,tMa,tOf: int64; konj: boolean; daO: tFFTDatenordnung);
{$INCLUDE werteunit.inc}
procedure tLLWerte.produkt(f1: pTLLWerteDouble; f2: pTLLWerteSingle; xMi,xMa,xOf,tMi,tMa,tOf: int64; konj: boolean; daO: tFFTDatenordnung);
{$INCLUDE werteunit.inc}
procedure tLLWerte.produkt(f1: pTLLWerteDouble; f2: pTLLWerteDouble; xMi,xMa,xOf,tMi,tMa,tOf: int64; konj: boolean; daO: tFFTDatenordnung);
{$INCLUDE werteunit.inc}
procedure tLLWerte.produkt(f1: pTLLWerteDouble; f2: pTLLWerteExtended; xMi,xMa,xOf,tMi,tMa,tOf: int64; konj: boolean; daO: tFFTDatenordnung);
{$INCLUDE werteunit.inc}
procedure tLLWerte.produkt(f1: pTLLWerteExtended; f2: pTLLWerteSingle; xMi,xMa,xOf,tMi,tMa,tOf: int64; konj: boolean; daO: tFFTDatenordnung);
{$INCLUDE werteunit.inc}
procedure tLLWerte.produkt(f1: pTLLWerteExtended; f2: pTLLWerteDouble; xMi,xMa,xOf,tMi,tMa,tOf: int64; konj: boolean; daO: tFFTDatenordnung);
{$INCLUDE werteunit.inc}
procedure tLLWerte.produkt(f1: pTLLWerteExtended; f2: pTLLWerteExtended; xMi,xMa,xOf,tMi,tMa,tOf: int64; konj: boolean; daO: tFFTDatenordnung);
{$INCLUDE werteunit.inc}
{$UNDEF tLLWerte_produkt}

procedure tLLWerte.wertAusUmgebungMitteln(x,y: longint);
begin
  werte[x+y*params.xSteps]:=(
    werte[x-1+y*params.xSteps] +
    werte[x+(y-1)*params.xSteps] +
    werte[x+1+y*params.xSteps] +
    werte[x+(y+1)*params.xSteps]
  )/4;
end;

procedure tLLWerte.extrahiereKanten(xMi,xMa,tMi,tMa: longint; vert: boolean; expo: int64);
var
  i,j: int64;
begin
  for j:=tMi to tMa-byte(vert) do
    for i:=xMi to xMa-byte(not vert) do
      werte[i+j*params.xSteps]:=intPower(
        werte[i+byte(not vert)+(j+byte(vert))*params.xSteps]-werte[i+j*params.xSteps],
        expo
      );
  if vert then begin
    for i:=xMi to xMa do
      werte[i+tMa*params.xSteps]:=0;
  end
  else
  for j:=tMi to tMa do
    werte[xMa+j*params.xSteps]:=0;
end;

procedure tLLWerte.nullenEinfuegen(hor,vert: boolean);
var
  i,siz2: int64;
begin
  // Der Speicher ist schon geholt, die Werte stehen aber noch an alter Stelle.
  // Es muss also (ggf.) noch umkopiert und Nullen eingefügt werden.

  if hor then begin
    siz2:=params.xSteps div 2;
    // im Falle von "vert" sparen wir uns die Hälfte, weil das eh nur Nullen sind
    // außerdem muss das nullte Element nicht verschoben werden
    for i:=params.tSiz div (1+byte(vert)) - 1 downto 1 do
      move(werte[i*siz2],werte[i*2*siz2],sizeOf(wGen)*siz2);
    for i:=0 to params.tSiz-1 do
      fillChar(werte[(2*i+1)*siz2],sizeOf(wGen)*siz2,0);
  end;
  if vert then
    fillChar(werte[params.xSteps*(params.tSiz div 2)],sizeOf(wGen)*params.xSteps*(params.tSiz div 2),0);
end;

procedure tLLWerte.skaliere(tMi,tMa: int64; skalierung: string; trafo: tTransformation; kvs: tKnownValues; cbgv: tCallBackGetValue);
var
  knownValues: tKnownValues;
  i,j:         int64;
begin
  knownValues:=tKnownValues.create(kvs);
  for j:=tMi to tMa do begin
    knownValues.add('y',trafo.positionAufAchseZuWert(lLinks,j/(params.tSiz-1),false));
    for i:=0 to params.xSteps-1 do begin
      knownValues.add('x',trafo.positionAufAchseZuWert(lOben,i/(params.xSteps-1),false));
      werte[i+j*params.xSteps]:=
        werte[i+j*params.xSteps]*exprToFloat(false,skalierung,knownValues,cbgv);
    end;
  end;
  knownValues.clear; // sonst werden unsere Werte dem Chef übertragen
  knownValues.free;
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!'#10'Ich kenne:'#10'''Frequenzfenster'''#10'''Sin2''',3);
  result:=false;
end;

function tWavelet.berechneWerte: boolean;
var
  i:          longint;
  tmp:        extended;
  nur0:       boolean;
  tmpFFTAlgo: tFFTAlgorithmus;
begin
  result:=false;
  werte.params.xSteps:=2;
  werte.params.transformationen.xStart:=0;
  werte.params.transformationen.xStop:=1;
  werte.params.transformationen.tStart:=0;
  werte.params.transformationen.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;
        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,nil,pvFehler);
        tmpFFTAlgo.free;
      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;
      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;
  werte.params.refreshKnownValues;
end;

constructor tWavelet.create(globaleWerte: tKnownValues);
var
  ps: tExtraInfos;
begin
  inherited create;
  ps:=tExtraInfos.create(globaleWerte);
  werte:=tLLWerteDouble.create(ps);
  pvFehler:=0;
end;

destructor tWavelet.destroy;
begin
  werte.params.free;
  werte.free;
  inherited destroy;
end;

end.