From 17ec74a07cacd08b0904b29a390d3b7f26ba7b1b Mon Sep 17 00:00:00 2001 From: Erich Eckner Date: Tue, 23 Sep 2014 16:12:24 +0200 Subject: Initialer Commit --- .gitignore | 13 + Bild.pdf | Bin 0 -> 74662 bytes ROM.lpi | 96 ++++++ ROM.lpr | 150 ++++++++ ROM.lps | 166 +++++++++ make.gnu | 26 ++ mathunit.pas | 46 +++ romunit.pas | 1086 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 8 files changed, 1583 insertions(+) create mode 100644 .gitignore create mode 100644 Bild.pdf create mode 100644 ROM.lpi create mode 100644 ROM.lpr create mode 100644 ROM.lps create mode 100644 make.gnu create mode 100644 mathunit.pas create mode 100644 romunit.pas diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..92b4f6f --- /dev/null +++ b/.gitignore @@ -0,0 +1,13 @@ +*.bmp +*.png +*.bak +*.ppu +*.o +*.zip +*.tar.gz +*.dat +Ergebnisse +lib +ROM +*~ +Log* diff --git a/Bild.pdf b/Bild.pdf new file mode 100644 index 0000000..52cdcb7 Binary files /dev/null and b/Bild.pdf differ diff --git a/ROM.lpi b/ROM.lpi new file mode 100644 index 0000000..d57d7be --- /dev/null +++ b/ROM.lpi @@ -0,0 +1,96 @@ + + + + + + + + + + + + + <UseAppBundle Value="False"/> + <ResourceType Value="res"/> + </General> + <i18n> + <EnableI18N LFM="False"/> + </i18n> + <VersionInfo> + <StringTable ProductVersion=""/> + </VersionInfo> + <BuildModes Count="1"> + <Item1 Name="Default" Default="True"/> + </BuildModes> + <PublishOptions> + <Version Value="2"/> + </PublishOptions> + <RunParams> + <local> + <FormatVersion Value="1"/> + </local> + </RunParams> + <Units Count="3"> + <Unit0> + <Filename Value="ROM.lpr"/> + <IsPartOfProject Value="True"/> + <UnitName Value="ROM"/> + </Unit0> + <Unit1> + <Filename Value="romunit.pas"/> + <IsPartOfProject Value="True"/> + <UnitName Value="romunit"/> + </Unit1> + <Unit2> + <Filename Value="mathunit.pas"/> + <IsPartOfProject Value="True"/> + <UnitName Value="mathunit"/> + </Unit2> + </Units> + </ProjectOptions> + <CompilerOptions> + <Version Value="11"/> + <Target> + <Filename Value="ROM"/> + </Target> + <SearchPaths> + <IncludeFiles Value="$(ProjOutDir)"/> + <UnitOutputDirectory Value="lib/$(TargetCPU)-$(TargetOS)"/> + </SearchPaths> + <CodeGeneration> + <Checks> + <IOChecks Value="True"/> + <RangeChecks Value="True"/> + <OverflowChecks Value="True"/> + <StackChecks Value="True"/> + </Checks> + <VerifyObjMethodCallValidity Value="True"/> + </CodeGeneration> + <Linking> + <Options> + <Win32> + <GraphicApplication Value="True"/> + </Win32> + </Options> + </Linking> + <Other> + <CompilerMessages> + <MsgFileName Value=""/> + </CompilerMessages> + <CompilerPath Value="$(CompPath)"/> + </Other> + </CompilerOptions> + <Debugging> + <Exceptions Count="3"> + <Item1> + <Name Value="EAbort"/> + </Item1> + <Item2> + <Name Value="ECodetoolError"/> + </Item2> + <Item3> + <Name Value="EFOpenError"/> + </Item3> + </Exceptions> + </Debugging> +</CONFIG> diff --git a/ROM.lpr b/ROM.lpr new file mode 100644 index 0000000..5863de0 --- /dev/null +++ b/ROM.lpr @@ -0,0 +1,150 @@ +program ROM; + +{$DEFINE UseCThreads} + +{$mode objfpc}{$H+} + +uses + {$IFDEF UNIX}{$IFDEF UseCThreads} + cthreads, + {$ENDIF}{$ENDIF} + Classes + { you can add units after this }, + SysUtils,ROMunit, mathunit, Math; + +var inPuls,refPuls,surTraj: tExtPointArray; + i,smooth: longint; + tmax,wmax,absShift: extended; + force,fourier: boolean; + +const Verwendung='Verwendung: ROM ($Einfallspuls_Datei $Ausfallspuls_Datei)/(- $trace-Datei-Prefix) $output_inPuls $output_refPuls $output_Trajektorie '+ + '[-s/--smooth $n] [-f/--force] [-t/--tmax $t] [-w/--wmax $w] [-F/--FFT] [-d/--dt $dt]'; + +begin + if paramcount<5 then Fehler(Verwendung); + i:=6; + force:=false; + smooth:=1; + tmax:=-1; + wmax:=-1; + fourier:=false; + absShift:=-1e9; + while i<=paramcount do begin + if (paramstr(i)='--force') or (paramstr(i)='-f') then begin + force:=true; + inc(i); + continue; + end; + if (paramstr(i)='--smooth') or (paramstr(i)='-s') then begin + inc(i); + smooth:=strtoint(paramstr(i)); + inc(i); + continue; + end; + if (paramstr(i)='--tmax') or (paramstr(i)='-t') then begin + inc(i); + tmax:=strtofloat(paramstr(i)); + inc(i); + continue; + end; + if (paramstr(i)='--wmax') or (paramstr(i)='-w') then begin + inc(i); + wmax:=strtofloat(paramstr(i)); + fourier:=true; + inc(i); + continue; + end; + if (paramstr(i)='--dt') or (paramstr(i)='-d') then begin + inc(i); + if paramstr(i)='auto' then begin + absShift:=-1e9; + inc(i); + continue; + end; + if paramstr(i)='input' then begin + if paramstr(1)<>'-' then Fehler('Ich brauche zur Bestimmung der Gesamtverschiebung die Inputdatei vom LPIC!'); + absShift:=-2e9; + inc(i); + continue; + end; + absShift:=strtofloat(paramstr(i)); + inc(i); + continue; + end; + if (paramstr(i)='--FFT') or (paramstr(i)='-F') then begin + fourier:=true; + inc(i); + continue; + end; + Fehler(Verwendung); + end; + if paramstr(1)='-' then + readRawInputs(paramstr(2),inPuls,refPuls,absShift) + else begin + readTextInput(paramstr(1),inPuls); + readTextInput(paramstr(2),refPuls); + end; + write(stderr,'Input sortieren ...'); + sort(inPuls); + sort(refPuls); + writeln(stderr,' fertig'); + uniq(inPuls,false); + uniq(refPuls,false); + write(stderr,'Input interpolieren ...'); + interpoliere(inPuls); + interpoliere(refPuls); + writeln(stderr,' fertig'); + flip(refPuls); + integrate(inPuls); + integrate(refPuls); + removeLinearOffset(inPuls); + removeLinearOffset(refPuls); + if smooth>1 then begin + write(stderr,'glätten ...'); + smoothen(inPuls,smooth); + smoothen(refPuls,smooth); + writeln(stderr,' fertig'); + end; + cut(inPuls,tmax); + cut(refPuls,tmax); + for i:=3 to 5 do + if (paramstr(i)<>'-') and fileexists(paramstr(i)) and not force then Fehler('Die Ausgabedatei '''+paramstr(i)+''' existiert bereits!'); + gesamtverschiebung(inPuls,refPuls,absShift); + write(stderr,'Trajektorie berechnen ...'); + berechneTrajektorie(inPuls,refPuls,surTraj,absShift*byte(not fourier)); + writeln(stderr,' fertig'); + write(stderr,'Ergebnis sortieren ...'); + sort(surTraj); + writeln(stderr,' fertig'); + uniq(surTraj,false); + if fourier then begin + write(stderr,'Ergebnis interpolieren ...'); + interpoliere(surTraj); + writeln(stderr,' fertig'); + fft(inPuls); + fft(refPuls); + fft(surTraj); + inPuls[0].y:=0; + refPuls[0].y:=0; + surTraj[0].y:=0; + if wmax<0 then begin + cut(surTraj,min(refPuls[length(refPuls)-1].x,inPuls[length(inPuls)-1].x)/2); + cut(inPuls,surTraj[length(surTraj)-1].x); + cut(refPuls,surTraj[length(surTraj)-1].x); + end + else begin + cut(surTraj,wmax); + cut(inPuls,wmax); + cut(refPuls,wmax); + end; + write(stderr,'alles normieren ...'); + normiere(inPuls); + normiere(refPuls); + normiere(surTraj); + writeln(stderr,' fertig'); + end; + if paramstr(3)<>'-' then writeOutput(paramstr(3),inPuls); + if paramstr(4)<>'-' then writeOutput(paramstr(4),refPuls); + if paramstr(5)<>'-' then writeOutput(paramstr(5),surTraj); +end. + diff --git a/ROM.lps b/ROM.lps new file mode 100644 index 0000000..b6e2d2d --- /dev/null +++ b/ROM.lps @@ -0,0 +1,166 @@ +<?xml version="1.0" encoding="UTF-8"?> +<CONFIG> + <ProjectSession> + <Version Value="9"/> + <BuildModes Active="Default"/> + <Units Count="3"> + <Unit0> + <Filename Value="ROM.lpr"/> + <IsPartOfProject Value="True"/> + <UnitName Value="ROM"/> + <IsVisibleTab Value="True"/> + <EditorIndex Value="0"/> + <WindowIndex Value="0"/> + <TopLine Value="1"/> + <CursorPos X="151" Y="20"/> + <UsageCount Value="90"/> + <Loaded Value="True"/> + </Unit0> + <Unit1> + <Filename Value="romunit.pas"/> + <IsPartOfProject Value="True"/> + <UnitName Value="romunit"/> + <EditorIndex Value="1"/> + <WindowIndex Value="0"/> + <TopLine Value="1004"/> + <CursorPos X="1" Y="998"/> + <FoldState Value=" T3iA041 pkRkZ0Y113A"/> + <UsageCount Value="90"/> + <Loaded Value="True"/> + </Unit1> + <Unit2> + <Filename Value="mathunit.pas"/> + <IsPartOfProject Value="True"/> + <UnitName Value="mathunit"/> + <EditorIndex Value="2"/> + <WindowIndex Value="0"/> + <TopLine Value="4"/> + <CursorPos X="1" Y="43"/> + <UsageCount Value="89"/> + <Loaded Value="True"/> + </Unit2> + </Units> + <General> + <ActiveWindowIndexAtStart Value="0"/> + </General> + <JumpHistory Count="29" HistoryIndex="28"> + <Position1> + <Filename Value="ROM.lpr"/> + <Caret Line="93" Column="58" TopLine="77"/> + </Position1> + <Position2> + <Filename Value="ROM.lpr"/> + <Caret Line="123" Column="66" TopLine="90"/> + </Position2> + <Position3> + <Filename Value="romunit.pas"/> + <Caret Line="577" Column="69" TopLine="393"/> + </Position3> + <Position4> + <Filename Value="romunit.pas"/> + <Caret Line="212" Column="2" TopLine="38"/> + </Position4> + <Position5> + <Filename Value="mathunit.pas"/> + <Caret Line="22" Column="1" TopLine="1"/> + </Position5> + <Position6> + <Filename Value="romunit.pas"/> + <Caret Line="528" Column="14" TopLine="476"/> + </Position6> + <Position7> + <Filename Value="romunit.pas"/> + <Caret Line="519" Column="37" TopLine="475"/> + </Position7> + <Position8> + <Filename Value="romunit.pas"/> + <Caret Line="542" Column="1" TopLine="512"/> + </Position8> + <Position9> + <Filename Value="romunit.pas"/> + <Caret Line="536" Column="34" TopLine="514"/> + </Position9> + <Position10> + <Filename Value="romunit.pas"/> + <Caret Line="537" Column="25" TopLine="508"/> + </Position10> + <Position11> + <Filename Value="romunit.pas"/> + <Caret Line="679" Column="28" TopLine="649"/> + </Position11> + <Position12> + <Filename Value="romunit.pas"/> + <Caret Line="333" Column="150" TopLine="305"/> + </Position12> + <Position13> + <Filename Value="romunit.pas"/> + <Caret Line="310" Column="3" TopLine="53"/> + </Position13> + <Position14> + <Filename Value="romunit.pas"/> + <Caret Line="26" Column="95" TopLine="1"/> + </Position14> + <Position15> + <Filename Value="romunit.pas"/> + <Caret Line="727" Column="6" TopLine="694"/> + </Position15> + <Position16> + <Filename Value="romunit.pas"/> + <Caret Line="26" Column="1" TopLine="1"/> + </Position16> + <Position17> + <Filename Value="romunit.pas"/> + <Caret Line="334" Column="31" TopLine="328"/> + </Position17> + <Position18> + <Filename Value="romunit.pas"/> + <Caret Line="311" Column="41" TopLine="288"/> + </Position18> + <Position19> + <Filename Value="romunit.pas"/> + <Caret Line="310" Column="3" TopLine="291"/> + </Position19> + <Position20> + <Filename Value="romunit.pas"/> + <Caret Line="22" Column="1" TopLine="16"/> + </Position20> + <Position21> + <Filename Value="romunit.pas"/> + <Caret Line="1018" Column="1" TopLine="891"/> + </Position21> + <Position22> + <Filename Value="romunit.pas"/> + <Caret Line="669" Column="16" TopLine="625"/> + </Position22> + <Position23> + <Filename Value="romunit.pas"/> + <Caret Line="621" Column="38" TopLine="527"/> + </Position23> + <Position24> + <Filename Value="romunit.pas"/> + <Caret Line="312" Column="1" TopLine="295"/> + </Position24> + <Position25> + <Filename Value="romunit.pas"/> + <Caret Line="657" Column="15" TopLine="626"/> + </Position25> + <Position26> + <Filename Value="romunit.pas"/> + <Caret Line="655" Column="16" TopLine="637"/> + </Position26> + <Position27> + <Filename Value="romunit.pas"/> + <Caret Line="285" Column="9" TopLine="254"/> + </Position27> + <Position28> + <Filename Value="romunit.pas"/> + <Caret Line="312" Column="49" TopLine="301"/> + </Position28> + <Position29> + <Filename Value="romunit.pas"/> + <Caret Line="21" Column="69" TopLine="1"/> + </Position29> + </JumpHistory> + </ProjectSession> + <EditorMacros Count="0"/> +</CONFIG> diff --git a/make.gnu b/make.gnu new file mode 100644 index 0000000..9a85259 --- /dev/null +++ b/make.gnu @@ -0,0 +1,26 @@ +set terminal pdf +set decimalsign "," + +set output "~/Dokumente/Prograemmchen/ROM/Bild.pdf" +set title "Oberflächentrajektorie" +set size .9,.9 +set format x "%g" +set format y "%g" + +if (system("echo $logscale")=="1") { + set xlabel "omega in omega_L" + set xrange [0:40] + set yrange [10**-11:*] + set logscale y +} +else { + set xlabel "t in T" +} + +# set ylabel sprintf("$\\left|\\Big(\\mathcal{FT}\\left(%s\\right)\\Big)(\\omega)\\right|$ in $%s$",groesze,einheit) +# set xtics out +# set ytics 10**-2 out +plot "~/Dokumente/Prograemmchen/ROM/surface.dat" using 1:2 with lines lt 1 title "surface", \ + "~/Dokumente/Prograemmchen/ROM/input.dat" using 1:2 with lines lt 2 title "input", \ + "~/Dokumente/Prograemmchen/ROM/output.dat" using 1:2 with lines lt 3 title "output" + diff --git a/mathunit.pas b/mathunit.pas new file mode 100644 index 0000000..3d5e111 --- /dev/null +++ b/mathunit.pas @@ -0,0 +1,46 @@ +unit mathunit; + +{$mode objfpc}{$H+} + +interface + +uses + Classes, SysUtils; + +type + tExtPoint = record + x,y: extended; + end; + tExtPointArray = array of tExtPoint; + pTExtPointArray = ^tExtPointArray; + tLongintArray = array of longint; + pTLongintArray = ^tLongintArray; + tExtendedArray = array of extended; + +function Plus(a,b: tExtPoint): tExtPoint; +function Durch(a: tExtPoint; b: extended): tExtPoint; +function myFrac(x: extended): extended; + +implementation + +function Plus(a,b: tExtPoint): tExtPoint; +begin + result.x:=a.x+b.x; + result.y:=a.y+b.y; +end; + +function Durch(a: tExtPoint; b: extended): tExtPoint; +begin + result.x:=a.x/b; + result.y:=a.y/b; +end; + +function myFrac(x: extended): extended; +begin + result:=frac(x); + while result<0 do + result:=result+1; +end; + +end. + 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. + -- cgit v1.2.3-70-g09d2