bei fft Winkel mit ausgeben

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 @@
       <Unit0>
         <Filename Value="ROM.lpr"/>
         <IsPartOfProject Value="True"/>
-        <TopLine Value="248"/>
-        <CursorPos X="26" Y="265"/>
-        <UsageCount Value="98"/>
+        <TopLine Value="18"/>
+        <CursorPos X="28" Y="29"/>
+        <UsageCount Value="99"/>
         <Loaded Value="True"/>
       </Unit0>
       <Unit1>
@@ -17,10 +17,10 @@
         <IsPartOfProject Value="True"/>
         <IsVisibleTab Value="True"/>
         <EditorIndex Value="1"/>
-        <TopLine Value="618"/>
-        <CursorPos X="24" Y="640"/>
-        <FoldState Value=" T3iD04042H p0-a0747"/>
-        <UsageCount Value="98"/>
+        <TopLine Value="834"/>
+        <CursorPos X="14" Y="852"/>
+        <FoldState Value=" T3iE040528"/>
+        <UsageCount Value="99"/>
         <Loaded Value="True"/>
       </Unit1>
       <Unit2>
@@ -28,7 +28,7 @@
         <IsPartOfProject Value="True"/>
         <EditorIndex Value="-1"/>
         <CursorPos Y="10"/>
-        <UsageCount Value="97"/>
+        <UsageCount Value="98"/>
       </Unit2>
       <Unit3>
         <Filename Value="../units/matheunit.pas"/>
@@ -52,123 +52,122 @@
     <JumpHistory Count="30" HistoryIndex="29">
       <Position1>
         <Filename Value="romunit.pas"/>
-        <Caret Line="856" Column="47" TopLine="838"/>
+        <Caret Line="1011" Column="58" TopLine="984"/>
       </Position1>
       <Position2>
-        <Filename Value="romunit.pas"/>
-        <Caret Line="858" Column="14" TopLine="841"/>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="13" Column="62"/>
       </Position2>
       <Position3>
-        <Filename Value="romunit.pas"/>
-        <Caret Line="867" Column="12" TopLine="850"/>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="282" Column="20" TopLine="262"/>
       </Position3>
       <Position4>
-        <Filename Value="romunit.pas"/>
-        <Caret Line="869" Column="25" TopLine="851"/>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="279" Column="16" TopLine="266"/>
       </Position4>
       <Position5>
         <Filename Value="romunit.pas"/>
-        <Caret Line="872" Column="29" TopLine="854"/>
+        <Caret Line="30" Column="14" TopLine="23"/>
       </Position5>
       <Position6>
-        <Filename Value="romunit.pas"/>
-        <Caret Line="885" Column="57" TopLine="867"/>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="279" Column="16" TopLine="266"/>
       </Position6>
       <Position7>
-        <Filename Value="romunit.pas"/>
-        <Caret Line="887" Column="96" TopLine="869"/>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="17" Column="28"/>
       </Position7>
       <Position8>
-        <Filename Value="romunit.pas"/>
-        <Caret Line="893" Column="119" TopLine="875"/>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="265" Column="26" TopLine="248"/>
       </Position8>
       <Position9>
         <Filename Value="romunit.pas"/>
-        <Caret Line="895" Column="106" TopLine="883"/>
+        <Caret Line="25" Column="19" TopLine="7"/>
       </Position9>
       <Position10>
-        <Filename Value="romunit.pas"/>
-        <Caret Line="899" Column="26" TopLine="879"/>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="295" Column="9" TopLine="274"/>
       </Position10>
       <Position11>
         <Filename Value="romunit.pas"/>
-        <Caret Line="966" Column="33" TopLine="930"/>
+        <Caret Line="32" Column="19" TopLine="14"/>
       </Position11>
       <Position12>
-        <Filename Value="romunit.pas"/>
-        <Caret Line="969" Column="40" TopLine="950"/>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="295" Column="9" TopLine="268"/>
       </Position12>
       <Position13>
         <Filename Value="romunit.pas"/>
-        <Caret Line="976" Column="47" TopLine="958"/>
+        <Caret Line="30" Column="14" TopLine="12"/>
       </Position13>
       <Position14>
-        <Filename Value="romunit.pas"/>
-        <Caret Line="981" Column="33" TopLine="961"/>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="335" Column="46" TopLine="302"/>
       </Position14>
       <Position15>
-        <Filename Value="romunit.pas"/>
-        <Caret Line="987" Column="71" TopLine="970"/>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="17" Column="21"/>
       </Position15>
       <Position16>
-        <Filename Value="romunit.pas"/>
-        <Caret Line="989" Column="32" TopLine="972"/>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="268" Column="64" TopLine="239"/>
       </Position16>
       <Position17>
-        <Filename Value="romunit.pas"/>
-        <Caret Line="991" Column="141" TopLine="974"/>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="334" Column="51" TopLine="302"/>
       </Position17>
       <Position18>
-        <Filename Value="romunit.pas"/>
-        <Caret Line="992" Column="108" TopLine="975"/>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="17" Column="58"/>
       </Position18>
       <Position19>
-        <Filename Value="romunit.pas"/>
-        <Caret Line="1000" Column="33" TopLine="983"/>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="276" Column="5" TopLine="258"/>
       </Position19>
       <Position20>
         <Filename Value="romunit.pas"/>
-        <Caret Line="1006" Column="22" TopLine="984"/>
+        <Caret Line="30" Column="14" TopLine="12"/>
       </Position20>
       <Position21>
         <Filename Value="romunit.pas"/>
-        <Caret Line="1010" Column="43" TopLine="984"/>
+        <Caret Line="737" Column="66" TopLine="708"/>
       </Position21>
       <Position22>
-        <Filename Value="romunit.pas"/>
-        <Caret Line="1011" Column="58" TopLine="984"/>
+        <Filename Value="ROM.lpr"/>
+        <Caret Line="312" Column="7" TopLine="294"/>
       </Position22>
       <Position23>
-        <Filename Value="ROM.lpr"/>
-        <Caret Line="13" Column="62"/>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="13" Column="22"/>
       </Position23>
       <Position24>
-        <Filename Value="ROM.lpr"/>
-        <Caret Line="282" Column="20" TopLine="262"/>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="14" Column="22"/>
       </Position24>
       <Position25>
-        <Filename Value="ROM.lpr"/>
-        <Caret Line="279" Column="16" TopLine="266"/>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="197" Column="80" TopLine="181"/>
       </Position25>
       <Position26>
         <Filename Value="romunit.pas"/>
-        <Caret Line="30" Column="14" TopLine="23"/>
       </Position26>
       <Position27>
-        <Filename Value="ROM.lpr"/>
-        <Caret Line="279" Column="16" TopLine="266"/>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="31" Column="14" TopLine="2"/>
       </Position27>
       <Position28>
         <Filename Value="ROM.lpr"/>
-        <Caret Line="17" Column="28"/>
+        <Caret Line="29" Column="60" TopLine="18"/>
       </Position28>
       <Position29>
-        <Filename Value="ROM.lpr"/>
-        <Caret Line="265" Column="26" TopLine="248"/>
+        <Filename Value="romunit.pas"/>
+        <Caret Line="294" Column="12" TopLine="277"/>
       </Position29>
       <Position30>
         <Filename Value="romunit.pas"/>
-        <Caret Line="25" Column="19" TopLine="7"/>
+        <Caret Line="31" Column="14" TopLine="13"/>
       </Position30>
     </JumpHistory>
   </ProjectSession>
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);