unit werteunit; {$mode objfpc}{$H+} interface uses classes, sysutils, typenunit, math, process, lowlevelunit, matheunit, fftunit; 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 procedure sortiereMaxima(var maxima: tIntPointArray); 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: extended; out pvFehler: extended): boolean; overload; inline; function fft(sMin,sMax: longint; senkrecht,invers: boolean; algo: tFFTAlgorithmus; fen: tFenster; hg: extended; 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 fft2dNachbearbeitungVerdoppeln(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,komplex: 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); procedure gauszFit(amplituden,breiten,positionen,ueberlappe,hintergruende: pTLLWerteExtended; von,bis: longint; senkrecht: boolean; fensterBreite,maxBreite,maxVerschiebung: extended; positionsMitten: tExtendedArray); function ermittleHintergrund: extended; procedure kantenFilter(betraege: tLLWerte; xFak,yFak: extended; filterTyp: tKantenFilterTyp); overload; procedure kantenFilter(betraege: tLLWerte; xFak,yFak: extended; filterTyp: tKantenFilterTyp; einseitig: boolean; out maxPos: tIntPoint); overload; procedure fenstereWerte(xMi,xMa,tMi,tMa: int64; xFen,tFen: tFenster; hg: extended); procedure verschiebe(richtung: tIntPoint; xV,xB,tV,tB: longint; komplex: boolean); 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; const anzSergeyFelder = 12; implementation uses systemunit; // tLLWerte ******************************************************************** procedure tLLWerte.sortiereMaxima(var maxima: tIntPointArray); var mins,maxs: tExtendedArray; pivot,wert,wertLi,wertRe: extended; vons,biss: tLongintArray; tmp: tIntPoint; i,li,re,cnt: int64; begin setLength(vons,1); vons[0]:=0; setLength(biss,1); biss[0]:=length(maxima)-1; setLength(mins,1); setLength(maxs,1); mins[0]:=0; maxs[0]:=0; for i:=vons[0] to biss[0] do begin wert:=werte[maxima[i]['x']+maxima[i]['y']*params.xSteps]; if (i=vons[0]) or (wert>maxs[0]) then maxs[0]:=wert; if (i=vons[0]) or (wert0 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; pivot:=sqrt(maxs[cnt-1]*mins[cnt-1]);//(maxs[cnt-1]+mins[cnt-1])/2; 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[maxima[li]['x']+maxima[li]['y']*params.xSteps]; if wertLi>=pivot then begin if wertLimaxs[cnt] then maxs[cnt]:=wertRe; dec(re); continue; end; if wertLi>maxs[cnt] then maxs[cnt]:=wertLi; if wertRe=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; procedure tLLWerte.kopiereLOVerzerrt(original: pTLLWerteSingle; xMin,xMax,tMin,tMax: longint; verhHo,verhVe: extended); overload; var i,j,hV,hB,vV,vB,h,v: int64; begin writeln(xMin,' .. ',xMax,' x ',tMin,' .. ',tMax); for i:=tMin to tMax do begin if verhVe>0 then begin vV:=max( 0, round( (params.tSiz-1-i+0.5)/ (1 + (i+0.5)/verhVe/(params.tSiz-1)) ) ); vB:=min( params.tSiz-1, round( (params.tSiz-1-i+0.5)/ (1 + (i-0.5)/verhVe/(params.tSiz-1)) ) ); end else begin vV:=i; vB:=i; end; for j:=xMin to xMax do begin if verhHo>0 then begin hV:=max( 0, round( (params.xSteps-1-(j+0.5))/ (1 + (j+0.5)/verhHo/(params.xSteps-1)) ) ); hB:=min( params.xSteps-1, round( (params.xSteps-1-(j-0.5))/ (1 + (j-0.5)/verhHo/(params.xSteps-1)) ) ); end else begin hV:=j; hB:=j; end; werte[j+i*params.xSteps]:=0; for h:=hV to hB do for v:=vV to vB do werte[j+i*params.xSteps]:= werte[j+i*params.xSteps]+ original^.werte[h+v*params.xSteps]; if (hB>=hV) and (vB>=vV) then werte[j+i*params.xSteps]:= werte[j+i*params.xSteps] / (hB+1-hV) / (vB+1-vV); end; end; end; procedure tLLWerte.kopiereLOVerzerrt(original: pTLLWerteDouble; xMin,xMax,tMin,tMax: longint; verhHo,verhVe: extended); overload; var i,j,hV,hB,vV,vB,h,v: int64; begin writeln(xMin,' .. ',xMax,' x ',tMin,' .. ',tMax); for i:=tMin to tMax do begin if verhVe>0 then begin vV:=max( 0, round( (params.tSiz-1-i+0.5)/ (1 + (i+0.5)/verhVe/(params.tSiz-1)) ) ); vB:=min( params.tSiz-1, round( (params.tSiz-1-i+0.5)/ (1 + (i-0.5)/verhVe/(params.tSiz-1)) ) ); end else begin vV:=i; vB:=i; end; for j:=xMin to xMax do begin if verhHo>0 then begin hV:=max( 0, round( (params.xSteps-1-(j+0.5))/ (1 + (j+0.5)/verhHo/(params.xSteps-1)) ) ); hB:=min( params.xSteps-1, round( (params.xSteps-1-(j-0.5))/ (1 + (j-0.5)/verhHo/(params.xSteps-1)) ) ); end else begin hV:=j; hB:=j; end; werte[j+i*params.xSteps]:=0; for h:=hV to hB do for v:=vV to vB do werte[j+i*params.xSteps]:= werte[j+i*params.xSteps]+ original^.werte[h+v*params.xSteps]; if (hB>=hV) and (vB>=vV) then werte[j+i*params.xSteps]:= werte[j+i*params.xSteps] / (hB+1-hV) / (vB+1-vV); end; end; end; procedure tLLWerte.kopiereLOVerzerrt(original: pTLLWerteExtended; xMin,xMax,tMin,tMax: longint; verhHo,verhVe: extended); overload; var i,j,hV,hB,vV,vB,h,v: int64; begin writeln(xMin,' .. ',xMax,' x ',tMin,' .. ',tMax); for i:=tMin to tMax do begin if verhVe>0 then begin vV:=max( 0, round( (params.tSiz-1-i+0.5)/ (1 + (i+0.5)/verhVe/(params.tSiz-1)) ) ); vB:=min( params.tSiz-1, round( (params.tSiz-1-i+0.5)/ (1 + (i-0.5)/verhVe/(params.tSiz-1)) ) ); end else begin vV:=i; vB:=i; end; for j:=xMin to xMax do begin if verhHo>0 then begin hV:=max( 0, round( (params.xSteps-1-(j+0.5))/ (1 + (j+0.5)/verhHo/(params.xSteps-1)) ) ); hB:=min( params.xSteps-1, round( (params.xSteps-1-(j-0.5))/ (1 + (j-0.5)/verhHo/(params.xSteps-1)) ) ); end else begin hV:=j; hB:=j; end; werte[j+i*params.xSteps]:=0; for h:=hV to hB do for v:=vV to vB do werte[j+i*params.xSteps]:= werte[j+i*params.xSteps]+ original^.werte[h+v*params.xSteps]; if (hB>=hV) and (vB>=vV) then werte[j+i*params.xSteps]:= werte[j+i*params.xSteps] / (hB+1-hV) / (vB+1-vV); end; 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; 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'' 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; 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; algo: tFFTAlgorithmus; fen: tFenster; hg: extended; 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: extended; 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:=byte(not invers)*hg; 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 (hg<>0) then for i:=sMin to sMax do // Hintergrund addieren for j:=0 to pMax do werte[i*sStep+j*pStep]:= werte[i*sStep+j*pStep]+hg; 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; 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; else raise exception.create('Nachbearbeitungsmethode '+fftDoToStr(nB)+' nicht implementiert!'); 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 doResIms,doResSmi: ; 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; else raise exception.create('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 raise exception.create('Nachbearbeitungsmethode '+fftDoToStr(nB)+' verdoppelt die Anzahl der Daten 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,komplex: longint; xZ,yZ: extended; pPWerte: pTExtendedArray; pPAnzahlen: pTLongintArray): boolean; var i,j,k,l, xV,xB,tV,tB: longint; b,imPart: boolean; 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; for imPart:=false to komplex>0 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(params.tSiz div (1+byte(komplex>0))-1,max(0,ceil((j-1/2)/yZ+tMi)))+byte(imPart)*komplex; tB:=min(params.tSiz div (1+byte(komplex>0))-1,max(0,ceil((j+1/2)/yZ+tMi-1)))+byte(imPart)*komplex; 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); pPWerte^[(j+wHoehe*byte(imPart))*wBreite+i]:= pPWerte^[(j+wHoehe*byte(imPart))*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; 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=100000) or (verbesserung<-10) or (fehlers[aParams]*10 params.tSiz) or (betraege.params.xSteps <> params.xSteps) or (not (istVollKomplex in [1,2])) then raise exception.create( '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 raise exception.create('Kann nur voll komplexe Werte einseitig per Kantenfilter filtern!'); if filterTyp<>kfHochpass then raise exception.create('Kann nur einen Hochpass als einseitigen Kantenfilter verwenden!'); end; setLength(maxima,0); mCnt:=0; for j:=0 to betraege.params.tSiz div 2 + 1 do begin jM:=j - 1 + betraege.params.tSiz*byte(j=0); 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+(j+1)*betraege.params.xSteps]) 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); writeln(length(maxima),' (von ',betraege.params.xSteps*betraege.params.tSiz,')'); betraege.sortiereMaxima(maxima); maxPos:=maxima[1]; maxWert:=0; for i:=1 to length(maxima)-1 do begin minWert:=0; for j:=0 to i-1 do begin wert:=sqr((maxima[j]['x']-maxima[i]['x'])*xFak)+sqr((maxima[j]['y']-maxima[i]['y'])*yFak); if (j=0) or (wert1) then break; end; end; if (i=1) 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; for j:=0 to jM do for i:=0 to iM do begin wert:=(sqr(i)*xFak+sqr(j)*yFak)/(sqr(maxPos['x'])*xFak+sqr(maxPos['y'])*yFak); 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 (i0) and (j0) and (i0) and (j 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; procedure tLLWerte.fenstereWerte(xMi,xMa,tMi,tMa: int64; xFen,tFen: tFenster; hg: extended); var i,j: int64; begin for j:=tMi to tMa do for i:=xMi to xMa do werte[i+j*params.xSteps]:= (werte[i+j*params.xSteps]-hg)*xFen.werte[i]*tFen.werte[j]; end; procedure tLLWerte.verschiebe(richtung: tIntPoint; xV,xB,tV,tB: longint; komplex: boolean); var xS,tS,x,t,xN,tN,xM,tM: longint; imPart: boolean; wert: wGen; begin xM:=params.xSteps; tM:=params.tSiz div (1+byte(komplex)); 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 raise exception.create('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 komplex do for xS:=xV to xB do for tS:=tV to tB do begin wert:=werte[xS+tS*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*xM]:=werte[xN+tN*xM]; until (xN=xS) and (tN=tS); werte[x+t*xM]:=wert; end; 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,0,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; 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.