From 807f9bc2be261d6ea1ac41777f93d0895a81f444 Mon Sep 17 00:00:00 2001 From: Erich Eckner Date: Thu, 12 Jul 2018 11:24:50 +0200 Subject: bei fft Winkel mit ausgeben --- ROM.lpr | 45 ++++++++++++++++++------ ROM.lps | 115 ++++++++++++++++++++++++++++++------------------------------ romunit.pas | 42 +++++++++++++++++++--- 3 files changed, 130 insertions(+), 72 deletions(-) diff --git a/ROM.lpr b/ROM.lpr index 941a700..30901d5 100644 --- a/ROM.lpr +++ b/ROM.lpr @@ -13,8 +13,8 @@ uses SysUtils,ROMunit, matheunit, Math, systemunit, lowlevelunit; var - inPulsO,inPuls,refPulsO,refPuls, - surTraj,cRefPuls,surVel: tExtPointArray; + inPulsO,inPuls,refPulsO,refPuls,surTraj,cRefPuls, + surVel,inPulsArg,refPulsArg,surTrajArg,surVelArg: tExtPointArray; smooth,betaSmooth,veloSmooth: longint; tmax,wmax,absShift,betaBound,veloBound: extended; force,fourier,mitAmplMod: boolean; @@ -273,23 +273,32 @@ begin if outVel<>'' then interpoliere(surVel); writeln(stderr,' fertig'); - fft(inPuls); - fft(refPuls); - fft(surTraj); - fft(surVel); + fft(inPuls,inPulsArg); + fft(refPuls,refPulsArg); + fft(surTraj,surTrajArg); + 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']); cut(surVel,surTraj[length(surTraj)-1]['x']); + cut(surVelArg,surTraj[length(surTraj)-1]['x']); end else begin cut(surTraj,wmax); + cut(surTrajArg,wmax); cut(inPuls,wmax); + cut(inPulsArg,wmax); cut(refPuls,wmax); + cut(refPulsArg,wmax); + cut(surVel,wmax); + cut(surVelArg,wmax); end; write(stderr,'alles normieren ...'); normiere(inPuls); @@ -299,14 +308,30 @@ begin end; if outIn<>'' then begin writeOutput(outIn+'.ori',inPulsO); - writeOutput(outIn,inPuls); + if fourier then + writeOutput(outIn,inPuls,inPulsArg) + else + writeOutput(outIn,inPuls); end; if outRef<>'' then begin writeOutput(outRef+'.ori',refPulsO); - writeOutput(outRef,refPuls); + 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); 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 8e5b84a..1381843 100644 --- a/ROM.lps +++ b/ROM.lps @@ -7,9 +7,9 @@ - - - + + + @@ -17,10 +17,10 @@ - - - - + + + + @@ -28,7 +28,7 @@ - + @@ -52,123 +52,122 @@ - + - - + + - - + + - - + + - + - - + + - - + + - - + + - + - - + + - + - - + + - + - - + + - - + + - - + + - - + + - - + + - - + + - + - + - - + + - - + + - - + + - - + + - - - + + - + - - + + - + diff --git a/romunit.pas b/romunit.pas index 6bce2e0..ff8ea7c 100644 --- a/romunit.pas +++ b/romunit.pas @@ -10,7 +10,8 @@ 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); +procedure writeOutput(nam: string; const dat: tExtPointArray); overload; +procedure writeOutput(nam: string; const dat,datArg: tExtPointArray); overload; procedure berechneTrajektorie(const inPuls,outPuls: tExtPointArray; out surTraj: tExtPointArray; absShift: extended); procedure smoothen(var dat: tExtPointArray; width: longint); procedure sort(var dat: tExtPointArray); overload; @@ -27,7 +28,7 @@ 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 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); @@ -284,6 +285,31 @@ 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; exList: array[0..3] of tLongintArray; @@ -734,14 +760,14 @@ begin setlength(dat,j); end; -procedure fft(var dat: tExtPointArray); +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 ... '); + 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))))); @@ -815,9 +841,17 @@ begin 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); -- cgit v1.2.3-70-g09d2