Changeset 40308

Timestamp:

Jan 1, 2018, 9:36:48 AM (9 years ago)

Author:

eugene

Message:

add options to periodogram-fm

Location:

trunk/Ohana/src/opihi/cmd.data

Files:

• periodogram-fm.c (modified) (4 diffs)
• test/periodogram-fm.sh (modified) (1 diff)

trunk/Ohana/src/opihi/cmd.data/periodogram-fm.c

@@ -8 +8 @@
   if ((N = get_argument (argc, argv, "-v"))) {
     VERBOSE = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
+  int LINEAR = FALSE;
+  if ((N = get_argument (argc, argv, "-linear"))) {
+    LINEAR = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
+  int Nperiods = 1024;
+  if ((N = get_argument (argc, argv, "-Nperiods"))) {
+    remove_argument (N, &argc, argv);
+    Nperiods = atoi(argv[N]);
     remove_argument (N, &argc, argv);
   }
@@ -54 +67 @@
 
   int Np = 0;
-  int NP = 100;
+  int NP = Nperiods;
   ResetVector (power,  OPIHI_FLT, NP);
   ResetVector (period, OPIHI_FLT, NP);
 
   // for testing, we are going to write out the terms explicitly 
-  // for optimization, some of the tri functions can be calculated more efficiently 
+  // for optimization, some of the trig functions can be calculated more efficiently 
   // by storing cos,sin(w t) and using some recurrence rules
+
+  // maxP = minP * dP^n [n = 0 -- Nperiods-1]
+  // log(maxP) = log(minP) + (Nperiods - 1) * log (dP)
+  // log(dP) = (log(maxP) - log(minP)) / (Nperiods - 1)
+  double dP = LINEAR ? (maxP - minP)/(Nperiods - 1) : pow(10.0, ((log10(maxP) - log10(minP))/(Nperiods - 1)));
 
   // XXX consider choice of period or frequency step
   double P = minP;
-  while (P <= maxP) {
+  for (Np = 0; Np < Nperiods; Np++) {
     double w = 2*M_PI/P;
     
@@ -126 +144 @@
     power[0].elements.Flt[Np] = Po;
     period[0].elements.Flt[Np] = P;
-    Np ++;
     if (Np >= NP) {
       NP += 100;
@@ -133 +150 @@
     }
 
-    double ratio = 1 + 0.1*P/dTime;
+    // double ratio = 1 + 0.1*P/dTime;
 
-    if (VERBOSE) gprint (GP_ERR, "tau: %f, P: %f, ratio: %f, dTime: %f, nextP: %f\n", tau, P, ratio, dTime, P*ratio);
+    if (VERBOSE) gprint (GP_ERR, "tau: %f, P: %f, ratio: %f, dTime: %f, nextP: %f\n", tau, P, dP, dTime, P*dP);
 
-    P *= ratio;
+    if (LINEAR) {
+      P += dP;
+    } else {
+      P *= dP;
+    }
   }
 

trunk/Ohana/src/opihi/cmd.data/test/periodogram-fm.sh

@@ -25 +25 @@
 
  periodogram_fm t f df 5 50 period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.5)
+   $PASS = 0
+   echo "OFFSET: {$peakpos - $P}"
+ end
+
+ if ($PLOT)
+  lim period power; clear; box; line -c red70 -lw 3 $P 0 to $P $peakval; plot period power -x line
+ end
+end
+
+# test using even samples
+macro test1a
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 15.0
+
+ delete -q t f period power
+
+ create t 0 100
+ set f = sin(2*$PI*t/$P)
+ set df = 0.01 + zero(f)
+
+ periodogram_fm t f df 5 50 period power -Nperiods 512
 
  peak -q period power