diff options
Diffstat (limited to 'romunit.pas')
-rw-r--r-- | romunit.pas | 202 |
1 files changed, 1 insertions, 201 deletions
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. |