Index: /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/Makefile
===================================================================
--- /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/Makefile	(revision 40145)
+++ /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/Makefile	(revision 40146)
@@ -102,4 +102,5 @@
 $(SRC)/peak.$(ARCH).o		\
 $(SRC)/periodogram.$(ARCH).o	\
+$(SRC)/periodogram-fm.$(ARCH).o	\
 $(SRC)/plot.$(ARCH).o		\
 $(SRC)/dot.$(ARCH).o		\
Index: /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/init.c
===================================================================
--- /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/init.c	(revision 40145)
+++ /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/init.c	(revision 40146)
@@ -91,4 +91,5 @@
 int peak             PROTO((int, char **));
 int periodogram      PROTO((int, char **));
+int periodogram_fm   PROTO((int, char **));
 int plot             PROTO((int, char **));
 int dot              PROTO((int, char **));
@@ -270,5 +271,6 @@
   {1, "parity",       parity,           "set image parity"},
   {1, "peak",         peak,             "find vector peak in range"},
-  {1, "periodogram",  periodogram,      "measure periods in unevenly sampled data"},
+  {1, "periodogram",    periodogram,    "measure periods in unevenly sampled data (Lomb-Scargle)"},
+  {1, "periodogram_fm", periodogram_fm, "measure periods in unevenly sampled data (generalized Lomb-Scargle; floating mean)"},
   {1, "plot",         plot,             "plot a pair of vectors"},
   {1, "png",          jpeg,             "convert display graphic to PNG"},
Index: /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/periodogram-fm.c
===================================================================
--- /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/periodogram-fm.c	(revision 40145)
+++ /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/periodogram-fm.c	(revision 40146)
@@ -1,14 +1,9 @@
 # include "data.h"
