summaryrefslogtreecommitdiff
path: root/romunit.pas
diff options
context:
space:
mode:
Diffstat (limited to 'romunit.pas')
-rw-r--r--romunit.pas1086
1 files changed, 1086 insertions, 0 deletions
diff --git a/romunit.pas b/romunit.pas
new file mode 100644
index 0000000..aaab446
--- /dev/null
+++ b/romunit.pas
@@ -0,0 +1,1086 @@
+unit romunit;
+
+{$mode objfpc}{$H+}
+
+interface
+
+uses
+ Classes, SysUtils, Math, mathunit;
+
+procedure Fehler(s: string);
+procedure readRawInputs(nam: string; out d1,d2: tExtPointArray; var absShift: extended);
+procedure readTextInput(nam: string; out dat: tExtPointArray);
+procedure writeOutput(nam: string; const dat: tExtPointArray);
+procedure berechneTrajektorie(const inPuls,outPuls: tExtPointArray; out surTraj: tExtPointArray; absShift: extended);
+procedure smoothen(var dat: tExtPointArray; width: longint);
+procedure sort(var dat: tExtPointArray); overload;
+procedure sort(var dat: tExtPointArray; start, stop, threads: longint); overload;
+procedure cut(var dat: tExtPointArray; tmax: extended);
+function findeExtrema(const dat: tExtPointArray; maxima: boolean): tLongintArray;
+procedure filtereExtrema(const dat: tExtPointArray; var extrema: tLongintArray; sollAbst, toleranz: extended);
+procedure monotonieHerstellen(const dat: tExtPointArray; var minima,maxima: tLongintArray);
+procedure gesamtverschiebung(var inPuls,outPuls: tExtPointArray; var absShift: extended);
+procedure integrate(var dat: tExtPointArray);
+procedure removeLinearOffset(var dat: tExtPointArray);
+procedure flip(var dat: tExtPointArray);
+procedure vereineExtrema(const dat1,dat2: tExtPointArray; var el1,el2: tLongintArray; toleranz: extended);
+procedure uniq(var dat: tExtPointArray; streng: boolean);
+procedure fft(var dat: tExtPointArray);
+procedure interpoliere(var dat: tExtPointArray);
+procedure normiere(var dat: tExtPointArray);
+
+type
+ tSortThread = class(tThread)
+ _pdat: pTExtPointArray;
+ _start,_stop: longint;
+ fertig: boolean;
+ constructor create(pd: pTExtPointArray; sta, sto: longint);
+ destructor destroy; override;
+ procedure execute; override;
+ end;
+ tPhaseShiftThread = class(tThread)
+ _dats: array[boolean] of tExtPointArray;
+ _ergebnis: tExtPointArray;
+ _starts,_stops: array[boolean] of longint;
+ _absShift: extended;
+ fertig: boolean;
+ constructor create(p1d,p2d: pTExtPointArray; sta1, sto1, sta2, sto2: longint; absShift: extended);
+ destructor destroy; override;
+ procedure execute; override;
+ end;
+
+implementation
+
+procedure Fehler(s: string);
+begin
+ writeln(s);
+ halt(1);
+end;
+
+procedure readRawInputs(nam: string; out d1,d2: tExtPointArray; var absShift: extended);
+var f: file;
+ tf: textfile;
+ s: string;
+ i,j,k,start,stop,traces,steps: longint;
+ i32: int32;
+ sr: tSearchRec;
+ fl: single;
+ buff: array of single;
+ factor,dx,cells_left,maxAmp: extended;
+begin
+ if not fileexists(nam) then Fehler('Die Datei '''+nam+''' existiert nicht!');
+ assignfile(tf,nam);
+ i:=0;
+ reset(tf);
+ factor:=1;
+ j:=0;
+ dx:=0;
+ cells_left:=0;
+ maxAmp:=0;
+ while not eof(tf) do begin
+ readln(tf,s);
+ if pos('output path :',s)=1 then begin
+ delete(s,1,pos(':',s));
+ while (length(s)>0) and (s[1]=' ') do
+ delete(s,1,1);
+ while (length(s)>0) and (s[length(s)]=' ') do
+ delete(s,length(s),1);
+ nam:=extractfilepath(nam)+s+'/trace-1-';
+ i:=i or 1;
+ end;
+ if s='pulse # 1' then
+ i:=i or 2;
+ if s='pulse # 2' then
+ i:=i and not $7E;
+ if s='pulse component # 0' then
+ i:=(i or $4) and not $78;
+ if s='pulse component # 1' then
+ i:=(i or $8) and not $74;
+ if s='pulse component # 2' then
+ i:=(i or $10) and not $7C;
+ if s='pulse component # 3' then
+ i:=(i or $20) and not $5C;
+ if odd(i shr 1) and (pos('.a0 :',s)=1) then begin
+ delete(s,1,pos(':',s));
+ while (length(s)>0) and (s[1]=' ') do
+ delete(s,1,1);
+ if strtofloat(s)>maxAmp then begin
+ i:=i or $40;
+ maxAmp:=strtofloat(s);
+ end;
+ end;
+ if odd(i shr 6) and (pos('.frequency :',s)=1) then begin
+ delete(s,1,pos(':',s));
+ while (length(s)>0) and (s[1]=' ') do
+ delete(s,1,1); writeln('bla'+s+'blu');
+ factor:=strtofloat(s);
+ end;
+ case j of
+ 0:
+ if s='domain: plasma' then
+ j:=1;
+ 1: begin
+ if pos('cells_per_wl :',s)=1 then begin
+ delete(s,1,pos(':',s));
+ while (length(s)>0) and (s[1]=' ') do
+ delete(s,1,1);
+ dx:=1/strtofloat(s);
+ end;
+ if pos('cells_left :',s)=1 then begin
+ delete(s,1,pos(':',s));
+ while (length(s)>0) and (s[1]=' ') do
+ delete(s,1,1);
+ cells_left:=strtofloat(s);
+ end;
+ if s='domain: particles' then
+ j:=2;
+ end;
+ end{of case};
+ end;
+ closefile(tf);
+ if (not odd(i)) or ((absShift<-1.5e9) and ((dx=0) or (cells_left=0))) then Fehler('Unerwartetes Dateiende in '''+nam+'''!');
+
+ setlength(d1,0);
+ setlength(d2,0);
+ i32:=0;
+ fl:=0;
+ start:=-1;
+ stop:=-1;
+ i:=findFirst(nam+'*',$3F,sr);
+ while i=0 do begin
+ s:=sr.Name;
+ while pos('-',s)>0 do
+ delete(s,1,pos('-',s));
+ i:=strtoint(s);
+ if (start=-1) or (start>i) then
+ start:=i;
+ if (stop=-1) or (stop<i) then
+ stop:=i;
+ i:=findNext(sr);
+ end;
+ findClose(sr);
+ if start=-1 then
+ Fehler('Keine Dateien '''+nam+'*'' gefunden!');
+ for i:=start to stop do
+ if not fileexists(nam+inttostr(i)) then
+ Fehler('Die Dateien beginnend mit '''+nam+''' bilden keine vollständige Abfolge! ('+inttostr(i)+' fehlt)');
+
+ for i:=start to stop do begin
+ write(stderr,'Datei '''+nam+inttostr(i)+''' einlesen ');
+ assignfile(f,nam+inttostr(i));
+ reset(f,1);
+ blockread(f,i32,sizeof(int32)); // Dateinummer
+ if i32<>i then begin
+ closefile(f);
+ Fehler('Die Datei '''+nam+inttostr(i)+''' behauptet, sie sei die '+inttostr(i32)+'-te!');
+ end;
+
+ blockread(f,i32,sizeof(int32)); // Anzahl x-Positionen
+ traces:=i32;
+ blockread(f,i32,sizeof(int32)); // Anzahl t-Positionen
+ steps:=i32;
+
+ setlength(buff,steps);
+ setlength(d1,length(d1)+steps);
+ setlength(d2,length(d2)+steps);
+
+ for j:=0 to traces-1 do begin
+ blockread(f,fl,sizeof(single)); // x-Position
+ if j=0 then begin
+ if i=start then cells_left:=cells_left-fl;
+ blockread(f,buff[0],sizeof(single)*steps); // fp
+ for k:=0 to steps-1 do begin
+ d1[k+length(d1)-steps].x:=(i-1+k/steps)*factor;
+ d1[k+length(d1)-steps].y:=buff[k];
+ end;
+ blockread(f,buff[0],sizeof(single)*steps); // fm
+ for k:=0 to steps-1 do begin
+ d2[k+length(d2)-steps].x:=(i-1+k/steps)*factor;
+ d2[k+length(d2)-steps].y:=buff[k];
+ end;
+ end;
+ for k:=0 to 7+2*byte(j<>0) do // (fp,fm,)gp,gm,ex,de,di,jx,jy,jz
+ blockread(f,buff[0],sizeof(single)*steps);
+ end;
+
+ if not eof(f) then
+ Fehler('Die Datei '''+nam+inttostr(i)+''' ist zu groß!');
+
+ closefile(f);
+ writeln(stderr,'fertig');
+ end;
+ if absShift<-1.5e9 then begin
+ absShift:=dx*cells_left*factor*2;
+ writeln(stderr,'Aus Inputdatei ausgelesene konstante Verschiebung des Outputpulses: '+floattostr(absShift)+' T');
+ end;
+end;
+
+procedure readTextInput(nam: string; out dat: tExtPointArray);
+var f: textfile;
+ i: longint;
+ s: string;
+begin
+ if not fileexists(nam) then
+ Fehler('Datei '''+nam+''' existiert nicht!');
+ write(stderr,'Datei '''+nam+''' einlesen ');
+ assignfile(f,nam);
+ reset(f);
+ i:=0;
+ while not eof(f) do begin
+ readln(f,s);
+ if pos('#',s)>0 then delete(s,pos('#',s),length(s));
+ while (length(s)>0) and (s[1]=' ') do delete(s,1,1);
+ if s='' then continue;
+ inc(i);
+ end;
+ closefile(f);
+ reset(f);
+ setlength(dat,i);
+ for i:=0 to length(dat)-1 do begin
+ readln(f,s);
+ if pos('#',s)>0 then delete(s,pos('#',s),length(s));
+ while (length(s)>0) and (s[1]=' ') do delete(s,1,1);
+ if s='' then continue;
+ dat[i].x:=strtofloat(copy(s,1,pos(' ',s)-1));
+ delete(s,1,pos(' ',s));
+ while (length(s)>0) and (s[1]=' ') do delete(s,1,1);
+ dat[i].y:=strtofloat(copy(s,1,pos(' ',s+' ')-1));
+ if i and 1023 = 0 then write(stderr,'.');
+ end;
+ closefile(f);
+ writeln(stderr,'fertig');
+end;
+
+procedure writeOutput(nam: string; const dat: tExtPointArray);
+var f: textfile;
+ i,ml,mr: longint;
+ leer: string;
+ strings: array of array[boolean] of string;
+begin
+ if length(dat)=0 then exit;
+ write(stderr,'Datei '''+nam+''' schreiben ');
+ assignfile(f,nam);
+ rewrite(f);
+ setlength(strings,length(dat));
+ ml:=0;
+ mr:=0;
+ for i:=0 to length(dat)-1 do begin
+ strings[i,false]:=floattostr(dat[i].x);
+ strings[i,true]:=floattostr(dat[i].y);
+ ml:=max(ml,length(strings[i,false]));
+ mr:=max(mr,pos('.',strings[i,true]+'.'));
+ if i and 65535 = 0 then write(stderr,'.');
+ end;
+ setlength(leer,ml+mr);
+ fillchar(leer[1],length(leer),' ');
+ for i:=0 to length(dat)-1 do begin
+ writeln(f,strings[i,false]+copy(leer,1,1+ml+mr-length(strings[i,false])-pos('.',strings[i,true]+'.'))+strings[i,true]);
+ if i and 65535 = 0 then write(stderr,'.');
+ end;
+ closefile(f);
+ writeln(stderr,' fertig');
+end;
+
+procedure berechneTrajektorie(const inPuls,outPuls: tExtPointArray; out surTraj: tExtPointArray; absShift: extended);
+var i,j: longint;
+ exList: array[0..3] of tLongintArray;
+ psThreads: array of tPhaseShiftThread;
+ stuecke: array of array[boolean,boolean] of longint;
+ erledigt: array of boolean;
+ b1,b2: boolean;
+begin
+ for i:=0 to length(inPuls)-2 do
+ if inPuls[i+1].x<=inPuls[i].x then
+ writeln('Der Inpuls ist nicht sortiert!');
+ for i:=0 to length(outPuls)-2 do
+ if outPuls[i+1].x<=outPuls[i].x then
+ writeln('Der Outpuls ist nicht sortiert!');
+
+ // Extrema identifizieren
+ exList[0]:=findeExtrema(inPuls,false);
+ exList[1]:=findeExtrema(inPuls,true);
+ exList[2]:=findeExtrema(outPuls,false);
+ exList[3]:=findeExtrema(outPuls,true);
+
+ // und filtern
+ filtereExtrema(inPuls,exList[0],1,2);
+ filtereExtrema(inPuls,exList[1],1,2);
+ filtereExtrema(outPuls,exList[2],1,2);
+ filtereExtrema(outPuls,exList[3],1,2);
+
+ monotonieHerstellen(inPuls,exList[0],exList[1]);
+ monotonieHerstellen(outPuls,exList[2],exList[3]);
+
+ writeln;
+ for i:=0 to 3 do begin
+ for j:=0 to length(exList[i])-1 do
+ write(exList[i,j],'(',j,') ');
+ writeln;
+ end;
+
+ vereineExtrema(inPuls,outPuls,exList[0],exList[2],0.5);
+ vereineExtrema(inPuls,outPuls,exList[1],exList[3],0.5);
+
+writeln(length(exList[0]),' ',length(exList[1]),' ',length(exList[2]),' ',length(exList[3]));
+
+{ for k:=0 to 1 do begin
+ i:=0;
+ while (i<length(exList[k])-1) and (exList[k,i+1]<exList[1-k,0]) do
+ inc(i);
+ for j:=i to length(exList[k])-1 do
+ exList[k,j-i]:=exList[k,j];
+ setlength(exList[k],length(exList[k])-i);
+ for j:=i to length(exList[2+k])-1 do
+ exList[2+k,j-i]:=exList[2+k,j];
+ setlength(exList[2+k],length(exList[2+k])-i);
+
+ i:=length(exList[k])-1;
+ while (i>0) and (exList[k,i-1]>exList[1-k,length(exList[1-k])-1]) do
+ dec(i);
+ setlength(exList[k],i+1);
+ setlength(exList[2+k],i+1);
+ end; }
+
+ writeln;
+ for i:=0 to 3 do begin
+ for j:=0 to length(exList[i])-1 do
+ write(exList[i,j],'(',j,') ');
+ writeln;
+ end;
+
+ setlength(stuecke,length(exList[0])+length(exList[1])-1);
+ writeln(length(stuecke),' ',length(exList[0]),' ',length(exList[1]),' ',length(exList[2]),' ',length(exList[3]));
+ for i:=0 to length(stuecke)-1 do
+ for b1:=false to true do // Anfang / Ende
+ for b2:=false to true do // input / output
+ stuecke[i,b1,b2]:=
+ exList[(i and 1) xor byte(exList[0,0]>exList[1,0]) xor (2*byte(b2)) xor byte(b1),(i+byte(b1)) div 2];
+
+ writeln;
+ for i:=0 to length(stuecke)-1 do
+ writeln(i,': ',stuecke[i,false,false],' ',stuecke[i,false,true],' ',stuecke[i,true,false],' ',stuecke[i,true,true]);
+ writeln;
+{ for i:=0 to length(stuecke)-1 do
+ writeln(inPuls[stuecke[i,false,false]].x,' ',outPuls[stuecke[i,false,true]].x,' ',inPuls[stuecke[i,true,false]].x,' ',outPuls[stuecke[i,true,true]].x);
+ writeln; }
+
+ setlength(erledigt,length(stuecke));
+ for i:=0 to length(erledigt)-1 do
+ erledigt[i]:=false;
+ setlength(surTraj,0);
+ setlength(psThreads,8);
+ for i:=0 to length(psThreads)-1 do
+ psThreads[i]:=nil;
+ b1:=true;
+ repeat
+ b1:=true;
+ for i:=0 to length(psThreads)-1 do begin
+ if (psThreads[i]<>nil) and (psThreads[i].fertig) then begin
+ setlength(surTraj,length(surTraj)+length(psThreads[i]._ergebnis));
+ for j:=0 to length(psThreads[i]._ergebnis)-1 do
+ surTraj[length(surTraj)-1-j]:=psThreads[i]._ergebnis[j];
+ psThreads[i].free;
+ psThreads[i]:=nil;
+ end;
+ if psThreads[i]=nil then begin
+ j:=0;
+ while (j<length(erledigt)) and erledigt[j] do
+ inc(j);
+ if j>=length(erledigt) then continue;
+ psThreads[i]:=
+ tPhaseShiftThread.create(
+ @inPuls,
+ @outPuls,
+ stuecke[j,false,false],
+ stuecke[j,true,false],
+ stuecke[j,false,true],
+ stuecke[j,true,true],
+ absShift);
+ psThreads[i].Suspended:=false;
+ erledigt[j]:=true;
+ end;
+ b1:=false;
+ end;
+ if not b1 then sleep(100);
+ for i:=0 to length(erledigt)-1 do
+ b1:=b1 and erledigt[i];
+ until b1;
+
+
+(* setlength(surTraj,min(length(inPuls),length(outPuls)));
+ for i:=0 to length(surTraj)-1 do begin
+ surTraj[i].x:=(inPuls[i].x+outPuls[i].x)/2;
+ surTraj[i].y:=inPuls[i].y+outPuls[i].y;
+ end;
+ for i:=0 to length(exList[0])-1 do
+ surTraj[exList[0,i]].y:=-1;
+ for i:=0 to length(exList[1])-1 do
+ surTraj[exList[1,i]].y:=1;
+ for i:=0 to length(exList[2])-1 do
+ surTraj[exList[2,i]].y:=0;
+ for i:=0 to length(exList[3])-1 do
+ surTraj[exList[3,i]].y:=0; *)
+end;
+
+procedure smoothen(var dat: tExtPointArray; width: longint);
+var i,j: longint;
+begin
+ for i:=0 to length(dat)-width do begin
+ for j:=1 to width-1 do
+ dat[i]:=Plus(dat[i],dat[i+j]);
+ dat[i]:=Durch(dat[i],width);
+ end;
+ setlength(dat,length(dat)-width+1);
+end;
+
+procedure sort(var dat: tExtPointArray);
+begin
+ sort(dat,0,length(dat)-1,8);
+end;
+
+procedure sort(var dat: tExtPointArray; start, stop, threads: longint);
+var sortThreads: array of tSortThread;
+ i,j,k: longint;
+ fertig: boolean;
+ apos,epos: array of longint;
+ tmp: tExtPointArray;
+begin
+ if start>=stop then exit;
+ if threads=1 then begin
+ sort(dat,start,(start+stop) div 2,threads);
+ sort(dat,(start+stop) div 2 + 1,stop,threads);
+ setlength(apos,2);
+ setlength(epos,2);
+ apos[0]:=start;
+ apos[1]:=(start+stop) div 2 + 1;
+ epos[0]:=(start+stop) div 2;
+ epos[1]:=stop;
+ end
+ else begin
+ threads:=min(threads,stop-start+1);
+ setlength(sortThreads,threads);
+ for i:=0 to length(sortThreads)-1 do
+ sortThreads[i]:=
+ tSortThread.create(
+ @dat,
+ start+round(i * (stop+1-start)/threads),
+ start+round((i+1)*(stop+1-start)/threads)-1);
+ for i:=0 to length(sortThreads)-1 do
+ sortThreads[i].Suspended:=false;
+ repeat
+ fertig:=true;
+ sleep(100);
+ for i:=0 to length(sortThreads)-1 do
+ fertig:=fertig and sortThreads[i].fertig;
+ until fertig;
+ for i:=0 to length(sortThreads)-1 do
+ sortThreads[i].free;
+ setlength(sortThreads,0);
+ setlength(apos,threads);
+ setlength(epos,threads);
+ for i:=0 to threads-1 do begin
+ apos[i]:=start+round(i * (stop+1-start)/threads);
+ epos[i]:=start+round((i+1)*(stop+1-start)/threads)-1;
+ end;
+ end;
+ // und jetzt nur noch mergen
+ setlength(tmp,stop-start+1);
+ k:=0;
+ repeat
+ fertig:=true;
+ j:=-1;
+ for i:=0 to length(apos)-1 do begin
+ if apos[i]>epos[i] then continue;
+ fertig:=false;
+ if (j=-1) or (dat[apos[i]].x<dat[apos[j]].x) then j:=i;
+ end;
+ if not fertig then begin
+ tmp[k]:=dat[apos[j]];
+ inc(k);
+ inc(apos[j]);
+ end;
+ until fertig;
+ for i:=start to stop do
+ dat[i]:=tmp[i-start];
+ setlength(tmp,0);
+end;
+
+procedure cut(var dat: tExtPointArray; tmax: extended);
+var i: longint;
+begin
+ if tmax<0 then exit;
+ i:=0;
+ while (i<length(dat)) and (dat[i].x<tmax) do
+ inc(i);
+ setlength(dat,i);
+end;
+
+function findeExtrema(const dat: tExtPointArray; maxima: boolean): tLongintArray;
+var i,j: longint;
+begin
+ j:=0;
+ for i:=1 to length(dat)-2 do
+ if (maxima and (dat[i].y>dat[i-1].y) and (dat[i].y>=dat[i+1].y) and (dat[i].y>0)) or
+ ((not maxima) and (dat[i].y<dat[i-1].y) and (dat[i].y<=dat[i+1].y) and (dat[i].y<0)) then
+ inc(j);
+ setlength(result,j);
+ j:=0;
+ for i:=1 to length(dat)-2 do
+ if (maxima and (dat[i].y>dat[i-1].y) and (dat[i].y>=dat[i+1].y) and (dat[i].y>0)) or
+ ((not maxima) and (dat[i].y<dat[i-1].y) and (dat[i].y<=dat[i+1].y) and (dat[i].y<0)) then begin
+ result[j]:=i;
+ inc(j);
+ end;
+end;
+
+procedure filtereExtrema(const dat: tExtPointArray; var extrema: tLongintArray; sollAbst, toleranz: extended);
+var mitte,i,j,k: longint;
+ mx,dist: extended;
+ behalten: array of boolean;
+ Toepfe: array of longint;
+begin
+ setlength(behalten,length(extrema));
+ for i:=0 to length(behalten)-1 do
+ behalten[i]:=false;
+
+ setlength(toepfe,round(sqrt(length(extrema))));
+ for i:=0 to length(toepfe)-1 do
+ toepfe[i]:=0;
+
+ for i:=0 to length(extrema)-1 do
+ inc(toepfe[round(dat[extrema[i]].x/sollAbst*length(toepfe)) mod length(toepfe)]);
+
+ mitte:=0;
+ for i:=1 to length(toepfe)-1 do
+ if toepfe[i]>toepfe[mitte] then
+ mitte:=i;
+
+ mx:=0;
+ for i:=0 to length(extrema)-1 do
+ mx:=mx + myfrac(dat[extrema[i]].x - mitte/length(toepfe) + 0.5);
+
+ mx:=(mx/length(extrema) + mitte/length(toepfe) - 0.5)*sollAbst;
+
+ for i:=round((dat[extrema[0]].x - mx)/sollAbst) to round((dat[extrema[length(extrema)-1]].x - mx)/sollAbst) do begin
+ dist:=abs(dat[extrema[0]].x - i*sollAbst - mx);
+ k:=0;
+ for j:=1 to length(behalten)-1 do
+ if abs(dat[extrema[j]].x - i*sollAbst - mx) < dist then begin
+ dist:=abs(dat[extrema[j]].x - i*sollAbst - mx);
+ k:=j;
+ end;
+ if 2*dist<toleranz then
+ behalten[k]:=true;
+ end;
+
+(*// Variante 1
+ mx:=(dat[0].x+dat[length(dat)-1].x)/2;
+ mitte:=0;
+ for i:=1 to length(extrema)-1 do
+ if abs(dat[extrema[mitte]].x-mx)>abs(dat[extrema[i]].x-mx) then
+ mitte:=i;
+ repeat
+ j:=mitte;
+ repeat
+ k:=-1;
+ for i:=j-1 downto 0 do
+ if (abs(dat[extrema[j]].x-dat[extrema[i]].x-sollAbst)<=toleranz) and
+ ((k=-1) or (abs(dat[extrema[j]].x-dat[extrema[i]].x-sollAbst)<abs(dat[extrema[j]].x-dat[extrema[k]].x-sollAbst))) then
+ k:=i;
+ if k<>-1 then begin
+ behalten[k]:=true;
+ behalten[j]:=true;
+ j:=k;
+ end;
+ until k=-1;
+ if j=mitte then
+ inc(mitte);
+ until (j<>mitte) or (mitte>length(extrema));
+ if j=mitte then Fehler('So, wie sich der Programmierer das gedacht hat, findet man die richtigen Maxima nicht ...');
+ j:=mitte;
+ repeat
+ k:=-1;
+ for i:=j+1 to length(extrema)-1 do
+ if (abs(dat[extrema[i]].x-dat[extrema[j]].x-sollAbst)<=toleranz) and
+ ((k=-1) or (abs(dat[extrema[i]].x-dat[extrema[j]].x-sollAbst)<abs(dat[extrema[k]].x-dat[extrema[j]].x-sollAbst))) then
+ k:=i;
+ if k<>-1 then begin
+ behalten[k]:=true;
+ behalten[j]:=true;
+ j:=k;
+ end;
+ until k=-1;
+ if j=mitte then Fehler('So, wie sich der Programmierer das gedacht hat, findet man die richtigen Maxima nicht ...');
+*)
+
+ j:=0;
+ for i:=0 to length(behalten)-1 do
+ if behalten[i] then begin
+ extrema[j]:=extrema[i];
+// write(extrema[j],' ');
+ inc(j);
+ end;
+ setlength(extrema,j);
+// writeln;
+end;
+
+procedure monotonieHerstellen(const dat: tExtPointArray; var minima,maxima: tLongintArray);
+var i,j,k: longint;
+ b: boolean;
+ behalten: array[boolean] of array of Boolean;
+begin
+ setlength(behalten[false],length(maxima));
+ setlength(behalten[true],length(minima));
+ for b:=false to true do
+ for i:=0 to length(behalten[b])-1 do
+ behalten[b,i]:=false;
+
+ i:=0;
+ j:=0;
+ b:=maxima[0]>minima[0];
+ repeat
+ if b then begin
+ k:=j;
+ repeat
+ if (j<length(minima)) and (dat[minima[j]].y<dat[minima[k]].y) then
+ k:=j;
+ inc(j);
+ until (j>=length(minima)) or ((i<length(maxima)) and (maxima[i]<minima[j]));
+ behalten[true,k]:=true;
+ b:=false;
+ end
+ else begin
+ k:=i;
+ repeat
+ if (i<length(maxima)) and (dat[maxima[i]].y>dat[maxima[k]].y) then
+ k:=i;
+ inc(i);
+ until (i>=length(maxima)) or ((j<length(minima)) and (maxima[i]>=minima[j]));
+ behalten[false,k]:=true;
+ b:=true;
+ end;
+ until (i>=length(maxima)) and (j>=length(minima));
+
+ i:=0;
+ j:=0;
+ while i<length(maxima) do begin
+ if behalten[false,i] then begin
+ maxima[j]:=maxima[i];
+ inc(j);
+ end;
+ inc(i);
+ end;
+ setlength(maxima,j);
+
+ i:=0;
+ j:=0;
+ while i<length(minima) do begin
+ if behalten[true,i] then begin
+ minima[j]:=minima[i];
+ inc(j);
+ end;
+ inc(i);
+ end;
+ setlength(minima,j); writeln(length(maxima),' ',length(minima));
+end;
+
+procedure gesamtverschiebung(var inPuls,outPuls: tExtPointArray; var absShift: extended);
+var iMax,oMax,i: longint;
+begin
+ iMax:=0;
+ for i:=1 to length(inPuls)-1 do
+ if inPuls[i].y>inPuls[iMax].y then iMax:=i;
+ oMax:=0;
+ for i:=1 to length(outPuls)-1 do
+ if outPuls[i].y>outPuls[oMax].y then oMax:=i;
+ if absShift<-0.9e9 then absShift:=outPuls[oMax].x-inPuls[iMax].x
+ else absShift:=(outPuls[oMax].x-inPuls[iMax].x)-round((outPuls[oMax].x-inPuls[iMax].x)-absShift);
+
+ for i:=0 to length(outPuls)-1 do
+ outPuls[i].x:=outPuls[i].x-absShift;
+ oMax:=0;
+ while (oMax<length(outPuls)) and (outPuls[oMax].x<inPuls[0].x) do
+ inc(oMax);
+ for i:=oMax to length(outPuls)-1 do
+ outPuls[i-oMax]:=outPuls[i];
+ setlength(outPuls,length(outPuls)-oMax);
+ cut(inPuls,outPuls[length(outPuls)-1].x);
+ writeln(stderr,'Die konstante Verschiebung wurde auf '+floattostr(absShift)+' T optimiert.');
+end;
+
+procedure integrate(var dat: tExtPointArray);
+var i: longint;
+begin
+ for i:=1 to length(dat)-1 do
+ dat[i].y:=dat[i].y * (dat[i].x-dat[i-1].x)+dat[i-1].y;
+end;
+
+procedure removeLinearOffset(var dat: tExtPointArray);
+var i: longint;
+ dx,dy: extended;
+begin
+ dx:=0;
+ dy:=0;
+ for i:=0 to 99 do begin
+ dy:=dy + dat[length(dat)-1-i].y - dat[i].y;
+ dx:=dx + dat[length(dat)-1-i].x - dat[i].x;
+ end;
+ for i:=0 to length(dat)-1 do
+ dat[i].y:=dat[i].y - dy/dx*(dat[i].x-dat[0].x);
+ dy:=0;
+ for i:=0 to 99 do
+ dy:=dy + dat[length(dat)-1-i].y + dat[i].y;
+ dy:=dy/200;
+ for i:=0 to length(dat)-1 do
+ dat[i].y:=dat[i].y - dy;
+end;
+
+procedure flip(var dat: tExtPointArray);
+var i: longint;
+begin
+ for i:=0 to length(dat)-1 do
+ dat[i].y:=-dat[i].y;
+end;
+
+procedure vereineExtrema(const dat1,dat2: tExtPointArray; var el1,el2: tLongintArray; toleranz: extended);
+var i,j,k: longint;
+ b: boolean;
+ behalten: array[boolean] of array of boolean;
+begin
+ setlength(behalten[false],length(el1));
+ setlength(behalten[true],length(el2));
+ for b:=false to true do
+ for i:=0 to length(behalten[b])-1 do
+ behalten[b,i]:=false;
+
+(* i:=0;
+ j:=0;
+ b:=dat1[el1[0]].x > dat2[el2[0]].x;
+ repeat
+ if b then begin
+ behalten[true,i]:=true;
+ while (i+1<length(el1)) and ((j>=length(el2)) or (dat1[el1[i+1]].x <= dat2[el2[j]].x)) do
+ inc(i);
+ inc(i);
+ b:=false;
+ end
+ else begin
+ behalten[false,j]:=true;
+ while (j+1<length(el2)) and ((i>=length(el1)) or (dat1[el1[i]].x > dat2[el2[j+1]].x)) do
+ inc(j);
+ inc(j);
+ b:=true;
+ end;
+ until (i>=length(el1)) and (j>=length(el2)); *)
+
+ for i:=0 to length(el1)-1 do begin
+ k:=-1;
+ for j:=0 to length(el2)-1 do
+ if (abs(dat1[el1[i]].x-dat2[el2[j]].x)<=toleranz) and
+ ((k=-1) or (abs(dat1[el1[i]].x-dat2[el2[j]].x)<abs(dat1[el1[i]].x-dat2[el2[k]].x))) then
+ k:=j;
+ if (k<>-1) and not behalten[true,k] then begin
+ behalten[false,i]:=true;
+ behalten[true,k]:=true;
+ end;
+ end;
+
+ j:=0;
+ for i:=0 to length(behalten[false])-1 do
+ if behalten[false,i] then begin
+ el1[j]:=el1[i];
+ inc(j);
+ end;
+ setlength(el1,j);
+ j:=0;
+ for i:=0 to length(behalten[true])-1 do
+ if behalten[true,i] then begin
+ el2[j]:=el2[i];
+ inc(j);
+ end;
+ setlength(el2,j);
+end;
+
+procedure uniq(var dat: tExtPointArray; streng: boolean);
+var i,j,k: longint;
+ tmp: extended;
+ s: string;
+begin
+ write(stderr,'Dopplungen entfernen ...');
+ j:=0;
+ k:=0;
+ tmp:=0;
+ for i:=0 to length(dat)-1 do begin
+ if (i<length(dat)-1) and (dat[i].x=dat[i+1].x) and ((dat[i].y=dat[i+1].y) or not streng) then begin
+ inc(k);
+ tmp:=tmp+dat[i].y;
+ continue;
+ end;
+ if (k>0) then begin
+ dat[j].y:=(tmp+dat[i].y)/(k+1);
+ dat[j].x:=dat[i].x;
+ k:=0;
+ tmp:=0;
+ end
+ else dat[j]:=dat[i];
+ inc(j);
+ end;
+ s:=floattostr(round(100*(length(dat)-j)/length(dat)*10)/10);
+ if pos('.',s)>0 then delete(s,pos('.',s)+2,length(s));
+ writeln(stderr,' fertig (es gab '+inttostr(length(dat)-j)+' Dopplungen - etwa '+s+'%)');
+ setlength(dat,j);
+end;
+
+procedure fft(var dat: tExtPointArray);
+var i,j,k,n,dist,absch,wnum,wstep,haL: longint;
+ in0,out0: boolean;
+ wRe,wIm: tExtendedArray;
+ t1,t2,vorher,nachher,fstep,pvFehler: extended;
+ umsortierung: tLongintArray;
+begin
+ write(stderr,'FFT ... ');
+ fstep:=dat[1].x-dat[0].x;
+ j:=length(dat);
+ setlength(dat,2*round(power(2,ceil(ln(length(dat))/ln(2)))));
+ k:=(length(dat)-j) div 2;
+ for i:=j-1 downto 0 do
+ dat[i+k]:=dat[i];
+ for i:=0 to k-1 do // weich auf 0 abklingen lassen
+ dat[i].y:=dat[k].y * sqr(sin(pi/2*i/k));
+ for i:=j+k to length(dat)-1 do // dito
+ dat[i].y:=dat[j+k-1].y * sqr(sin(pi/2*(length(dat)-1-i)/(length(dat)-j-k)));
+ for i:=0 to length(dat)-1 do begin
+ dat[i].x:=dat[i].y;
+ dat[i].y:=0;
+ end;
+
+ haL:=length(dat) div 2;
+ n:=round(ln(2*hal)/ln(2));
+
+ fstep:=1/(2*haL)/fstep;
+
+ vorher:=0;
+ for i:=0 to 2*haL-1 do
+ vorher:=vorher + sqr(dat[i].x);
+
+ setlength(umsortierung,2*haL);
+ for j:=0 to 2*haL-1 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(wRe,haL);
+ setlength(wIm,haL);
+ for i:=0 to haL-1 do begin
+ wRe[i]:=cos(pi*i/haL);
+ wIm[i]:=sin(pi*i/haL);
+ end;
+
+ in0:=true;
+ out0:=true;
+ for i:=0 to 2*haL-1 do begin
+ if umsortierung[i]>i then begin
+ t1:=dat[i].x;
+ dat[i].x:=dat[umsortierung[i]].x;
+ dat[umsortierung[i]].x:=t1;
+ end;
+ in0:=in0 and (dat[i].x=0);
+ end;
+
+ dist:=1;
+ wstep:=haL;
+ while dist<2*haL do begin
+ absch:=0;
+ while absch<2*haL do begin
+ wnum:=0;
+ for i:=0 to dist-1 do begin
+ // x_links: [absch+j]
+ // x_rechts: [absch+j+dist]
+ t1:=wRe[wnum]*dat[absch+i+dist].x - wIm[wnum]*dat[absch+i+dist].y;
+ t2:=wRe[wnum]*dat[absch+i+dist].y + wIm[wnum]*dat[absch+i+dist].x;
+
+ dat[absch+i+dist].x:=dat[absch+i].x-t1;
+ dat[absch+i+dist].y:=dat[absch+i].y-t2;
+ dat[absch+i].x:=dat[absch+i].x+t1;
+ dat[absch+i].y:=dat[absch+i].y+t2;
+ wnum:=wnum+wstep;
+ end;
+ absch:=absch+2*dist;
+ end;
+ wstep:=wstep div 2;
+ dist:=dist*2;
+ end;
+
+ for i:=0 to 2*haL-1 do begin
+ dat[i].y:=(sqr(dat[i].x)+sqr(dat[i].y))/(2*haL);
+ dat[i].x:=fStep*i;
+ end;
+ for j:=0 to 2*haL-1 do
+ out0:=out0 and (dat[i].y=0);
+
+ nachher:=0;
+ for i:=0 to 2*haL-1 do
+ nachher:=nachher + dat[i].y;
+
+ if (nachher=0) and (vorher=0) then pvFehler:=0
+ else pvFehler:=abs(nachher-vorher)/(nachher+vorher);
+
+ writeln(stderr,' fertig (Parseval-Fehler = '+floattostr(pvFehler)+' ... nämlich von '+floattostr(vorher)+' zu '+floattostr(nachher)+')');
+ if in0 then writeln(stderr,'Nur Nullen im Input der FFT!');
+ if out0 then writeln(stderr,'Nur Nullen im Output der FFT!');
+end;
+
+procedure interpoliere(var dat: tExtPointArray);
+var i,j: longint;
+ tmp: extended;
+ tdat: tExtPointArray;
+begin
+ tmp:=dat[1].x-dat[0].x;
+ for i:=2 to length(dat)-1 do
+ tmp:=min(tmp,dat[i].x-dat[i-1].x);
+ if tmp<=0 then Fehler('Die Daten müssen sortiert sein und dürfen keine doppelten x-Werte enthalten! ('+floattostr(tmp)+')');
+ setlength(tdat,length(dat));
+ for i:=0 to length(dat)-1 do
+ tdat[i]:=dat[i];
+ setlength(dat,max(2*length(tdat),round(min(power(length(tdat),1.3),(tdat[length(tdat)-1].x-tdat[0].x)/tmp+1))));
+ j:=0;
+ for i:=0 to length(dat)-1 do begin
+ dat[i].x:=tdat[0].x+(tdat[length(tdat)-1].x-tdat[0].x)*i/(length(dat)-1);
+ while (j<length(tdat)-1) and (tdat[j+1].x <= dat[i].x) do
+ inc(j);
+ if j=length(tdat)-1 then tmp:=tdat[j].y
+ else begin
+ tmp:=(dat[i].x-tdat[j].x)/(tdat[j+1].x-tdat[j].x);
+ tmp:=tdat[j].y*(1-tmp) + tdat[j+1].y*tmp;
+ end;
+ dat[i].y:=tmp;
+ end;
+end;
+
+procedure normiere(var dat: tExtPointArray);
+var i: longint;
+ m: extended;
+begin
+ if length(dat)=0 then exit;
+ m:=dat[0].y;
+ for i:=1 to length(dat)-1 do
+ m:=max(m,dat[i].y);
+ if m=0 then exit;
+ for i:=0 to length(dat)-1 do
+ dat[i].y:=dat[i].y/m;
+end;
+
+// tSortThread *****************************************************************
+
+constructor tSortThread.create(pd: pTExtPointArray; sta, sto: longint);
+begin
+ inherited create(true);
+ _pdat:=pd;
+ _start:=sta;
+ _stop:=sto;
+ fertig:=false;
+end;
+
+destructor tSortThread.destroy;
+begin
+ inherited destroy;
+end;
+
+procedure tSortThread.execute;
+begin
+ sort(_pdat^,_start,_stop,1);
+ fertig:=true;
+end;
+
+// tPhaseShiftThread ***********************************************************
+
+constructor tPhaseShiftThread.create(p1d,p2d: pTExtPointArray; sta1, sto1, sta2, sto2: longint; absShift: extended);
+var i: longint;
+ b: boolean;
+begin
+ inherited create(true);
+ _starts[false]:=sta1;
+ _stops[false]:=sto1;
+ _starts[true]:=sta2;
+ _stops[true]:=sto2;
+ _absShift:=absShift;
+ for b:=false to true do
+ setlength(_dats[b],_stops[b]-_starts[b]+1);
+ for i:=0 to length(_dats[false])-1 do
+ _dats[false,i]:=p1d^[i+_starts[false]];
+ for i:=0 to length(_dats[true])-1 do
+ _dats[true,i]:=p2d^[i+_starts[true]];
+ setlength(_ergebnis,2+sto1+sto2-sta1-sta2);
+ fertig:=false;
+end;
+
+destructor tPhaseShiftThread.destroy;
+begin
+ setlength(_ergebnis,0);
+ inherited destroy;
+end;
+
+procedure tPhaseShiftThread.execute;
+var i,j,offset: longint;
+ b: boolean;
+ t1,t2,t3,bestErr,bestPos: extended;
+begin
+ for b:=false to true do begin // normieren (auf 0..1)
+ t1:=0;
+ t2:=0;
+ for i:=0 to length(_dats[b])-1 do begin
+ t1:=max(t1,_dats[b,i].y);
+ t2:=min(t2,_dats[b,i].y);
+ end;
+ for i:=0 to length(_dats[b])-1 do
+ _dats[b,i].y:=(_dats[b,i].y-t2)/(t1-t2);
+ end;
+
+ for b:=false to true do begin
+ offset:=byte(b)*(_stops[false]-_starts[false]+1);
+ for i:=0 to _stops[b]-_starts[b] do begin
+ _ergebnis[offset + i].x:=_dats[b,i].x;
+ if (i=0) then bestPos:=_dats[not b,0].x
+ else begin
+ t1:=_dats[b,i].y;
+ j:=0; // floor(_ergebnis[offset + i - 1].y); <- das war Murks
+(* while (j<_stops[not b]-_starts[not b]) and
+ (_dats[not b,j].x<_ergebnis[offset + i - 1].y) do
+ inc(j);*)
+ bestPos:=-1;
+ bestErr:=-1;
+ repeat
+ while (j<_stops[not b]-_starts[not b]) and
+ ((_dats[not b,j].y - t1) * (_dats[not b,j + 1].y - t1) > 0) do
+ inc(j);
+ t2:=_dats[not b,j].x;
+ if j<_stops[not b]-_starts[not b] then // x1 = dx * y1 / (y1 - y2) ... Berechnung des gebrochenen Anteils
+ t2:=t2 + (_dats[not b,j+1].x-_dats[not b,j].x) * (_dats[not b,j].y - t1) / (_dats[not b,j].y - _dats[not b,j+1].y);
+ t3:=abs(t2-_ergebnis[offset + i].x-_ergebnis[offset + i - 1].y+_ergebnis[offset + i - 1].x);
+ if (bestPos<0) or (bestErr>t3) then begin
+ bestErr:=t3;
+ bestPos:=t2;
+ end;
+ inc(j);
+ until j>_stops[not b]-_starts[not b];
+ end;
+ _ergebnis[offset + i].y:=bestPos;
+ end;
+ end;
+ for i:=_stops[false]-_starts[false]+1 to length(_ergebnis)-1 do begin
+ t1:=_ergebnis[i].x;
+ _ergebnis[i].x:=_ergebnis[i].y;
+ _ergebnis[i].y:=t1;
+ end;
+ for i:=0 to length(_ergebnis)-1 do begin
+ _ergebnis[i].y:=_ergebnis[i].y+_absShift;
+ t1:=(_ergebnis[i].y+_ergebnis[i].x)/2;
+ _ergebnis[i].y:=(_ergebnis[i].y-_ergebnis[i].x)/2;
+ _ergebnis[i].x:=t1;
+ end;
+ fertig:=true;
+end;
+
+end.
+
generated by cgit v1.2.3-70-g09d2 (git 2.47.0) at 2024-12-04 08:21:26 +0000