diff options
| -rw-r--r-- | .gitignore | 13 | ||||
| -rw-r--r-- | Bild.pdf | bin | 0 -> 74662 bytes | |||
| -rw-r--r-- | ROM.lpi | 96 | ||||
| -rw-r--r-- | ROM.lpr | 150 | ||||
| -rw-r--r-- | ROM.lps | 166 | ||||
| -rw-r--r-- | make.gnu | 26 | ||||
| -rw-r--r-- | mathunit.pas | 46 | ||||
| -rw-r--r-- | romunit.pas | 1086 |
8 files changed, 1583 insertions, 0 deletions
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 Binary files differnew file mode 100644 index 0000000..52cdcb7 --- /dev/null +++ b/Bild.pdf @@ -0,0 +1,96 @@ +<?xml version="1.0" encoding="UTF-8"?> +<CONFIG> + <ProjectOptions> + <Version Value="9"/> + <General> + <Flags> + <MainUnitHasCreateFormStatements Value="False"/> + <MainUnitHasTitleStatement Value="False"/> + </Flags> + <SessionStorage Value="InProjectDir"/> + <MainUnit Value="0"/> + <Title Value="ROM"/> + <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> @@ -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. + @@ -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. + |