+# define MIN_VAR 1e-8
+int periodogram_fm (int argc, char **argv) {
+  
+  int N;
 
-int periodogram (int argc, char **argv) {
-  
-  int i, N, Npt, Np, NP, VERBOSE;
-  opihi_flt *tv, *fv;
-  float minP, maxP, minT, maxT, dTime;
-  float mean, var, w, tau, P, Pc, Ps, Po;
-  float C, S, cs, sn, cs2, sn2, ratio;
-  Vector *time, *flux, *dflux, *power, *period;
-
-  VERBOSE = FALSE;
+  int VERBOSE = FALSE;
   if ((N = get_argument (argc, argv, "-v"))) {
     VERBOSE = TRUE;
@@ -17,5 +12,5 @@
 
   if (argc != 8) {
-    gprint (GP_ERR, "USAGE: periodogram-fm (time) (flux) (dflux) (minP) (maxP) (period) (power)\n");
+    gprint (GP_ERR, "USAGE: periodogram_fm (time) (flux) (dflux) (minP) (maxP) (period) (power)\n");
     return (FALSE);
   }
@@ -23,36 +18,34 @@
   // XXX allow dflux to be dropped?
 
-  if ((time  = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((flux  = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((dflux = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
-  minP = atof(argv[4]);
-  maxP = atof(argv[5]);
+  Vector *time, *flux, *dflux, *power, *period;
+  if ((time   = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((flux   = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((dflux  = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
   if ((period = SelectVector (argv[6], ANYVECTOR, TRUE)) == NULL) return (FALSE);
   if ((power  = SelectVector (argv[7], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  double minP = atof(argv[4]);
+  double maxP = atof(argv[5]);
 
   REQUIRE_VECTOR_FLT (time, FALSE); 
   REQUIRE_VECTOR_FLT (flux, FALSE); 
+  REQUIRE_VECTOR_FLT (dflux, FALSE); 
 
-  /* find the max baseline, mean, and variance */
-  minT = maxT = time[0].elements.Flt[0];
-  Npt = time[0].Nelements;
-  tv = time[0].elements.Flt;
-  fv = flux[0].elements.Flt;
-  mean = var = 0;
-  for (i = 0; i < Npt; i++, tv++, fv++) {
+  /* find the max baseline, sum the inverse variances */
+  double minT = time[0].elements.Flt[0];
+  double maxT = time[0].elements.Flt[0];
+  opihi_flt Weight = 0;
+  int Npt = time[0].Nelements;
+
+  opihi_flt *tv =  time[0].elements.Flt;
+  opihi_flt *df = dflux[0].elements.Flt;
+  for (int i = 0; i < Npt; i++, tv++, df++) {
     minT = MIN (minT, *tv);
     maxT = MAX (maxT, *tv);
-    mean += *fv;
+    // skip points with 0.0 error? or add minimum variance?
+    Weight += 1.0 / (SQ(*df) + MIN_VAR); // MIN_VAR : added in quadrature to avoid Inf (XXX make this a user parameter?)
   }
-  mean = mean / Npt;
-  fv = flux[0].elements.Flt;
-  for (i = 0; i < Npt; i++, fv++) {
-    var += SQ(*fv - mean);
-  }
-  var = var / (Npt - 1);
 
-  if (VERBOSE) gprint (GP_ERR, "mean: %f, var: %f, minT: %f, maxT: %f\n", mean, var, minT, maxT);
-
-  dTime = maxT - minT;
+  double WeightInv = 1.0 / Weight;
+  double dTime = maxT - minT;
   if (dTime == 0) {
     gprint (GP_ERR, "ERROR: time range is zero\n");
@@ -60,6 +53,6 @@
   }
 
-  Np = 0;
-  NP = 100;
+  int Np = 0;
+  int NP = 100;
   ResetVector (power,  OPIHI_FLT, NP);
   ResetVector (period, OPIHI_FLT, NP);
@@ -70,7 +63,7 @@
 
   // XXX consider choice of period or frequency step
-  P = minP;
-  while (P < maxP) {
-    w = 2*M_PI/P;
+  double P = minP;
+  while (P <= maxP) {
+    double w = 2*M_PI/P;
     
     /* find the period offset tau  */
@@ -78,30 +71,56 @@
     df = dflux[0].elements.Flt;
 
-    cs = sn = cs2 = sn2 = 0;
+    double cs = 0.0, sn = 0.0, sn2 = 0.0, cs2 = 0.0;
 
-    for (i = 0; i < Npt; i++, tv++) {
-      opihi_flt wt = 
-      cs += cos (*tv*w*2);
-      sn += sin (*tv*2*2);
+    for (int i = 0; i < Npt; i++, tv++, df++) {
+      opihi_flt wt = WeightInv / (SQ(*df) + MIN_VAR);
+      cs  += wt*cos (*tv*w);
+      sn  += wt*sin (*tv*w);
+      cs2 += wt*cos (*tv*w*2);
+      sn2 += wt*sin (*tv*w*2);
     }
-    tau = 0.5*atan2 (sn, cs) / w;
+    double ytan = sn2 - 2*cs*sn;
+    double xtan = cs2 - (SQ(cs) - SQ(sn));
+    double tau = 0.5*atan2 (ytan, xtan) / w;
       
     /* find the power at this period */
     tv = time[0].elements.Flt;
-    fv = flux[0].elements.Flt;
-    cs = sn = cs2 = sn2 = 0;
-    for (i = 0; i < Npt; i++, tv++, fv++) {
-      C = cos (w*(*tv-tau));
-      S = sin (w*(*tv-tau));
-      // C = cos (w**tv);
-      // S = sin (w**tv);
-      cs += (*fv - mean) * C;
-      sn += (*fv - mean) * S;
-      cs2 += SQ(C);
-      sn2 += SQ(S);
+    df = dflux[0].elements.Flt;
+    opihi_flt *fv = flux[0].elements.Flt;
+
+    // YY = YY_s - Y_s*Y_s
+    // CC = CC_s - C_s*C_s
+    // SS = SS_s - S_s*S_s
+    // YC = YC_s - Y_s*C_s
+    // YS = YS_s - Y_s*S_s
+    
+    double YY_s = 0.0, CC_s = 0.0, SS_s = 0.0, YC_s = 0.0, YS_s = 0.0;
+    double Y_s = 0.0, C_s = 0.0, S_s = 0.0;
+    for (int i = 0; i < Npt; i++, tv++, fv++, df++) {
+      opihi_flt wt = WeightInv / (SQ(*df) + MIN_VAR);
+      cs = cos (w*(*tv-tau));
+      sn = sin (w*(*tv-tau));
+
+      double wtf = *fv*wt;
+      Y_s += wtf;
+      C_s += wt*cs;
+      S_s += wt*sn;
+
+      YY_s += wtf*(*fv);
+      YC_s += wtf*cs;
+      YS_s += wtf*sn;
+
+      CC_s += wt*cs*cs;
+      SS_s += wt*sn*sn;
     }
-    Pc = SQ(cs) / cs2;
-    Ps = SQ(sn) / sn2;
-    Po = (Pc + Ps) / (2*var);
+    double YY = YY_s - SQ(Y_s);
+    double YC = YC_s - Y_s*C_s;
+    double YS = YS_s - Y_s*S_s;
+    double CC = CC_s - SQ(C_s);
+    double SS = SS_s - SQ(S_s);
+
+    double Pa = SQ(YC) / CC;
+    double Pb = SQ(YS) / SS;
+    double Po = (Pa + Pb) / YY;
 
     power[0].elements.Flt[Np] = Po;
@@ -114,5 +133,5 @@
     }
 
-    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);
Index: /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/test/periodogram-fm.sh
===================================================================
--- /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/test/periodogram-fm.sh	(revision 40146)
+++ /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/test/periodogram-fm.sh	(revision 40146)
@@ -0,0 +1,313 @@
+
+if (not($?PLOT)) set PLOT = 0
+
+list tests
+ test1
+ test2
+ test3
+ memtest1
+end
+
+# test using even samples
+macro test1
+ $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
+
+ 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 random samples
+macro test2
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 15.0
+
+ delete -q x t f period power
+
+ create x 0 100
+ set t = 100 * rnd(x) 
+ set f = sin(2*$PI*t/$P)
+ set df = 0.01 + zero(f)
+
+ periodogram_fm t f df 5 50 period power
+#periodogram t f 5 50 period power
+
+ # lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
+ # lim -n 1 period power; clear; box; plot period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.5)
+   $PASS = 0
+ 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 random samples, higher frequency
+macro test3
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 2.0
+
+ delete -q x t f period power
+
+ create x 0 100
+ set t = 100 * rnd(x) 
+ set f = sin(2*$PI*t/$P)
+ set df = 0.01 + zero(f)
+
+ periodogram_fm t f df 1 10 period power
+#periodogram t f 1 10 period power
+
+ # lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
+ # lim -n 1 period power; clear; box; plot period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.05)
+   $PASS = 0
+ end
+end
+
+# test using random samples, offset start
+macro test4
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 15.0
+
+ delete -q x t f period power
+
+ create x 500 800
+ set t = 300 * rnd(x) + 500
+ set f = sin(2*$PI*t/$P)
+ set df = 0.01 + zero(f)
+
+ periodogram_fm t f df 2 30 period power
+#periodogram t f 2 30 period power
+
+#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
+#  lim -n 1 period power; clear; box; plot period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.05)
+   $PASS = 0
+ end
+end
+
+# test using random samples, offset start, non-zero DC
+macro test5
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 15.0
+
+ delete -q x t f period power
+
+ create x 500 800
+ set t = 300 * rnd(x) + 500
+ set f = sin(2*$PI*t/$P) + 0.5
+ set df = 0.01 + zero(f)
+
+ periodogram_fm t f 2 30 period power
+
+#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
+#  lim -n 1 period power; clear; box; plot period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.05)
+   $PASS = 0
+ end
+end
+
+# test using random samples, offset start, non-zero DC, some noise
+macro test6
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 15.0
+
+ delete -q x t f period power
+
+ create x 500 800
+ set t = 300 * rnd(x) + 500
+ set fraw = sin(2*$PI*t/$P) + 0.5
+
+ # 0.05 : peakpos = 14.95
+ # 0.10 : peakpos = 15.04 (
+ gaussdev df t[] 0.0 0.25
+ set f = fraw + df
+ set df = 0.01 + zero(f)
+
+ periodogram_fm t f 2 30 period power
+
+#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
+#  lim -n 1 period power; clear; box; plot period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.05)
+   $PASS = 0
+ end
+end
+
+# test using fewer random samples, offset start, non-zero DC, some noise
+macro test7
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 15.0
+
+ delete -q x t f period power
+
+ create x 0 100
+ set t = 100 * rnd(x)
+ set fraw = sin(2*$PI*t/$P) + 0.5
+
+ # 0.05 : peakpos = 14.95
+ # 0.10 : peakpos = 15.04 (
+ gaussdev df t[] 0.0 0.25
+ set f = fraw + df
+ set df = 0.01 + zero(f)
+
+ periodogram_fm t f 2 30 period power
+
+#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
+#  lim -n 1 period power; clear; box; plot period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.05)
+   $PASS = 0
+ end
+end
+
+# test using fewer random samples, high frequency, non-zero DC, some noise
+macro test8
+ if ($0 != 2)
+   echo "USAGE: test8: Ndays")
+   break
+ end
+ 
+ local Ndays 
+ $Ndays = $1
+
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 0.8*rnd(0) + 0.2
+ $trueP = $P
+
+ delete -q x t f period power
+
+ create x 0 $Ndays
+
+ # t is a time in days, but we always have 4 within 1 hour:
+ set t0 = int(100 * rnd(x))
+ set dtx = (3/24) * rnd(x)
+ set t0 = t0 + dtx
+
+ set dt1 = (15.0 / 1440) * rnd(x) + ( 0 + 7.5) / 1440
+ set dt2 = (15.0 / 1440) * rnd(x) + (15 + 7.5) / 1440
+ set dt3 = (15.0 / 1440) * rnd(x) + (30 + 7.5) / 1440
+
+ delete -q t
+ concat t0 t
+ set tmp = t0 + dt1; concat tmp t
+ set tmp = t0 + dt2; concat tmp t
+ set tmp = t0 + dt3; concat tmp t
+
+ set fraw = sin(2*$PI*t/$P) + 0.5
+
+ # 0.05 : peakpos = 14.95
+ # 0.10 : peakpos = 15.04 (
+ gaussdev df t[] 0.0 0.25
+ set f = fraw + df
+ set df = 0.01 + zero(f)
+
+ periodogram_fm t f 0.1 2.0 period power
+
+#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
+#  lim -n 1 period power; clear; box; plot period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.05)
+   $PASS = 0
+ end
+end
+
+# Memory test
+macro memtest1
+
+ local i
+ local P PI
+ $PI = 3.14159265359
+ $P  = 15.0
+
+ delete -q x t f period power
+
+ create x 500 800
+ set t = 300 * rnd(x) + 500
+ set f = sin(2*$PI*t/$P)
+
+ list word -x "ps -p $PID -o rss"
+ $startmem = $word:1
+
+ for i 0 100
+  periodogram_fm t f 2 30 period power
+ end
+  
+ list word -x "ps -p $PID -o rss"
+ $endmem = $word:1
+
+ $PASS = 1
+
+ if ($endmem - $startmem > 180)
+   $PASS = 0
+   echo "growth: {$endmem-$startmem}"
+   echo "kB/loop: {($endmem-$startmem)/100}"
+ end
+end
Index: /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/test/periodogram.sh
===================================================================
--- /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/test/periodogram.sh	(revision 40145)
+++ /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/test/periodogram.sh	(revision 40146)
@@ -1,2 +1,4 @@
+
+if (not($?PLOT)) set PLOT = 0
 
 list tests
@@ -28,4 +30,8 @@
    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