unit werteunit; {$mode objfpc}{$H+} interface uses Classes, SysUtils, typenunit, math, process, lowlevelunit, matheunit; type // tLLWerte ******************************************************************** pTLLWerteSingle = ^tLLWerteSingle; pTLLWerteDouble = ^tLLWerteDouble; pTLLWerteExtended = ^tLLWerteExtended; generic tLLWerte = class(tObject) { Diese Klasse stellt nur die Berechnungsroutinen und ähnliches bereit, nicht jedoch ein (wie auch immer geartetes) Userinterface. Auch die korrekte Parallelisierung obliegt dem übergeordneten Programmteil. } private public werte: array of wgen; params: 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; destructor destroy; override; function liesDateien(dateien: tGenerischeInputDateiInfoArray): boolean; function fft(senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fen: tFenster; out pvFehler: extended): boolean; overload; function fft(xmin,xmax,tmin,tmax: longint; senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fen: tFenster; out pvFehler: extended): boolean; overload; procedure spiegle; overload; procedure spiegle(tmin,tmax: longint); overload; procedure fft2dNachbearbeitungA(nb: tFFTDatenordnung); procedure fft2dNachbearbeitungB(xmin,xmax: longint; nb: tFFTDatenordnung); procedure schreibeWert(var f: textfile; x,y: longint); overload; procedure schreibeWert(var f: textfile; x,y,wert: extended); overload; procedure schreibeWertIntegriert(var f: textfile; i: longint; hor: boolean); 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; xmin,xmax,tmin,tmax: longint); function zuPixelWerten(whoehe,wbreite,xpmi,xmi,tmi: longint; xz,yz: extended; pPWerte: pTExtendedArray; pPAnzahlen: pTLongintArray): boolean; function findeSchwellwerte(xmi,xma,tmi,tma: longint; Schw: extended): tExtPointArray; procedure integriereSingle(qu: pTLLWerteSingle; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung); procedure integriereDouble(qu: pTLLWerteDouble; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung); procedure integriereExtended(qu: pTLLWerteDouble; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung); end; tLLWerteSingle = specialize tLLWerte; tLLWerteDouble = specialize tLLWerte; tLLWerteExtended = specialize tLLWerte; // andere Typen **************************************************************** tWavelet = class(tObject) freq,tfwhm,pvFehler: extended; hlen: longint; typ: tWaveletTyp; werte: tLLWerteDouble; mitFFT: boolean; constructor create; destructor destroy; override; function setzeTyp(s: string): boolean; function berechneWerte: boolean; end; implementation uses systemunit; constructor tLLWerte.create(ps: tExtrainfos); begin inherited create; params:=ps; setlength(werte,0); end; constructor tLLWerte.create(original: pTLLWerteSingle; ps: tExtrainfos; xmin,xmax: longint); begin inherited create; params:=ps; kopiereVon(false,original,xmin,xmax); end; constructor tLLWerte.create(original: pTLLWerteDouble; ps: tExtrainfos; xmin,xmax: longint); begin inherited create; params:=ps; kopiereVon(false,original,xmin,xmax); end; constructor tLLWerte.create(original: pTLLWerteExtended; ps: tExtrainfos; xmin,xmax: longint); begin inherited create; params:=ps; kopiereVon(false,original,xmin,xmax); end; procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteSingle); begin kopiereVon(st,original,0,original^.params.xsteps-1); end; procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteDouble); begin kopiereVon(st,original,0,original^.params.xsteps-1); end; procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteExtended); begin kopiereVon(st,original,0,original^.params.xsteps-1); end; procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteSingle; xmin,xmax: longint); begin kopiereVon(st,original,xmin,xmax,0,original^.params.tsiz-1); end; procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteDouble; xmin,xmax: longint); begin kopiereVon(st,original,xmin,xmax,0,original^.params.tsiz-1); end; procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteExtended; xmin,xmax: longint); begin kopiereVon(st,original,xmin,xmax,0,original^.params.tsiz-1); end; procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteSingle; xmin,xmax,tmin,tmax: longint); var i,j: longint; begin inherited create; tmax:=min(tmax,original^.params.tsiz-1); tmin:=max(tmin,0); params.tsiz:=tmax+1-tmin; xmax:=min(xmax,original^.params.xsteps-1); xmin:=max(xmin,0); params.xsteps:=xmax+1-xmin; if not params.transformationen.hatNachfolger then params.transformationen.free; params.transformationen:=tKoordinatenAusschnitt.create(original^.params.transformationen,xmin,xmax,tmin,tmax); params.maxW:=0; params.minW:=0; params.np:=original^.params.np; params.beta:=original^.params.beta; params.refreshKnownValues; if not st then begin holeRam(0); for i:=xmin to xmax do for j:=tmin to tmax do werte[i-xmin+(j-tmin)*params.xsteps]:=original^.werte[i+j*original^.params.xsteps]; end; end; procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteDouble; xmin,xmax,tmin,tmax: longint); var i,j: longint; begin inherited create; tmax:=min(tmax,original^.params.tsiz-1); tmin:=max(tmin,0); params.tsiz:=tmax+1-tmin; xmax:=min(xmax,original^.params.xsteps-1); xmin:=max(xmin,0); params.xsteps:=xmax+1-xmin; if not params.transformationen.hatNachfolger then params.transformationen.free; params.transformationen:=tKoordinatenAusschnitt.create(original^.params.transformationen,xmin,xmax,tmin,tmax); params.maxW:=0; params.minW:=0; params.np:=original^.params.np; params.beta:=original^.params.beta; params.refreshKnownValues; if not st then begin holeRam(0); for i:=xmin to xmax do for j:=tmin to tmax do werte[i-xmin+(j-tmin)*params.xsteps]:=original^.werte[i+j*original^.params.xsteps]; end; end; procedure tLLWerte.kopiereVon(st: boolean; original: pTLLWerteExtended; xmin,xmax,tmin,tmax: longint); var i,j: longint; begin inherited create; tmax:=min(tmax,original^.params.tsiz-1); tmin:=max(tmin,0); params.tsiz:=tmax+1-tmin; xmax:=min(xmax,original^.params.xsteps-1); xmin:=max(xmin,0); params.xsteps:=xmax+1-xmin; if not params.transformationen.hatNachfolger then params.transformationen.free; params.transformationen:=tKoordinatenAusschnitt.create(original^.params.transformationen,xmin,xmax,tmin,tmax); params.maxW:=0; params.minW:=0; params.np:=original^.params.np; params.beta:=original^.params.beta; params.refreshKnownValues; if not st then begin holeRam(0); for i:=xmin to xmax do for j:=tmin to tmax do werte[i-xmin+(j-tmin)*params.xsteps]:=original^.werte[i+j*original^.params.xsteps]; end; end; procedure tLLWerte.kopiereVonNach(original: pTLLWerteSingle; qxmin,qxmax,qtmin,qtmax,zxmin,ztmin: longint); var i,j: longint; begin inherited create; for i:=qxmin to qxmax do for j:=qtmin to qtmax do werte[i-qxmin+zxmin + (j-qtmin+ztmin)*params.xsteps]:= original^.werte[i+j*original^.params.xsteps]; end; procedure tLLWerte.kopiereVonNach(original: pTLLWerteDouble; qxmin,qxmax,qtmin,qtmax,zxmin,ztmin: longint); var i,j: longint; begin inherited create; for i:=qxmin to qxmax do for j:=qtmin to qtmax do werte[i-qxmin+zxmin + (j-qtmin+ztmin)*params.xsteps]:= original^.werte[i+j*original^.params.xsteps]; end; procedure tLLWerte.kopiereVonNach(original: pTLLWerteExtended; qxmin,qxmax,qtmin,qtmax,zxmin,ztmin: longint); var i,j: longint; begin inherited create; for i:=qxmin to qxmax do for j:=qtmin to qtmax do werte[i-qxmin+zxmin + (j-qtmin+ztmin)*params.xsteps]:= original^.werte[i+j*original^.params.xsteps]; end; procedure tLLWerte.kopiereVerzerrt(original: pTLLWerteSingle; ZPs: tIntPointArray; ZGs: tExtPointArray; ZAs: tExtendedArray; xmin,xmax,tmin,tmax: longint; vb,nb: tTransformation; va,na: longint); var i,j,k: longint; tmp: extended; begin for i:=tmin to tmax do for j:=xmin to xmax do werte[j+i*params.xsteps]:=0; for i:=0 to length(ZPs)-1 do for j:=0 to 1 do for k:=0 to 1 do if (ZPs[i]['x']+j>=xmin) and (ZPs[i]['x']+j<=xmax) and (ZPs[i]['y']+k>=tmin) and (ZPs[i]['y']+k<=tmax) then begin tmp:=original^.werte[i]; if (va>0) or (na>0) then tmp:=(tmp-original^.params.minW)/(original^.params.maxW-original^.params.minW); if va>0 then vb.transformiereWert(tmp,va-1); tmp:=tmp * (ZGs[i]['x'] * (2*j-1) + 1-j) * (ZGs[i]['y'] * (2*k-1) + 1-k); werte[ZPs[i]['x']+j + (ZPs[i]['y']+k)*params.xsteps]:= werte[ZPs[i]['x']+j + (ZPs[i]['y']+k)*params.xsteps] + tmp; end; for i:=tmin to tmax do for j:=xmin to xmax do begin tmp:=werte[j + i*params.xsteps] / ZAs[j + i*params.xsteps]; if na>0 then tmp:=nb.transformiereWert(tmp,na-1); werte[j + i*params.xsteps]:=tmp; end; end; procedure tLLWerte.kopiereVerzerrt(original: pTLLWerteDouble; ZPs: tIntPointArray; ZGs: tExtPointArray; ZAs: tExtendedArray; xmin,xmax,tmin,tmax: longint; vb,nb: tTransformation; va,na: longint); var i,j,k: longint; tmp: extended; begin for i:=tmin to tmax do for j:=xmin to xmax do werte[j+i*params.xsteps]:=0; for i:=0 to length(ZPs)-1 do for j:=0 to 1 do for k:=0 to 1 do if (ZPs[i]['x']+j>=xmin) and (ZPs[i]['x']+j<=xmax) and (ZPs[i]['y']+k>=tmin) and (ZPs[i]['y']+k<=tmax) then begin tmp:=original^.werte[i]; if (va>0) or (na>0) then tmp:=(tmp-original^.params.minW)/(original^.params.maxW-original^.params.minW); if va>0 then vb.transformiereWert(tmp,va-1); tmp:=tmp * (ZGs[i]['x'] * (2*j-1) + 1-j) * (ZGs[i]['y'] * (2*k-1) + 1-k); werte[ZPs[i]['x']+j + (ZPs[i]['y']+k)*params.xsteps]:= werte[ZPs[i]['x']+j + (ZPs[i]['y']+k)*params.xsteps] + tmp; end; for i:=tmin to tmax do for j:=xmin to xmax do begin tmp:=werte[j + i*params.xsteps] / ZAs[j + i*params.xsteps]; if na>0 then tmp:=nb.transformiereWert(tmp,na-1); werte[j + i*params.xsteps]:=tmp; end; end; procedure tLLWerte.kopiereVerzerrt(original: pTLLWerteExtended; ZPs: tIntPointArray; ZGs: tExtPointArray; ZAs: tExtendedArray; xmin,xmax,tmin,tmax: longint; vb,nb: tTransformation; va,na: longint); var i,j,k: longint; tmp: extended; begin for i:=tmin to tmax do for j:=xmin to xmax do werte[j+i*params.xsteps]:=0; for i:=0 to length(ZPs)-1 do for j:=0 to 1 do for k:=0 to 1 do if (ZPs[i]['x']+j>=xmin) and (ZPs[i]['x']+j<=xmax) and (ZPs[i]['y']+k>=tmin) and (ZPs[i]['y']+k<=tmax) then begin tmp:=original^.werte[i]; if (va>0) or (na>0) then tmp:=(tmp-original^.params.minW)/(original^.params.maxW-original^.params.minW); if va>0 then vb.transformiereWert(tmp,va-1); tmp:=tmp * (ZGs[i]['x'] * (2*j-1) + 1-j) * (ZGs[i]['y'] * (2*k-1) + 1-k); werte[ZPs[i]['x']+j + (ZPs[i]['y']+k)*params.xsteps]:= werte[ZPs[i]['x']+j + (ZPs[i]['y']+k)*params.xsteps] + tmp; end; for i:=tmin to tmax do for j:=xmin to xmax do begin tmp:=werte[j + i*params.xsteps] / ZAs[j + i*params.xsteps]; if na>0 then tmp:=nb.transformiereWert(tmp,na-1); werte[j + i*params.xsteps]:=tmp; end; end; destructor tLLWerte.destroy; begin setlength(werte,0); inherited destroy; end; function tLLWerte.liesDateien(dateien: tGenerischeInputDateiInfoArray): boolean; var i,j,k,l,tmpi,etsiz,spAnz,br: longint; f: file; tmps: single; tmpd: double; tmpe,Zeit: extended; sa: tSingleArray; da: tDoubleArray; ea: tExtendedArray; ipp: tProcess; buf: tByteArray; etwasGelesen: boolean; 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); writeln(tmpi); close(f); exit; end; if dateien[i] is tTraceInputDateiInfo then begin blockread(f,tmpi,sizeof(integer)); // #Traces spAnz:=tmpi; end; blockread(f,etsiz,sizeof(integer)); if dateien[i] is tSpaceTimeInputDateiInfo then begin for j:=0 to etsiz-1 do begin case Dateien[i].Genauigkeit of gSingle: begin blockread(f,tmps,sizeof(single)); // xstart tmpe:=tmps; end; 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) = sizeof(single)) and (Dateien[i].Genauigkeit=gSingle)) or ((sizeof(wgen) = sizeof(double)) and (Dateien[i].Genauigkeit=gDouble)) or ((sizeof(wgen) = sizeof(extended)) and (Dateien[i].Genauigkeit=gExtended)) then blockread(f,werte[(j+Dateien[i].t0abs)*params.xsteps],params.xsteps*sizeof(wgen)) else begin setlength(sa,params.xsteps); blockread(f,sa[0],params.xsteps*sizeof(single)); for k:=0 to params.xsteps-1 do werte[(j+Dateien[i].t0abs)*params.xsteps+k]:=sa[k]; end; if 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 not etwasGelesen then begin gibAus('Ich habe diese Runde keine Daten gelesen!',3); exit; end; end; params.refreshKnownValues; gibAus('... fertig '+timetostr(now-Zeit),1); result:=true; end; procedure tLLWerte.gibMinMaxDichten(out wMi,wMa: extended; xmin,xmax,tmin,tmax: longint); var i,j: longint; begin wMi:=werte[xmin+tmin*params.xsteps]; wMa:=Werte[xmin+tmin*params.xsteps]; for i:=xmin to xmax do for j:=tmin to tmax do begin wMi:=min(wMi,werte[i+j*params.xsteps]); wMa:=max(wMa,werte[i+j*params.xsteps]); end; end; function tLLWerte.fft(senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fen: tFenster; out pvFehler: extended): boolean; var len: longint; begin len:=1; if senkrecht then begin while 2*len<=params.tsiz do len:=len*2; result:=fft(0,params.xsteps-1,(params.tsiz-len) div 2,((params.tsiz+len) div 2) - 1,true,invers,vor,nach,fen,pvFehler); end else begin while 2*len<=params.xsteps do len:=len*2; result:=fft((params.xsteps-len) div 2,((params.xsteps+len) div 2) - 1,0,params.tsiz-1,false,invers,vor,nach,fen,pvFehler); end; end; function tLLWerte.fft(xmin,xmax,tmin,tmax: longint; senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fen: tFenster; out pvFehler: extended): boolean; var i,j,k,n,dist,absch,wnum,wstep,haL, pmax,pmin,smax,smin: longint; in0,out0: boolean; ims,wRe,wIm: tExtendedArray; t1,t2,vorher,nachher,fenavg: extended; umsortierung: tLongintArray; const faktoren: array[tFFTDatenordnung,0..2] of longint = // (doResIms,doResSmi,doRes,doBetr,doBetrQdr); ((2,1,1),(2,1,1),(1,1,1),(1,1,1),(1,1,1)); begin result:=false; if senkrecht then begin pmax:=tmax; pmin:=tmin; smax:=xmax; smin:=xmin; end else begin pmax:=xmax; pmin:=xmin; smax:=tmax; smin:=tmin; end; if length(fen.Werte)<>pmax+1-pmin then fen.berechneWerte(pmax+1-pmin); n:=round(ln(pmax+1-pmin)/ln(2)); if round(power(2,n))<>pmax+1-pmin then begin gibAus(inttostr(pmax+1-pmin)+' ist keine Zweierpotenz, ich bin zu dumm um Werte dieser Anzahl fourierzutransformieren!',3); exit; faktoren[doRes,0]:=faktoren[doRes,1]; // um so eine dämliche Warnung loszuwerden ... end; haL:=(pmax+1-pmin) div 2; if fen.aktiv then begin if senkrecht then begin for i:=xmin to xmax do // Werte fenstern for j:=0 to tmax-tmin do werte[i+(j+tmin)*params.xsteps]:=werte[i+(j+tmin)*params.xsteps]*fen.Werte[j]; end else begin for i:=0 to xmax-xmin do for j:=tmin to tmax do // Werte fenstern werte[i+xmin+j*params.xsteps]:=werte[i+xmin+j*params.xsteps]*fen.Werte[i]; end; fenavg:=0; for i:=0 to pmax-pmin do fenavg:=fenavg+fen.Werte[i]; fenavg:=fenavg/(pmax-pmin+1); if fenavg<0.5 then begin gibAus('Sehr geringer Fensterdurchschnitt: '+floattostr(fenavg)+' ('+inttostr(fen.Breite)+' '+inttostr(length(fen.werte))+' '+inttostr(pmax-pmin+1)+')!',3); end; end; vorher:=0; if vor=doBetrQdr then begin for i:=xmin to xmax do for j:=tmin to tmax do vorher:=vorher + werte[i+j*params.xsteps]; end else begin for i:=xmin to xmax do for j:=tmin to tmax do vorher:=vorher + werte[i+j*params.xsteps]*werte[i+j*params.xsteps] *faktoren[vor, Byte(((not senkrecht) and (i=xmin)) or (senkrecht and (j=tmin))) +2*Byte(((not senkrecht) and (2*i=xmin+xmax+1)) or (senkrecht and (2*j=tmin+tmax+1)))]; end; setlength(umsortierung,pmax+1-pmin); for j:=0 to pmax-pmin do begin k:=0; for i:=0 to n-1 do // Bits umkehren k:=2*k + byte(odd(j shr i)); umsortierung[j]:=k; end; setlength(ims,pmax+1-pmin); setlength(wRe,haL); setlength(wIm,haL); for i:=0 to haL-1 do begin wRe[i]:=cos(pi*2*i/(pmax+1-pmin)); wIm[i]:=sin(pi*2*i/(pmax+1-pmin)); end; in0:=true; out0:=true; for i:=smin to smax do begin k:=1-2*byte(invers); // Eingang komplex konjugieren case vor of doResIms: begin // Re_0 Re_1=Re_n-1, ... Re_n/2-1=Re_n/2+1 Re_n/2 Im_1=-Im_n-1, ... Im_n/2-1=-Im_n/2+1 if pmin<>0 then begin // pmin != 0 wird bei dieser Eingangscodierung nicht geduldet (warum auch?) gibAus('Bei getrennten Imaginär- und Realteilen wird eine Verschiebung der Eingangsdaten in Transformationsrichtung nicht geduldet - wie kommt soetwas denn zustande?',3); exit; end; if senkrecht then begin for j:=1 to haL-1 do begin Ims[j]:=k*werte[i+(haL+j)*params.xsteps]; Ims[2*haL-j]:=-k*werte[i+(haL+j)*params.xsteps]; werte[i+(haL+j)*params.xsteps]:=werte[i+(haL-j)*params.xsteps]; end; end else begin for j:=1 to haL-1 do begin Ims[j]:=k*werte[haL+j+i*params.xsteps]; Ims[2*haL-j]:=-k*werte[haL+j+i*params.xsteps]; werte[haL+j+i*params.xsteps]:=werte[haL-j+i*params.xsteps]; end; end; ims[0]:=0; ims[haL]:=0; end; doResSmi: begin // Re_0 Re_1=Re_n-1, ... Re_n/2-1=Re_n/2+1 Re_n/2 -Im_n/2+1=Im_n/2-1, ... -Im_n-1=Im_1 if pmin<>0 then begin // pmin != 0 wird bei dieser Eingangscodierung nicht geduldet (warum auch?) gibAus('Bei getrennten Imaginär- und Realteilen wird eine Verschiebung der Eingangsdaten in Transformationsrichtung nicht geduldet - wie kommt soetwas denn zustande?',3); exit; end; if senkrecht then begin for j:=1 to haL-1 do begin Ims[j]:=k*werte[i+(2*haL-j)*params.xsteps]; Ims[2*haL-j]:=-Ims[j]; werte[i+(2*haL-j)*params.xsteps]:=werte[i+(j)*params.xsteps]; end; end else begin for j:=1 to haL-1 do begin Ims[j]:=k*werte[2*haL-j+i*params.xsteps]; Ims[2*haL-j]:=-Ims[j]; werte[2*haL-j+i*params.xsteps]:=werte[j+i*params.xsteps]; end; end; ims[0]:=0; ims[haL]:=0; end; doRes: // im = 0 for j:=0 to pmax-pmin do ims[j]:=0; doBetr,doBetrQdr: begin gibAus('Ich brauche mehr als Beträge oder Betragsquadrate um eine FFT auzurechnen!',3); exit; end; end{of case}; if senkrecht then begin for j:=0 to tmax-tmin do begin if umsortierung[j]>j then begin t1:=werte[i+(j+tmin)*params.xsteps]; werte[i+(j+tmin)*params.xsteps]:=werte[i+(umsortierung[j]+tmin)*params.xsteps]; werte[i+(umsortierung[j]+tmin)*params.xsteps]:=t1; t1:=ims[j]; ims[j]:=ims[umsortierung[j]]; ims[umsortierung[j]]:=t1; end; in0:=in0 and (werte[i+(j+tmin)*params.xsteps]=0) and (ims[j]=0); end; end else begin for j:=0 to xmax-xmin do begin if umsortierung[j]>j then begin t1:=werte[j+xmin+i*params.xsteps]; werte[j+xmin+i*params.xsteps]:=werte[umsortierung[j]+xmin+i*params.xsteps]; werte[umsortierung[j]+xmin+i*params.xsteps]:=t1; t1:=ims[j]; ims[j]:=ims[umsortierung[j]]; ims[umsortierung[j]]:=t1; end; in0:=in0 and (werte[j+xmin+i*params.xsteps]=0) and (ims[j]=0); end; end; dist:=1; wstep:=haL; if senkrecht then begin while dist werte[ i + {0..tmax-tmin}*xsteps ] doResIms: begin // Re_0 Re_1 Re_2 ... Re_n/2 Im_1 Im_2 ... Im_n/2-1 for j:=0 to haL do // Realteile zusammensuchen werte[i+j*params.xsteps]:=werte[i+(j+tmin)*params.xsteps]/sqrt(tmax+1-tmin); for j:=1 to haL - 1 do // Imaginärteile zusammensuchen werte[i+(j+haL)*params.xsteps]:=ims[j]/sqrt(tmax+1-tmin); end; doResSmi: begin // Re_0 Re_1 Re_2 ... Re_n/2 Im_n/2-1 Im_n/2-2 ... Im_1 for j:=0 to haL do // Realteile zusammensuchen werte[i+j*params.xsteps]:=werte[i+(j+tmin)*params.xsteps]/sqrt(tmax+1-tmin); for j:=1 to haL - 1 do // Imaginärteile zusammensuchen werte[i+(2*haL-j)*params.xsteps]:=ims[j]/sqrt(tmax+1-tmin); end; doBetr: // Abs(Re + i*Im) = Sqrt(Re^2 + Im^2) for j:=0 to tmax-tmin do werte[i+j*params.xsteps]:=sqrt((werte[i+(j+tmin)*params.xsteps]*werte[i+(j+tmin)*params.xsteps]+sqr(ims[j]))/(tmax+1-tmin)); doBetrQdr: // Re^2 + Im^2 for j:=0 to tmax-tmin do werte[i+j*params.xsteps]:=(werte[i+(j+tmin)*params.xsteps]*werte[i+(j+tmin)*params.xsteps]+sqr(ims[j]))/(tmax+1-tmin); doRes: // Re for j:=0 to tmax-tmin do werte[i+j*params.xsteps]:=werte[i+(j+tmin)*params.xsteps]/sqrt(tmax+1-tmin); else begin gibAus('Flieht, ihr Narren! (diese Meldung sollte eigentlich nicht auftreten können)',3); exit; end; end; for j:=0 to tmax-tmin do out0:=out0 and (werte[i+j*params.xsteps]=0); end else begin case nach of // werte[ {xtmin..xmax}+i*xsteps] & ims[{0..xmax-xmin}] -> werte[ {0..xmax-xmin}+i*xsteps ] doResIms: begin // Re_0 Re_1 Re_2 ... Re_n/2 Im_1 Im_2 ... Im_n/2-1 for j:=0 to haL do // Realteile zusammensuchen werte[j+i*params.xsteps]:=werte[j+xmin+i*params.xsteps]/sqrt(xmax+1-xmin); for j:=1 to haL - 1 do // Imaginärteile zusammensuchen werte[j+haL+i*params.xsteps]:=ims[j]/sqrt(xmax+1-xmin); end; doResSmi: begin // Re_0 Re_1 Re_2 ... Re_n/2 Im_n/2-1 Im_n/2-2 ... Im_1 for j:=0 to haL do // Realteile zusammensuchen werte[j+i*params.xsteps]:=werte[j+xmin+i*params.xsteps]/sqrt(xmax+1-xmin); for j:=1 to haL - 1 do // Imaginärteile zusammensuchen werte[2*haL-j+i*params.xsteps]:=ims[j]/sqrt(xmax+1-xmin); end; doBetr: // Abs(Re + i*Im) = Sqrt(Re^2 + Im^2) for j:=0 to xmax-xmin do werte[j+i*params.xsteps]:=sqrt((werte[j+xmin+i*params.xsteps]*werte[j+xmin+i*params.xsteps]+sqr(ims[j]))/(xmax+1-xmin)); doBetrQdr: // Re^2 + Im^2 for j:=0 to xmax-xmin do werte[j+i*params.xsteps]:=(werte[j+xmin+i*params.xsteps]*werte[j+xmin+i*params.xsteps]+sqr(ims[j]))/(xmax+1-xmin); doRes: // Re for j:=0 to xmax-xmin do werte[j+i*params.xsteps]:=werte[j+xmin+i*params.xsteps]/sqrt(xmax+1-xmin); else begin gibAus('Flieht, ihr Narren! (diese Meldung sollte eigentlich nicht auftreten können)',3); exit; end; end; for j:=0 to xmax-xmin do out0:=out0 and (werte[j+i*params.xsteps]=0); end; end; gibAus(inttostr(xmin)+'-'+inttostr(xmax)+'x'+inttostr(tmin)+'-'+inttostr(tmax),1); if senkrecht then begin tmax:=tmax-tmin; tmin:=0; end else begin xmax:=xmax-xmin; xmin:=0; end; nachher:=0; if nach=doBetrQdr then begin for i:=xmin to xmax do for j:=tmin to tmax do nachher:=nachher + werte[i+j*params.xsteps]; end else begin for i:=xmin to xmax do for j:=tmin to tmax do nachher:=nachher + werte[i+j*params.xsteps]*werte[i+j*params.xsteps] *faktoren[nach, Byte(((not senkrecht) and (i=xmin)) or (senkrecht and (j=tmin))) +2*Byte(((not senkrecht) and (2*i=xmin+xmax+1)) or (senkrecht and (2*j=tmin+tmax+1)))]; end; if (nachher=0) and (vorher=0) then pvFehler:=0 else pvFehler:=abs(nachher-vorher)/(nachher+vorher); gibAus(inttostr(byte(senkrecht))+' '+inttostr(byte(vor))+' '+inttostr(byte(nach))+' '+inttostr(byte(invers))+' (Parseval-Fehler = '+floattostr(pvFehler)+') ... nämlich von '+floattostr(vorher)+' zu '+floattostr(nachher),1); if in0 then gibAus('Nur Nullen im Input der FFT!',1); if out0 then gibAus('Nur Nullen im Output der FFT!',1); result:=not (in0 or out0); end; procedure tLLWerte.schreibeWert(var f: textfile; x,y: longint); begin schreibeWert(f,x,y,werte[x+y*params.xsteps]); end; procedure tLLWerte.schreibeWert(var f: textfile; x,y,wert: extended); var xa,ta: extended; begin if params.xstop=params.xstart then xa:=params.xstop else xa:=params.xstart+x/(params.xsteps-1)*(params.xstop-params.xstart); if params.tstop=params.tstart then ta:=params.tstop else ta:=params.tstart+y/(params.tsiz-1)*(params.tstop-params.tstart); writeln(f,floattostr(xa)+' '+floattostr(ta)+' '+floattostr(wert)); end; procedure tLLWerte.schreibeWertIntegriert(var f: textfile; i: longint; hor: boolean); 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]; schreibeWert(f,(params.xsteps-1)/2,i,tmp); end else begin for j:=0 to params.tsiz-1 do tmp:=tmp+werte[i+j*params.xsteps]; schreibeWert(f,i,(params.tsiz-1)/2,tmp); end; 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: longint; begin case NB of doResIms,doResSmi: ; doBetr: begin werte[0]:=abs(extended(werte[0])); // rein reell for i:=1 to params.xsteps div 2 -1 do begin // unterste Zeile ist reell in t werte[i]:=sqrt(sqr(extended(werte[i]))+sqr(extended(werte[params.xsteps-i]))); werte[params.xsteps-i]:=werte[i]; end; for i:=1 to params.xsteps div 2 -1 do begin // mittlere Zeile ist reell in t werte[i+(params.tsiz div 2)*params.xsteps]:=sqrt(sqr(extended(werte[i+(params.tsiz div 2)*params.xsteps]))+sqr(extended(werte[params.xsteps-i+(params.tsiz div 2)*params.xsteps]))); werte[params.xsteps-i+(params.tsiz div 2)*params.xsteps]:=werte[i+(params.tsiz div 2)*params.xsteps]; end; werte[params.xsteps div 2]:=abs(extended(werte[params.xsteps div 2])); // rein reell werte[params.tsiz div 2]:=abs(extended(werte[params.tsiz div 2])); // rein reell for i:=1 to params.tsiz div 2 -1 do begin // linkeste Spalte ist reell in x werte[i*params.xsteps]:=sqrt(sqr(extended(werte[i*params.xsteps]))+sqr(extended(werte[(params.tsiz-i)*params.xsteps]))); werte[(params.tsiz-i)*params.xsteps]:=werte[i*params.xsteps]; end; for i:=1 to params.tsiz div 2 -1 do begin // mittlere Spalte ist reell in x werte[params.xsteps div 2 + i*params.xsteps]:=sqrt(sqr(extended(werte[params.xsteps div 2 + i*params.xsteps]))+sqr(extended(werte[params.xsteps div 2 + (params.tsiz-i)*params.xsteps]))); werte[params.xsteps div 2 + (params.tsiz-i)*params.xsteps]:=werte[params.xsteps div 2 + i*params.xsteps]; end; end; doBetrQdr: begin werte[0]:=sqr(extended(werte[0])); // rein reell for i:=1 to params.xsteps div 2 -1 do begin // unterste Zeile ist reell in t werte[i]:=sqr(extended(werte[i]))+sqr(extended(werte[params.xsteps-i])); werte[params.xsteps-i]:=werte[i]; end; for i:=1 to params.xsteps div 2 -1 do begin // mittlere Zeile ist reell in t werte[i+(params.tsiz div 2)*params.xsteps]:=sqr(extended(werte[i+(params.tsiz div 2)*params.xsteps]))+sqr(extended(werte[params.xsteps-i+(params.tsiz div 2)*params.xsteps])); werte[params.xsteps-i+(params.tsiz div 2)*params.xsteps]:=werte[i+(params.tsiz div 2)*params.xsteps]; end; werte[params.xsteps div 2]:=sqr(extended(werte[params.xsteps div 2])); // rein reell werte[params.tsiz div 2]:=sqr(extended(werte[params.tsiz div 2])); // rein reell for i:=1 to params.tsiz div 2 -1 do begin // linkeste Spalte ist reell in x werte[i*params.xsteps]:=sqr(extended(werte[i*params.xsteps]))+sqr(extended(werte[(params.tsiz-i)*params.xsteps])); werte[(params.tsiz-i)*params.xsteps]:=werte[i*params.xsteps]; end; for i:=1 to params.tsiz div 2 -1 do begin // mittlere Spalte ist reell in x werte[params.xsteps div 2 + i*params.xsteps]:=sqr(extended(werte[params.xsteps div 2 + i*params.xsteps]))+sqr(extended(werte[params.xsteps div 2 + (params.tsiz-i)*params.xsteps])); werte[params.xsteps div 2 + (params.tsiz-i)*params.xsteps]:=werte[params.xsteps div 2 + i*params.xsteps]; end; end; end{of case}; end; procedure tLLWerte.fft2dNachbearbeitungB(xmin,xmax: longint; nb: tFFTDatenordnung); var i,j: longint; begin // bearbeitet nur den Hauptteil (außer erster und mittlerer Zeile/Spalte) nach! case nb of doBetr: begin for i:=xmin to xmax do for j:=1 to params.tsiz div 2 -1 do begin werte[i+j*params.xsteps]:= sqrt(sqr(extended(werte[i+j*params.xsteps]-werte[params.xsteps-i+(params.tsiz-j)*params.xsteps])) // Re^2 +sqr(extended(werte[params.xsteps-i+j*params.xsteps]+werte[i+(params.tsiz-j)*params.xsteps]))); // Im^2 werte[params.xsteps-i+j*params.xsteps]:=werte[i+j*params.xsteps]; werte[i+(params.tsiz-j)*params.xsteps]:=werte[i+j*params.xsteps]; werte[params.xsteps-i+(params.tsiz-j)*params.xsteps]:=werte[i+j*params.xsteps]; end; end; doBetrQdr: begin for i:=xmin to xmax do for j:=1 to params.tsiz div 2 -1 do begin werte[i+j*params.xsteps]:= sqr(extended(werte[i+j*params.xsteps]-werte[params.xsteps-i+(params.tsiz-j)*params.xsteps])) // Re^2 +sqr(extended(werte[params.xsteps-i+j*params.xsteps]+werte[i+(params.tsiz-j)*params.xsteps])); // Im^2 werte[params.xsteps-i+j*params.xsteps]:=werte[i+j*params.xsteps]; werte[i+(params.tsiz-j)*params.xsteps]:=werte[i+j*params.xsteps]; werte[params.xsteps-i+(params.tsiz-j)*params.xsteps]:=werte[i+j*params.xsteps]; end; end; end{of case} end; 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; pPWerte: pTExtendedArray; pPAnzahlen: pTLongintArray): boolean; var i,j,k,l, xv,xb,tv,tb: longint; b: boolean; begin result:=false; for i:=0 to length(pPWerte^)-1 do begin pPWerte^[i]:=0; pPAnzahlen^[i]:=0; end; b:=false; for j:=0 to whoehe-1 do for i:=0 to wbreite-1 do begin xv:=min(params.xsteps-1,max(0,ceil((i+max(0,xpmi)-1/2)/xz+xmi))); xb:=min(params.xsteps-1,max(0,ceil((i+max(0,xpmi)+1/2)/xz+xmi-1))); tv:=min(params.tsiz-1,max(0,ceil((j-1/2)/yz+tmi))); tb:=min(params.tsiz-1,max(0,ceil((j+1/2)/yz+tmi-1))); if xv>xb then begin xv:=xv-byte((i>0) or (xpmi>0)); xb:=xb+byte((xpmi<0) or ((i=0) and (xpmi=0))); end; if tv>tb then begin tv:=tv-byte(j>0); tb:=tb+byte(j=0); end; if (xv>xb) or (tv>tb) then begin gibAus('Keine Inputwerte für Position '+inttostr(i)+':'+inttostr(j)+'!',1); exit; end; pPanzahlen^[j*wbreite+i]:=(xb-xv+1)*(tb-tv+1); for k:=xv to xb do for l:=tv to tb do begin b:=b or (werte[k+l*params.xsteps]<>0); pPWerte^[j*wbreite+i]:=pPWerte^[j*wbreite+i]+min(params.maxW,max(params.minW,werte[k+l*params.xsteps])); end; end; if not b then gibAus('Habe nur Nullen im Input!',1); result:=true; end; function tLLWerte.findeSchwellwerte(xmi,xma,tmi,tma: longint; Schw: extended): tExtPointArray; var i,j,k,l,m,vz: longint; dx,dy,x0,y0: extended; begin setlength(result,0); gibAus('Schwellwerte finden ('+inttostr(xmi)+'-'+inttostr(xma)+') '+floattostr(Schw)+' ...',1); i:=0; dx:=(params.xstop-params.xstart)/(params.xsteps-1); x0:=params.xstart; dy:=(params.tstop-params.tstart)/(params.tsiz-1); y0:=params.tstart; for j:=xmi to xma do begin for k:=tmi to tma do begin vz:=0; for l:=0 to 1 do for m:=0 to 1 do vz:= vz or (byte(werte[j-l+(k-m)*params.xsteps]>=Schw) shl 1) or byte(werte[j-l+(k-m)*params.xsteps]<=Schw); if vz=3 then begin if i>=length(result) then setlength(result,length(result)+Speicherhappen); result[i]['x']:=j*dx+x0; result[i]['y']:=k*dy+y0; inc(i); end; end; if (xma-j) and $ff = 0 then gibAus('x = '+inttostr(j)+' ('+inttostr(xmi)+'-'+inttostr(xma)+')',1); end; gibAus('... fertig!',1); setlength(result,i); end; procedure tLLWerte.integriereSingle(qu: pTLLWerteSingle; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung); var i,j: longint; int,faktor: extended; begin case richtung of irHorizontal: begin faktor:=(qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1); for i:=tmi to tma do begin int:=0; for j:=0 to xmi-1 do int:=int+qu^.werte[j + i*qu^.params.xsteps]; for j:=xmi to xma do begin int:=int+qu^.werte[j + i*qu^.params.xsteps]; werte[j-xof + (i-tof)*params.xsteps]:=int*faktor; end; end; end; irEinfall: begin faktor:= sqrt( sqr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)) + sqr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1))); gibAus('dx = '+floattostr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)),1); gibAus('dt = '+floattostr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1)),1); for i:=tmi to tma do begin // von links eintretendes (inkl. Ecke links unten) int:=0; for j:=1 to min(xmi,i) do int:=int+qu^.werte[xmi-j + (i-j)*qu^.params.xsteps]; for j:=0 to min(tma-i,xma-xmi) do begin int:=int+qu^.werte[xmi+j + (i+j)*qu^.params.xsteps]; werte[j+xmi-xof + (i+j-tof)*params.xsteps]:=int*faktor; end; end; for i:=xmi+1 to xma do begin // von unten eintretendes (exkl. Ecke links unten) int:=0; for j:=1 to min(tmi,i) do int:=int+qu^.werte[i-j + (tmi-j)*qu^.params.xsteps]; for j:=0 to min(tma-tmi,xma-i) do begin int:=int+qu^.werte[i+j + (tmi+j)*qu^.params.xsteps]; werte[i+j-xof + (tmi+j-tof)*params.xsteps]:=int*faktor; end; end; end; irAusfall: begin faktor:= sqrt( sqr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)) + sqr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1))); gibAus('dx = '+floattostr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)),1); gibAus('dt = '+floattostr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1)),1); for i:=tmi to tma do begin // nach links austretendes (inkl. Ecke links oben) int:=0; for j:=1 to min(xmi,qu^.params.tsiz-1-i) do int:=int+qu^.werte[xmi-j + (i+j)*qu^.params.xsteps]; for j:=0 to min(i-tmi,xma-xmi) do begin int:=int+qu^.werte[xmi+j + (i-j)*qu^.params.xsteps]; werte[j+xmi-xof + (i-j-tof)*params.xsteps]:=int*faktor; end; end; for i:=xmi+1 to xma do begin // nach oben austretendes (exkl. Ecke links oben) int:=0; for j:=1 to min(qu^.params.tsiz-1-tma,i) do int:=int+qu^.werte[i-j + (tma+j)*qu^.params.xsteps]; for j:=0 to min(tma-tmi,xma-i) do begin int:=int+qu^.werte[i+j + (tma-j)*qu^.params.xsteps]; werte[i+j-xof + (tma-j-tof)*params.xsteps]:=int*faktor; end; end; end; end{of case}; end; procedure tLLWerte.integriereDouble(qu: pTLLWerteDouble; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung); var i,j: longint; int,faktor: extended; begin case richtung of irHorizontal: begin faktor:=(qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1); for i:=tmi to tma do begin int:=0; for j:=0 to xmi-1 do int:=int+qu^.werte[j + i*qu^.params.xsteps]; for j:=xmi to xma do begin int:=int+qu^.werte[j + i*qu^.params.xsteps]; werte[j-xof + (i-tof)*params.xsteps]:=int*faktor; end; end; end; irEinfall: begin faktor:= sqrt( sqr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)) + sqr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1))); gibAus('dx = '+floattostr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)),1); gibAus('dt = '+floattostr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1)),1); for i:=tmi to tma do begin // von links eintretendes (inkl. Ecke links unten) int:=0; for j:=1 to min(xmi,i) do int:=int+qu^.werte[xmi-j + (i-j)*qu^.params.xsteps]; for j:=0 to min(tma-i,xma-xmi) do begin int:=int+qu^.werte[xmi+j + (i+j)*qu^.params.xsteps]; werte[j+xmi-xof + (i+j-tof)*params.xsteps]:=int*faktor; end; end; for i:=xmi+1 to xma do begin // von unten eintretendes (exkl. Ecke links unten) int:=0; for j:=1 to min(tmi,i) do int:=int+qu^.werte[i-j + (tmi-j)*qu^.params.xsteps]; for j:=0 to min(tma-tmi,xma-i) do begin int:=int+qu^.werte[i+j + (tmi+j)*qu^.params.xsteps]; werte[i+j-xof + (tmi+j-tof)*params.xsteps]:=int*faktor; end; end; end; irAusfall: begin faktor:= sqrt( sqr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)) + sqr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1))); gibAus('dx = '+floattostr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)),1); gibAus('dt = '+floattostr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1)),1); for i:=tmi to tma do begin // nach links austretendes (inkl. Ecke links oben) int:=0; for j:=1 to min(xmi,qu^.params.tsiz-1-i) do int:=int+qu^.werte[xmi-j + (i+j)*qu^.params.xsteps]; for j:=0 to min(i-tmi,xma-xmi) do begin int:=int+qu^.werte[xmi+j + (i-j)*qu^.params.xsteps]; werte[j+xmi-xof + (i-j-tof)*params.xsteps]:=int*faktor; end; end; for i:=xmi+1 to xma do begin // nach oben austretendes (exkl. Ecke links oben) int:=0; for j:=1 to min(qu^.params.tsiz-1-tma,i) do int:=int+qu^.werte[i-j + (tma+j)*qu^.params.xsteps]; for j:=0 to min(tma-tmi,xma-i) do begin int:=int+qu^.werte[i+j + (tma-j)*qu^.params.xsteps]; werte[i+j-xof + (tma-j-tof)*params.xsteps]:=int*faktor; end; end; end; end{of case}; end; procedure tLLWerte.integriereExtended(qu: pTLLWerteDouble; xmi,xma,tmi,tma,xof,tof: longint; richtung: tIntegrationsRichtung); var i,j: longint; int,faktor: extended; begin case richtung of irHorizontal: begin faktor:=(qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1); for i:=tmi to tma do begin int:=0; for j:=0 to xmi-1 do int:=int+qu^.werte[j + i*qu^.params.xsteps]; for j:=xmi to xma do begin int:=int+qu^.werte[j + i*qu^.params.xsteps]; werte[j-xof + (i-tof)*params.xsteps]:=int*faktor; end; end; end; irEinfall: begin faktor:= sqrt( sqr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)) + sqr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1))); gibAus('dx = '+floattostr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)),1); gibAus('dt = '+floattostr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1)),1); for i:=tmi to tma do begin // von links eintretendes (inkl. Ecke links unten) int:=0; for j:=1 to min(xmi,i) do int:=int+qu^.werte[xmi-j + (i-j)*qu^.params.xsteps]; for j:=0 to min(tma-i,xma-xmi) do begin int:=int+qu^.werte[xmi+j + (i+j)*qu^.params.xsteps]; werte[j+xmi-xof + (i+j-tof)*params.xsteps]:=int*faktor; end; end; for i:=xmi+1 to xma do begin // von unten eintretendes (exkl. Ecke links unten) int:=0; for j:=1 to min(tmi,i) do int:=int+qu^.werte[i-j + (tmi-j)*qu^.params.xsteps]; for j:=0 to min(tma-tmi,xma-i) do begin int:=int+qu^.werte[i+j + (tmi+j)*qu^.params.xsteps]; werte[i+j-xof + (tmi+j-tof)*params.xsteps]:=int*faktor; end; end; end; irAusfall: begin faktor:= sqrt( sqr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)) + sqr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1))); gibAus('dx = '+floattostr((qu^.params.xstop-qu^.params.xstart)/(qu^.params.xsteps-1)),1); gibAus('dt = '+floattostr((qu^.params.tstop-qu^.params.tstart)/(qu^.params.tsiz-1)),1); for i:=tmi to tma do begin // nach links austretendes (inkl. Ecke links oben) int:=0; for j:=1 to min(xmi,qu^.params.tsiz-1-i) do int:=int+qu^.werte[xmi-j + (i+j)*qu^.params.xsteps]; for j:=0 to min(i-tmi,xma-xmi) do begin int:=int+qu^.werte[xmi+j + (i-j)*qu^.params.xsteps]; werte[j+xmi-xof + (i-j-tof)*params.xsteps]:=int*faktor; end; end; for i:=xmi+1 to xma do begin // nach oben austretendes (exkl. Ecke links oben) int:=0; for j:=1 to min(qu^.params.tsiz-1-tma,i) do int:=int+qu^.werte[i-j + (tma+j)*qu^.params.xsteps]; for j:=0 to min(tma-tmi,xma-i) do begin int:=int+qu^.werte[i+j + (tma-j)*qu^.params.xsteps]; werte[i+j-xof + (tma-j-tof)*params.xsteps]:=int*faktor; end; end; end; end{of case}; end; // tWavelet ******************************************************************** function tWavelet.setzeTyp(s: string): boolean; begin result:=true; if s='Sin2' then begin typ:=wtSin2; exit; end; if s='Frequenzfenster' then begin typ:=wtFrequenzfenster; exit; end; gibAus('Kenne Wavelettyp '''+s+''' nicht!',3); result:=false; end; function tWavelet.berechneWerte: boolean; var i: longint; tmp: extended; fenster: tFenster; nur0: boolean; begin result:=false; werte.params.xsteps:=2; werte.params.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; fenster.aktiv:=false; // fenster.Breite:=0; // fenster.berechneWerte(werte.tsiz); gibAus(inttostr(werte.params.xsteps)+' '+inttostr(werte.params.tsiz)+' '+inttostr(length(werte.werte)),1); werte.fft(0,1,0,werte.params.tsiz-1,true,false,doRes,doResIms,fenster,pvFehler); end; wtFrequenzfenster: begin hlen:=werte.params.tsiz div 2; for i:=0 to hlen do begin werte.werte[2*i]:=Byte(tfwhm*Abs(i/werte.params.tsiz-freq)<=1); // Re=1 <=> |f-f_Mitte| < Laenge/T_fwhm werte.werte[2*i+1]:=0; // Dummy-Wavelet end; for i:=hlen+1 to 2*hlen-1 do begin werte.werte[2*i]:=0; // Ims = 0 werte.werte[2*i+1]:=0; // Dummy-Wavelet end; 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; var ps: tExtrainfos; begin inherited create; ps:=tExtrainfos.create; werte:=tLLWerteDouble.create(ps); pvFehler:=0; end; destructor tWavelet.destroy; begin werte.params.Free; werte.free; inherited destroy; end; end.