Changeset 40146 for branches/eam_branches/ohana.20170822/src/opihi/cmd.data/periodogram-fm.c

Timestamp:

Sep 26, 2017, 7:03:10 PM (9 years ago)

Author:

eugene

Message:

adding the periodogram - floating means

File:

• branches/eam_branches/ohana.20170822/src/opihi/cmd.data/periodogram-fm.c (modified) (7 diffs)

branches/eam_branches/ohana.20170822/src/opihi/cmd.data/periodogram-fm.c

@@ -1 +1 @@
 # 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 +12 @@
 
   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 +18 @@
   // 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 +53 @@
   }
 
-  Np = 0;
-  NP = 100;
+  int Np = 0;
+  int NP = 100;
   ResetVector (power,  OPIHI_FLT, NP);
   ResetVector (period, OPIHI_FLT, NP);
@@ -70 +63 @@
 
   // 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 +71 @@
     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 +133 @@
     }
 
-    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);