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

Timestamp:

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

Author:

eugene

Message:

add options to periodogram-fm

File:

• trunk/Ohana/src/opihi/cmd.data/periodogram-fm.c (modified) (4 diffs)

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;
+    }
   }