From 058b8a39480c87a3b31876dce87bc383b506661a Mon Sep 17 00:00:00 2001 From: Erich Eckner Date: Tue, 31 Jul 2018 14:59:24 +0200 Subject: überflüssiges entfernt - wir können jetzt "nur noch" x_arp rekonstruieren (keine fft, kein fenstern, ...) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- ROM.lpr | 138 ++++++++++------------------------------- ROM.lps | 102 +++++++++++++++--------------- romunit.pas | 202 +----------------------------------------------------------- 3 files changed, 84 insertions(+), 358 deletions(-) diff --git a/ROM.lpr b/ROM.lpr index 0dc1eba..6045da6 100644 --- a/ROM.lpr +++ b/ROM.lpr @@ -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. diff --git a/ROM.lps b/ROM.lps index 222e634..6d5de5b 100644 --- a/ROM.lps +++ b/ROM.lps @@ -7,9 +7,9 @@ - - - + + + @@ -17,10 +17,10 @@ - - - - + + + + @@ -28,7 +28,7 @@ - + @@ -59,123 +59,121 @@ - + - + - - + - - + + - + - - + + - - + + - - + + - - + - - + + - - + + - - + + - - + + - + - + - - + + - + - - + + - + - - + + - + - + - + - + - - + + - - + + - + - + - + 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=0) and - (pps[j]= 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. -- cgit v1.2.3-54-g00ecf