diff options
-rw-r--r-- | ROM.lpr | 138 | ||||
-rw-r--r-- | ROM.lps | 102 | ||||
-rw-r--r-- | romunit.pas | 202 |
3 files changed, 84 insertions, 358 deletions
@@ -14,10 +14,10 @@ uses var inPulsO,inPuls,refPulsO,refPuls,surTraj,cRefPuls, - surVel,inPulsArg,refPulsArg,surTrajArg,surVelArg: tExtPointArray; + surVel: tExtPointArray; smooth,betaSmooth,veloSmooth: longint; - tmax,wmax,absShift,betaBound,veloBound,fftBreite: extended; - force,fourier,mitAmplMod: boolean; + tmax,absShift,betaBound,veloBound: extended; + force,mitAmplMod: boolean; f,bekannteBefehle: tMyStringList; s,lpicIn,rohIn,rohRef,outIn, outRef,outRefC,outSur,outVel: string; @@ -34,11 +34,8 @@ begin smooth:=1; betaSmooth:=1; tmax:=-1; - wmax:=-1; absShift:=-1e9; betaBound:=0.95; - fourier:=false; - fftBreite:=-1; mitAmplMod:=true; veloSmooth:=1; veloBound:=1; @@ -99,11 +96,6 @@ begin tmax:=strtofloat(s); continue; end; - if istDasBefehl('wmax:',s,bekannteBefehle,true) or - istDasBefehl('ωmax:',s,bekannteBefehle,true) then begin - wmax:=strtofloat(s); - continue; - end; if istDasBefehl('Absolutverschiebung:',s,bekannteBefehle,true) then begin if s='auto' then begin absShift:=-1e9; @@ -116,18 +108,6 @@ begin absShift:=strtofloat(s); continue; end; - if s='mit FFT' then begin - fourier:=true; - continue; - end; - if istDasBefehl('FFT-Breite:',s,bekannteBefehle,true) then begin - fftBreite:=strToFloat(s); - continue; - end; - if s='ohne FFT' then begin - fourier:=false; - continue; - end; if istDasBefehl('lpic-Quelle:',s,bekannteBefehle,true) then begin lpicIn:=s; continue; @@ -224,97 +204,45 @@ begin end; cut(inPuls,tmax); cut(refPuls,tmax); - 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); - write(stderr,'Geschwindigkeit berechnen ...'); - ableiten(surTraj,surVel,veloSmooth,veloBound); - writeln(stderr,' fertig'); - write(stderr,'Reflektierten Puls berechnen ...'); - berechneRefPuls(inPulsO,surTraj,betaSmooth,betaBound,cRefPuls); - writeln(stderr,' fertig'); - if fourier then begin - write(stderr,'Ergebnis interpolieren ...'); - interpoliere(surTraj); - if outVel<>'' then - interpoliere(surVel); + if (outVel<>'') or (outRefC<>'') or (outSur<>'') then begin + gesamtverschiebung(inPuls,refPuls,absShift); + write(stderr,'Trajektorie berechnen ...'); + berechneTrajektorie(inPuls,refPuls,surTraj,absShift); writeln(stderr,' fertig'); - if fftBreite>0 then begin - fenstern(inPuls,fftBreite); - fenstern(refPuls,fftBreite); - fenstern(surTraj,fftBreite); - if outVel<>'' then - fenstern(surVel,fftBreite,true); - end; - fft(inPuls,inPulsArg); - fft(refPuls,refPulsArg); - fft(surTraj,surTrajArg); - if outVel<>'' then - fft(surVel,surVelArg); - 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(surTrajArg,min(refPuls[length(refPuls)-1]['x'],inPuls[length(inPuls)-1]['x'])/2); - cut(inPuls,surTraj[length(surTraj)-1]['x']); - cut(inPulsArg,surTraj[length(surTraj)-1]['x']); - cut(refPuls,surTraj[length(surTraj)-1]['x']); - cut(refPulsArg,surTraj[length(surTraj)-1]['x']); - if outVel<>'' then begin - cut(surVel,surTraj[length(surTraj)-1]['x']); - cut(surVelArg,surTraj[length(surTraj)-1]['x']); - end; - end - else begin - cut(surTraj,wmax); - cut(surTrajArg,wmax); - cut(inPuls,wmax); - cut(inPulsArg,wmax); - cut(refPuls,wmax); - cut(refPulsArg,wmax); - if outVel<>'' then begin - cut(surVel,wmax); - cut(surVelArg,wmax); - end; - end; - write(stderr,'alles normieren ...'); - normiere(inPuls); - normiere(refPuls); - normiere(surTraj); + write(stderr,'Ergebnis sortieren ...'); + sort(surTraj); writeln(stderr,' fertig'); + uniq(surTraj,false); end; + if (outVel<>'') or (outRefC<>'') then begin + write(stderr,'Geschwindigkeit berechnen ...'); + ableiten(surTraj,surVel,veloSmooth,veloBound); + writeln(stderr,' fertig'); + end; + if outRefC<>'' then begin + write(stderr,'Reflektierten Puls berechnen ...'); + berechneRefPuls(inPulsO,surTraj,betaSmooth,betaBound,cRefPuls); + writeln(stderr,' fertig'); + end; + write(stderr,'Ergebnis interpolieren ...'); + if outSur<>'' then + interpoliere(surTraj); + if outVel<>'' then + interpoliere(surVel); + writeln(stderr,' fertig'); + if outIn<>'' then begin writeOutput(outIn+'.ori',inPulsO); - if fourier then - writeOutput(outIn,inPuls,inPulsArg) - else - writeOutput(outIn,inPuls); + writeOutput(outIn,inPuls); end; if outRef<>'' then begin writeOutput(outRef+'.ori',refPulsO); - if fourier then - writeOutput(outRef,refPuls,refPulsArg) - else - writeOutput(outRef,refPuls); - end; - if outSur<>'' then begin - if fourier then - writeOutput(outSur,surTraj,surTrajArg) - else - writeOutput(outSur,surTraj); - end; - if outVel<>'' then begin - if fourier then - writeOutput(outVel,surVel,surVelArg) - else - writeOutput(outVel,surVel); + writeOutput(outRef,refPuls); end; + if outSur<>'' then + writeOutput(outSur,surTraj); + if outVel<>'' then + writeOutput(outVel,surVel); if outRefC<>'' then writeOutput(outRefC,cRefPuls); end. @@ -7,9 +7,9 @@ <Unit0> <Filename Value="ROM.lpr"/> <IsPartOfProject Value="True"/> - <TopLine Value="91"/> - <CursorPos X="5" Y="31"/> - <UsageCount Value="108"/> + <TopLine Value="177"/> + <CursorPos X="19" Y="206"/> + <UsageCount Value="110"/> <Loaded Value="True"/> </Unit0> <Unit1> @@ -17,10 +17,10 @@ <IsPartOfProject Value="True"/> <IsVisibleTab Value="True"/> <EditorIndex Value="1"/> - <TopLine Value="891"/> - <CursorPos X="59" Y="900"/> - <FoldState Value=" T3iF04053 pidrN0A1="/> - <UsageCount Value="108"/> + <TopLine Value="286"/> + <CursorPos Y="306"/> + <FoldState Value=" T3iC04043 pidqc0A1I"/> + <UsageCount Value="110"/> <Loaded Value="True"/> </Unit1> <Unit2> @@ -28,7 +28,7 @@ <IsPartOfProject Value="True"/> <EditorIndex Value="-1"/> <CursorPos Y="10"/> - <UsageCount Value="107"/> + <UsageCount Value="109"/> </Unit2> <Unit3> <Filename Value="../units/matheunit.pas"/> @@ -59,123 +59,121 @@ <JumpHistory Count="30" HistoryIndex="29"> <Position1> <Filename Value="ROM.lpr"/> - <Caret Line="26" Column="19"/> + <Caret Line="320" Column="33" TopLine="298"/> </Position1> <Position2> <Filename Value="ROM.lpr"/> - <Caret Line="33" Column="9" TopLine="4"/> + <Caret Line="322" Column="33" TopLine="298"/> </Position2> <Position3> <Filename Value="ROM.lpr"/> - <Caret Line="93" Column="13" TopLine="64"/> </Position3> <Position4> <Filename Value="ROM.lpr"/> - <Caret Line="256" Column="21" TopLine="234"/> + <Caret Line="16" Column="42"/> </Position4> <Position5> - <Filename Value="romunit.pas"/> - <Caret Line="23" Column="29" TopLine="5"/> + <Filename Value="ROM.lpr"/> + <Caret Line="243" Column="9" TopLine="242"/> </Position5> <Position6> <Filename Value="ROM.lpr"/> - <Caret Line="258" Column="14" TopLine="234"/> + <Caret Line="16" Column="42"/> </Position6> <Position7> - <Filename Value="romunit.pas"/> - <Caret Line="15" Column="30"/> + <Filename Value="ROM.lpr"/> + <Caret Line="238" Column="45" TopLine="209"/> </Position7> <Position8> - <Filename Value="romunit.pas"/> - <Caret Line="32" Column="23" TopLine="14"/> + <Filename Value="ROM.lpr"/> + <Caret Line="241" Column="15" TopLine="212"/> </Position8> <Position9> - <Filename Value="romunit.pas"/> - <Caret Line="876" Column="31"/> + <Filename Value="ROM.lpr"/> + <Caret Line="225" Column="3" TopLine="204"/> </Position9> <Position10> - <Filename Value="romunit.pas"/> - <Caret Line="888" TopLine="865"/> + <Filename Value="ROM.lpr"/> </Position10> <Position11> - <Filename Value="romunit.pas"/> - <Caret Line="36"/> + <Filename Value="ROM.lpr"/> + <Caret Line="284" Column="42" TopLine="257"/> </Position11> <Position12> - <Filename Value="romunit.pas"/> - <Caret Line="1063" Column="40" TopLine="1041"/> + <Filename Value="ROM.lpr"/> + <Caret Line="19" Column="58" TopLine="3"/> </Position12> <Position13> - <Filename Value="romunit.pas"/> - <Caret Line="1062" Column="26" TopLine="1041"/> + <Filename Value="ROM.lpr"/> + <Caret Line="40" TopLine="22"/> </Position13> <Position14> - <Filename Value="romunit.pas"/> - <Caret Line="1068" Column="72" TopLine="1042"/> + <Filename Value="ROM.lpr"/> + <Caret Line="117" TopLine="100"/> </Position14> <Position15> <Filename Value="ROM.lpr"/> - <Caret Line="286" Column="22" TopLine="139"/> + <Caret Line="216" Column="56" TopLine="198"/> </Position15> <Position16> <Filename Value="ROM.lpr"/> - <Caret Line="280" Column="25" TopLine="254"/> + <Caret Line="250" Column="22" TopLine="219"/> </Position16> <Position17> - <Filename Value="ROM.lpr"/> - <Caret Line="17" Column="11"/> + <Filename Value="romunit.pas"/> + <Caret Line="30" TopLine="28"/> </Position17> <Position18> <Filename Value="ROM.lpr"/> - <Caret Line="271" Column="26" TopLine="257"/> + <Caret Line="19" Column="58"/> </Position18> <Position19> - <Filename Value="romunit.pas"/> - <Caret Line="26" Column="19" TopLine="8"/> + <Filename Value="ROM.lpr"/> + <Caret Line="37" TopLine="19"/> </Position19> <Position20> <Filename Value="ROM.lpr"/> - <Caret Line="286" Column="37" TopLine="257"/> + <Caret Line="99" TopLine="83"/> </Position20> <Position21> - <Filename Value="romunit.pas"/> - <Caret Line="35" Column="61" TopLine="17"/> + <Filename Value="ROM.lpr"/> + <Caret Line="19" Column="7"/> </Position21> <Position22> <Filename Value="ROM.lpr"/> - <Caret Line="151" Column="15" TopLine="122"/> + <Caret Line="36" Column="7" TopLine="7"/> </Position22> <Position23> <Filename Value="ROM.lpr"/> - <Caret Line="21" Column="58" TopLine="16"/> + <Caret Line="95" Column="26" TopLine="66"/> </Position23> <Position24> <Filename Value="ROM.lpr"/> - <Caret Line="183" Column="25" TopLine="154"/> + <Caret Line="96" Column="11" TopLine="67"/> </Position24> <Position25> <Filename Value="ROM.lpr"/> - <Caret Line="165" Column="6" TopLine="152"/> + <Caret Line="205" Column="18" TopLine="176"/> </Position25> <Position26> - <Filename Value="ROM.lpr"/> - <Caret Line="59" Column="29" TopLine="35"/> + <Filename Value="romunit.pas"/> + <Caret Line="30"/> </Position26> <Position27> - <Filename Value="ROM.lpr"/> - <Caret Line="166" Column="9" TopLine="147"/> + <Filename Value="romunit.pas"/> + <Caret Line="888" Column="63" TopLine="871"/> </Position27> <Position28> <Filename Value="romunit.pas"/> - <Caret Line="1036" Column="116" TopLine="1015"/> + <Caret Line="1009" TopLine="975"/> </Position28> <Position29> <Filename Value="romunit.pas"/> - <Caret Line="337" Column="62" TopLine="319"/> + <Caret Line="13" Column="64"/> </Position29> <Position30> <Filename Value="romunit.pas"/> - <Caret Line="31" Column="14" TopLine="13"/> + <Caret Line="272" Column="22" TopLine="243"/> </Position30> </JumpHistory> </ProjectSession> diff --git a/romunit.pas b/romunit.pas index c1c0a7c..761bd68 100644 --- a/romunit.pas +++ b/romunit.pas @@ -10,8 +10,7 @@ uses 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); overload; -procedure writeOutput(nam: string; const dat,datArg: tExtPointArray); overload; +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; @@ -28,11 +27,9 @@ 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; out datArg: tExtPointArray); procedure interpoliere(var dat: tExtPointArray); procedure normiere(var dat: tExtPointArray); procedure berechneRefPuls(inPuls,surTraj: tExtPointArray; betaGlaette: longint; betaBound: extended; out cRefPuls: tExtPointArray); -procedure fenstern(var dat: tExtPointArray; breite: extended; alteMitte: boolean = false); type tSortThread = class(tThread) @@ -306,34 +303,6 @@ begin writeln(stderr,' fertig'); end; -procedure writeOutput(nam: string; const dat,datArg: tExtPointArray); overload; -var - f: textfile; - i: longint; -begin - if length(dat)=0 then - exit; - if length(dat)<>length(datArg) then - Fehler('dat und datArg haben unterschiedlich viele Werte!'); - write(stderr,'Datei '''+nam+''' schreiben '); - assignfile(f,nam); - rewrite(f); - for i:=0 to length(dat)-1 do begin - if dat[i]['x']<>datArg[i]['x'] then - Fehler('dat und datArg haben unterschiedliche x-Werte an Position '+inttostr(i)+' ('+floattostr(dat[i]['x'])+' vs. '+floattostr(datArg[i]['x'])+')!'); - writeln( - f, - floattostr(dat[i]['x'])+#9+ - floattostr(dat[i]['y'])+#9+ - floattostr(datArg[i]['y']) - ); - 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; @@ -811,119 +780,6 @@ begin setlength(dat,j); end; -procedure fft(var dat: tExtPointArray; out datArg: 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 ('+inttostr(length(dat))+') ... '); - 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[2*k-1-i]['y'] * sqr(sin(pi/2*i/k)); - for i:=j+k to length(dat)-1 do // dito - dat[i]['y']:=dat[2*j+2*k-1-i]['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; - - setlength(datArg,length(dat)); - - t1:=0; - for i:=0 to 2*haL-1 do begin - t2:=arctan2(dat[i]['y'],dat[i]['x']); - t2:=t2-2*pi*round((t2-t1)/2/pi); - datArg[i]['y']:=t2; - t1:=t2; - dat[i]['y']:=(sqr(dat[i]['x'])+sqr(dat[i]['y']))/(2*haL); - dat[i]['x']:=fStep*i; - datArg[i]['x']:=dat[i]['x']; - end; - for i:=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; @@ -1121,61 +977,5 @@ begin fertig:=true; end; -var - fensterMitte: longint; - -procedure fenstern(var dat: tExtPointArray; breite: extended; alteMitte: boolean = false); -var - i,j,len: longint; - extrema: tLongintArray; - pps: tExtendedArray; - mpp: extended; -begin - if not alteMitte then begin - setlength(extrema,0); - len:=0; - for i:=1 to length(dat)-2 do - if (dat[i]['y']>dat[i-1]['y']) xor - (dat[i+1]['y']>dat[i]['y']) then begin - if length(extrema)<=len then - setlength(extrema,len+65536); - extrema[len]:=i; - inc(len); - end; - setlength(extrema,len); - setlength(pps,len); - for i:=0 to length(extrema)-1 do begin - pps[i]:=-1; - for j:=0 to length(extrema)-1 do - if (abs(dat[extrema[j]]['x']-dat[extrema[i]]['x'])<2) and - (abs(dat[extrema[j]]['y']-dat[extrema[i]]['y'])>pps[i]) then - pps[i]:=abs(dat[extrema[j]]['y']-dat[extrema[i]]['y']); - end; - mpp:=pps[0]; - for i:=0 to length(pps)-1 do - if pps[i]>mpp then - mpp:=pps[i]; - i:=0; - while (i<length(pps)) and - (pps[i]<mpp/2) do - inc(i); - j:=length(pps)-1; - while (j>=0) and - (pps[j]<mpp/2) do - dec(j); - mpp:=(dat[extrema[i]]['x']+dat[extrema[j]]['x'])/2; - fensterMitte:=0; - while (fensterMitte<length(dat)) and - (dat[fensterMitte]['x']<mpp) do - inc(fensterMitte); - end; - - for i:=0 to length(dat)-1 do - if abs(dat[i]['x']-dat[fensterMitte]['x']) >= breite then - dat[i]['y']:=0 - else - dat[i]['y']:=dat[i]['y'] * (1+cos((dat[i]['x']-dat[fensterMitte]['x'])/breite*pi))/2; -end; - end. |