From 59cecd0385f90e25ce6c73782d5b01bba5bc5f7c Mon Sep 17 00:00:00 2001 From: Erich Eckner Date: Mon, 30 Jul 2018 15:26:47 +0200 Subject: fft nun fensterbar --- ROM.lpr | 30 ++++++++++++--- ROM.lps | 122 +++++++++++++++++++++++++++++++----------------------------- romunit.pas | 44 ++++++++++++++++++++++ 3 files changed, 131 insertions(+), 65 deletions(-) diff --git a/ROM.lpr b/ROM.lpr index 30901d5..034fb35 100644 --- a/ROM.lpr +++ b/ROM.lpr @@ -16,7 +16,7 @@ var inPulsO,inPuls,refPulsO,refPuls,surTraj,cRefPuls, surVel,inPulsArg,refPulsArg,surTrajArg,surVelArg: tExtPointArray; smooth,betaSmooth,veloSmooth: longint; - tmax,wmax,absShift,betaBound,veloBound: extended; + tmax,wmax,absShift,betaBound,veloBound,fftBreite: extended; force,fourier,mitAmplMod: boolean; f: textfile; s,t,u,lpicIn,rohIn,rohRef,outIn, @@ -37,6 +37,7 @@ begin absShift:=-1e9; betaBound:=0.95; fourier:=false; + fftBreite:=-1; mitAmplMod:=true; veloSmooth:=1; veloBound:=1; @@ -147,6 +148,11 @@ begin fourier:=true; continue; end; + if pos('FFT-Breite:',s)=1 then begin + delete(s,1,pos(':',s)); + fftBreite:=strToFloat(trim(s)); + continue; + end; if s='ohne FFT' then begin fourier:=false; continue; @@ -273,10 +279,18 @@ begin if outVel<>'' then interpoliere(surVel); 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); - fft(surVel,surVelArg); + if outVel<>'' then + fft(surVel,surVelArg); inPuls[0]['y']:=0; refPuls[0]['y']:=0; surTraj[0]['y']:=0; @@ -287,8 +301,10 @@ begin 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']); + 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); @@ -297,8 +313,10 @@ begin cut(inPulsArg,wmax); cut(refPuls,wmax); cut(refPulsArg,wmax); - cut(surVel,wmax); - cut(surVelArg,wmax); + if outVel<>'' then begin + cut(surVel,wmax); + cut(surVelArg,wmax); + end; end; write(stderr,'alles normieren ...'); normiere(inPuls); diff --git a/ROM.lps b/ROM.lps index 8c82481..90b1a22 100644 --- a/ROM.lps +++ b/ROM.lps @@ -7,9 +7,9 @@ - - - + + + @@ -17,10 +17,10 @@ - - - - + + + + @@ -28,7 +28,7 @@ - + @@ -40,131 +40,135 @@ - + - + - + - + - - + + - - + - - + + - - + + - - + - - + + - - + + - - + + - - + + - - + + - + - + + - + - - + + - - + + - + + - - + + - - + + - - + + - - + + - - + + - - + + - + - - + + - + - + - - + + + + + + diff --git a/romunit.pas b/romunit.pas index 41fa0c3..28e0df5 100644 --- a/romunit.pas +++ b/romunit.pas @@ -32,6 +32,7 @@ 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) @@ -1055,5 +1056,48 @@ begin fertig:=true; end; +var + fensterMitte: longint; + +procedure fenstern(var dat: tExtPointArray; breite: extended; alteMitte: boolean = false); +var + i,j,len: longint; + extrema: tLongintArray; + pp,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); + fensterMitte:=-1; + mpp:=-1; + for i:=0 to length(extrema)-1 do begin + pp:=-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'])>pp) then + pp:=abs(dat[extrema[j]]['y']-dat[extrema[i]]['y']); + if pp>mpp then begin + mpp:=pp; + fensterMitte:=extrema[i]; + end; + end; + 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. -- cgit v1.2.3-70-g09d2