[xiph-commits] r17928 - trunk/ghost/monty/chirp

xiphmont at svn.xiph.org xiphmont at svn.xiph.org
Thu Apr 21 03:09:00 PDT 2011


Author: xiphmont
Date: 2011-04-21 03:09:00 -0700 (Thu, 21 Apr 2011)
New Revision: 17928

Modified:
   trunk/ghost/monty/chirp/chirp.c
   trunk/ghost/monty/chirp/chirpgraph.c
   trunk/ghost/monty/chirp/chirpgraph.h
   trunk/ghost/monty/chirp/chirptest.c
Log:
Newest graphing code; sweep any dimension



Modified: trunk/ghost/monty/chirp/chirp.c
===================================================================
--- trunk/ghost/monty/chirp/chirp.c	2011-04-20 04:54:38 UTC (rev 17927)
+++ trunk/ghost/monty/chirp/chirp.c	2011-04-21 10:09:00 UTC (rev 17928)
@@ -283,9 +283,9 @@
 
         /* guard overflow; if we're this far out, assume we're never
            coming back. drop out now. */
-        if((aP*aP + bP*bP)>1e20 ||
-           (cP*cP + dP*dP)>1e20 ||
-           (eP*eP + fP*fP)>1e20){
+        if((aP*aP + bP*bP)>1e10 ||
+           (cP*cP + dP*dP)>1e10 ||
+           (eP*eP + fP*fP)>1e10){
           iter_limit=0;
           i=n;
         }

Modified: trunk/ghost/monty/chirp/chirpgraph.c
===================================================================
--- trunk/ghost/monty/chirp/chirpgraph.c	2011-04-20 04:54:38 UTC (rev 17927)
+++ trunk/ghost/monty/chirp/chirpgraph.c	2011-04-21 10:09:00 UTC (rev 17928)
@@ -79,8 +79,8 @@
 /********* draw everything in the graph except the graph data itself *******/
 
 static float fontsize=12;
-static int x0s,x1s,xmajor;
-static int y0s,y1s,ymajor;
+static int x0s,x1s,xmajor,xmajorf;
+static int y0s,y1s,ymajor,ymajorf;
 
 /* computed configuration globals */
 int leftpad=0;
@@ -100,10 +100,12 @@
 void setup_graphs(int start_x_step,
                   int end_x_step, /* inclusive; not one past */
                   int x_major_d,
+                  float x_major_f,
 
                   int start_y_step,
                   int end_y_step, /* inclusive; not one past */
                   int y_major_d,
+                  float y_major_f,
 
                   int subtitles,
                   float fs){
@@ -121,16 +123,18 @@
   x0s=start_x_step;
   x1s=end_x_step;
   xmajor=x_major_d;
+  xmajorf=x_major_f;
   y0s=start_y_step;
   y1s=end_y_step;
   ymajor=y_major_d;
+  ymajorf=y_major_f;
 
   /* y labels */
   cairo_set_font_size(ct, fontsize);
   for(y=y0s+fontsize;y<=y1s;y++){
     if(y%ymajor==0){
       char buf[80];
-      snprintf(buf,80,"%.0f",(float)y/ymajor);
+      snprintf(buf,80,"%.0f",(float)y/ymajor*ymajorf);
       cairo_text_extents(ct, buf, &extents);
       if(extents.width + extents.height*3.5>leftpad)leftpad=extents.width + extents.height*3.5;
     }
@@ -229,7 +233,7 @@
       int y = toppad+y_n-i+y0s;
       char buf[80];
 
-      snprintf(buf,80,"%.0f",(float)i/ymajor);
+      snprintf(buf,80,"%.0f",(float)i/ymajor*ymajorf);
       cairo_text_extents(c, buf, &extents);
       cairo_move_to(c,leftpad - fontsize*.5 - extents.width,y+extents.height*.5);
       cairo_show_text(c,buf);
@@ -237,28 +241,22 @@
   }
 
   /* X axis labels */
-  {
-    int xadv=0;
-    for(i=x0s;i<=x1s;i++){
-      if(i%xmajor==0){
-        char buf[80];
-        int x = leftpad + i - x0s;
-        if(i==0){
-          snprintf(buf,80,"DC");
-        }else{
-          snprintf(buf,80,"%.0f",(float)i/xmajor);
-        }
-        cairo_text_extents(c, buf, &extents);
-        if(x-extents.width/2 - fontsize > xadv){
-          cairo_move_to(c,x - extents.width/2,y_n+toppad+extents.height+fontsize*.5);
-          cairo_show_text(c,buf);
-          xadv = x+extents.width/2;
-        }
+  for(i=x0s;i<=x1s;i++){
+    if(i%xmajor==0){
+      char buf[80];
+      int x = leftpad + i - x0s;
+      if(i==0){
+        snprintf(buf,80,"DC");
+      }else{
+        snprintf(buf,80,"%.0f",(float)i/xmajor*xmajorf);
       }
+      cairo_text_extents(c, buf, &extents);
+      cairo_move_to(c,x - extents.width/2,y_n+toppad+extents.height+fontsize*.5);
+      cairo_show_text(c,buf);
     }
   }
 
-  /* Y axis label */
+  /* Y axis caption */
   {
     cairo_matrix_t a;
     cairo_matrix_t b = {0.,-1., 1.,0., 0.,0.}; // account for border!

Modified: trunk/ghost/monty/chirp/chirpgraph.h
===================================================================
--- trunk/ghost/monty/chirp/chirpgraph.h	2011-04-20 04:54:38 UTC (rev 17927)
+++ trunk/ghost/monty/chirp/chirpgraph.h	2011-04-21 10:09:00 UTC (rev 17928)
@@ -37,10 +37,12 @@
 extern void setup_graphs(int start_x_step,
                          int end_x_step, /* inclusive; not one past */
                          int x_major_d,
+                         float x_major_f,
 
                          int start_y_step,
                          int end_y_step, /* inclusive; not one past */
                          int y_major_d,
+                         float y_major_f,
 
                          int subtitles,
                          float fs);

Modified: trunk/ghost/monty/chirp/chirptest.c
===================================================================
--- trunk/ghost/monty/chirp/chirptest.c	2011-04-20 04:54:38 UTC (rev 17927)
+++ trunk/ghost/monty/chirp/chirptest.c	2011-04-21 10:09:00 UTC (rev 17928)
@@ -35,39 +35,113 @@
   return ret;
 }
 
+#define DIM_ESTIMATE_A 1
+#define DIM_ESTIMATE_P 2
+#define DIM_ESTIMATE_W 3
+#define DIM_ESTIMATE_dA 4
+#define DIM_ESTIMATE_dW 5
+#define DIM_ESTIMATE_ddA 6
+#define DIM_ESTIMATE_MASK 0xf
+
+#define DIM_CHIRP_A (1<<4)
+#define DIM_CHIRP_P (2<<4)
+#define DIM_CHIRP_W (3<<4)
+#define DIM_CHIRP_dA (4<<4)
+#define DIM_CHIRP_dW  (5<<4)
+#define DIM_CHIRP_ddA (6<<4)
+#define DIM_CHIRP_MASK 0xf0
+
 /* ranges are inclusive */
-void set_chirp(chirp *c,int rand_p, int i, int n,
+void set_chirp(chirp *c,
+               int xi, int xn, int xdim,
+               int yi, int yn, int ydim,
+               int stepi, int stepn, int rand_p,
                float A0, float A1,
                float P0, float P1,
                float W0, float W1,
                float dA0, float dA1,
                float dW0, float dW1,
                float ddA0, float ddA1){
-  if(n<=1){
-    c->A=A0;
-    c->P=P0;
-    c->W=W0;
-    c->dA=dA0;
-    c->dW=dW0;
-    c->ddA=ddA0;
-  }else if(rand_p){
-    c->A = A0 + (A1-A0)*drand48();
-    c->P = P0 + (P1-P0)*drand48();
-    c->W = W0 + (W1-W0)*drand48();
-    c->dA = dA0 + (dA1-dA0)*drand48();
-    c->dW = dW0 + (dW1-dW0)*drand48();
-    c->ddA = ddA0 + (ddA1-ddA0)*drand48();
-  }else{
-    c->A = A0 + (A1-A0)/(n-1)*i;
-    c->P = P0 + (P1-P0)/(n-1)*i;
-    c->W = W0 + (W1-W0)/(n-1)*i;
-    c->dA = dA0 + (dA1-dA0)/(n-1)*i;
-    c->dW = dW0 + (dW1-dW0)/(n-1)*i;
-    c->ddA = ddA0 + (ddA1-ddA0)/(n-1)*i;
+
+  int An,Pn,Wn,dAn,dWn,ddAn;
+  int Ai,Pi,Wi,dAi,dWi,ddAi;
+
+  xdim = (xdim | (xdim>>4)) & DIM_ESTIMATE_MASK;
+  ydim = (ydim | (ydim>>4)) & DIM_ESTIMATE_MASK;
+  if(stepn<2)stepn=2;
+
+  An=Pn=Wn=dAn=dWn=ddAn=stepn-1;
+  Ai = (rand_p ? drand48() : stepi);
+  Pi = (rand_p ? drand48() : stepi);
+  Wi = (rand_p ? drand48() : stepi);
+  dAi = (rand_p ? drand48() : stepi);
+  dWi = (rand_p ? drand48() : stepi);
+  ddAi = (rand_p ? drand48() : stepi);
+
+  switch(xdim){
+  case DIM_ESTIMATE_A:
+    An = xn-1;
+    Ai = xi;
+    break;
+  case DIM_ESTIMATE_P:
+    Pn = xn-1;
+    Pi = xi;
+    break;
+  case DIM_ESTIMATE_W:
+    Wn = xn-1;
+    Wi = xi;
+    break;
+  case DIM_ESTIMATE_dA:
+    dAn = xn-1;
+    dAi = xi;
+    break;
+  case DIM_ESTIMATE_dW:
+    dWn = xn-1;
+    dWi = xi;
+    break;
+  case DIM_ESTIMATE_ddA:
+    ddAn = xn-1;
+    ddAi = xi;
+    break;
   }
+
+  switch(ydim){
+  case DIM_ESTIMATE_A:
+    An = yn-1;
+    Ai = yi;
+    break;
+  case DIM_ESTIMATE_P:
+    Pn = yn-1;
+    Pi = yi;
+    break;
+  case DIM_ESTIMATE_W:
+    Wn = yn-1;
+    Wi = yi;
+    break;
+  case DIM_ESTIMATE_dA:
+    dAn = yn-1;
+    dAi = yi;
+    break;
+  case DIM_ESTIMATE_dW:
+    dWn = yn-1;
+    dWi = yi;
+    break;
+  case DIM_ESTIMATE_ddA:
+    ddAn = yn-1;
+    ddAi = yi;
+    break;
+  }
+
+  c->A = A0 + (A1-A0) / An * Ai;
+  c->P = P0 + (P1-P0) / Pn * Pi;
+  c->W = W0 + (W1-W0) / Wn * Wi;
+  c->dA = dA0 + (dA1-dA0) / dAn * dAi;
+  c->dW = dW0 + (dW1-dW0) / dWn * dWi;
+  c->ddA = ddA0 + (ddA1-ddA0) / ddAn * ddAi;
+
 }
 
-/*********************** Plot w estimate error against w **********************/
+/*********************** Plot single estimate vs. chirp **********************/
 
 typedef struct {
   float *in;
@@ -112,7 +186,8 @@
   float *chirp = localinit ? malloc(sizeof(*chirp)*blocksize) : arg->in;
 
   while(1){
-    float s=0.;
+    float rms_acc=0.;
+    float e_acc=0.;
 
     pthread_mutex_lock(&ymutex);
     y=next_y;
@@ -129,8 +204,8 @@
     if(localinit){
       for(i=0;i<blocksize;i++){
         double jj = i - blocksize/2 + .5;
-        double A = arg->chirp[y].A + arg->chirp[y].dA*jj + arg->chirp[y].ddA*jj*jj;
-        double P = arg->chirp[y].P + arg->chirp[y].W *jj + arg->chirp[y].dW *jj*jj;
+        double A = arg->chirp[y].A + (arg->chirp[y].dA + arg->chirp[y].ddA*jj)*jj;
+        double P = arg->chirp[y].P + (arg->chirp[y].W  + arg->chirp[y].dW *jj)*jj;
         chirp[i] = A*cos(P);
       }
     }
@@ -153,10 +228,12 @@
                         arg->fit_bound_zero);
 
     for(i=0;i<blocksize;i++){
-      float v = (chirp[i]-rec[i])*arg->window[i];
-      s += v*v;
+      float ce = chirp[i]*arg->window[i];
+      float ee = (chirp[i]-rec[i])*arg->window[i];
+      e_acc += ce*ce;
+      rms_acc += ee*ee;
     }
-    arg->rms_error[y] = sqrt(s/blocksize);
+    arg->rms_error[y] = sqrt(rms_acc)/sqrt(e_acc);
     arg->iterations[y] = arg->max_iterations-ret;
     feclearexcept(FE_ALL_EXCEPT);
     feenableexcept(except);
@@ -172,19 +249,12 @@
   char *subtitle1;
   char *subtitle2;
   char *subtitle3;
+  char *xaxis_label;
+  char *yaxis_label;
 
   int blocksize;
   int threads;
 
-  int x0;
-  int x1;
-  int xmajor;
-  int xminor;
-  int y0;
-  int y1;
-  int ymajor;
-  int yminor;
-
   void (*window)(float *,int n);
   float fit_tolerance;
 
@@ -199,30 +269,45 @@
   int fit_symm_norm;
   int fit_bound_zero;
 
+  int x_dim;
+  int x_steps;
+  float x_major;
+  float x_minor;
+  int y_dim;
+  int y_steps;
+  float y_major;
+  float y_minor;
+  int sweep_steps;
+
   /* If the randomize flag is unset and min!=max, a param is swept
      from min to max (inclusive) divided into <sweep_steps>
      increments.  If the rand flag is set, <sweep_steps> random values
      in the range min to max instead. */
-  int sweep_steps;
   int sweep_or_rand_p;
 
   float min_est_A;
   float max_est_A;
+  int rel_est_A;
 
   float min_est_P;
   float max_est_P;
+  int rel_est_P;
 
   float min_est_W;
   float max_est_W;
+  int rel_est_W;
 
   float min_est_dA;
   float max_est_dA;
+  int rel_est_dA;
 
   float min_est_dW;
   float max_est_dW;
+  int rel_est_dW;
 
   float min_est_ddA;
   float max_est_ddA;
+  int rel_est_ddA;
 
   float min_chirp_A;
   float max_chirp_A;
@@ -242,6 +327,31 @@
   float min_chirp_ddA;
   float max_chirp_ddA;
 
+  /* generate which graphs? */
+  int graph_convergence_max;
+  int graph_convergence_delta;
+
+  int graph_Aerror_max;
+  int graph_Aerror_delta;
+
+  int graph_Perror_max;
+  int graph_Perror_delta;
+
+  int graph_Werror_max;
+  int graph_Werror_delta;
+
+  int graph_dAerror_max;
+  int graph_dAerror_delta;
+
+  int graph_dWerror_max;
+  int graph_dWerror_delta;
+
+  int graph_ddAerror_max;
+  int graph_ddAerror_delta;
+
+  int graph_RMSerror_max;
+  int graph_RMSerror_delta;
+
 }  graph_run;
 
 /* performs a W initial estimate error vs chirp W plot.  Ignores the
@@ -254,8 +364,8 @@
   float in[blocksize];
   int i,xi,yi;
 
-  int x_n = arg->x1-arg->x0+1;
-  int y_n = arg->y1-arg->y0+1;
+  int x_n = arg->x_steps;
+  int y_n = arg->y_steps;
 
   /* graphs:
 
@@ -294,8 +404,6 @@
   pthread_t threadlist[threads];
   colarg targ[threads];
 
-  char *yaxis_label = "initial distance from W (cycles/block)";
-  char *xaxis_label = "W (cycles/block)";
   int swept=0;
   int chirp_swept=0;
   int est_swept=0;
@@ -318,182 +426,336 @@
   float ret_maxRMS[y_n];
   float ret_maxiter[y_n];
 
+  float minX=0,maxX=0;
+  float minY=0,maxY=0;
+  int x0s,x1s;
+  int y0s,y1s;
+  int xmajori,ymajori;
+  int xminori,yminori;
+
   struct timeval last;
   gettimeofday(&last,NULL);
 
+  switch(arg->x_dim){
+  case DIM_ESTIMATE_A:
+    minX = arg->min_est_A;
+    maxX = arg->max_est_A;
+    break;
+  case DIM_ESTIMATE_P:
+    minX = arg->min_est_P;
+    maxX = arg->max_est_P;
+    break;
+  case DIM_ESTIMATE_W:
+    minX = arg->min_est_W;
+    maxX = arg->max_est_W;
+    break;
+  case DIM_ESTIMATE_dA:
+    minX = arg->min_est_dA;
+    maxX = arg->max_est_dA;
+    break;
+  case DIM_ESTIMATE_dW:
+    minX = arg->min_est_dW;
+    maxX = arg->max_est_dW;
+    break;
+  case DIM_ESTIMATE_ddA:
+    minX = arg->min_est_ddA;
+    maxX = arg->max_est_ddA;
+    break;
+  case DIM_CHIRP_A:
+    minX = arg->min_chirp_A;
+    maxX = arg->max_chirp_A;
+    break;
+  case DIM_CHIRP_P:
+    minX = arg->min_chirp_P;
+    maxX = arg->max_chirp_P;
+    break;
+  case DIM_CHIRP_W:
+    minX = arg->min_chirp_W;
+    maxX = arg->max_chirp_W;
+    break;
+  case DIM_CHIRP_dA:
+    minX = arg->min_chirp_dA;
+    maxX = arg->max_chirp_dA;
+    break;
+  case DIM_CHIRP_dW:
+    minX = arg->min_chirp_dW;
+    maxX = arg->max_chirp_dW;
+    break;
+  case DIM_CHIRP_ddA:
+    minX = arg->min_chirp_ddA;
+    maxX = arg->max_chirp_ddA;
+    break;
+  }
+
+  switch(arg->y_dim){
+  case DIM_ESTIMATE_A:
+    minY = arg->min_est_A;
+    maxY = arg->max_est_A;
+    break;
+  case DIM_ESTIMATE_P:
+    minY = arg->min_est_P;
+    maxY = arg->max_est_P;
+    break;
+  case DIM_ESTIMATE_W:
+    minY = arg->min_est_W;
+    maxY = arg->max_est_W;
+    break;
+  case DIM_ESTIMATE_dA:
+    minY = arg->min_est_dA;
+    maxY = arg->max_est_dA;
+    break;
+  case DIM_ESTIMATE_dW:
+    minY = arg->min_est_dW;
+    maxY = arg->max_est_dW;
+    break;
+  case DIM_ESTIMATE_ddA:
+    minY = arg->min_est_ddA;
+    maxY = arg->max_est_ddA;
+    break;
+  case DIM_CHIRP_A:
+    minY = arg->min_chirp_A;
+    maxY = arg->max_chirp_A;
+    break;
+  case DIM_CHIRP_P:
+    minY = arg->min_chirp_P;
+    maxY = arg->max_chirp_P;
+    break;
+  case DIM_CHIRP_W:
+    minY = arg->min_chirp_W;
+    maxY = arg->max_chirp_W;
+    break;
+  case DIM_CHIRP_dA:
+    minY = arg->min_chirp_dA;
+    maxY = arg->max_chirp_dA;
+    break;
+  case DIM_CHIRP_dW:
+    minY = arg->min_chirp_dW;
+    maxY = arg->max_chirp_dW;
+    break;
+  case DIM_CHIRP_ddA:
+    minY = arg->min_chirp_ddA;
+    maxY = arg->max_chirp_ddA;
+    break;
+  }
+
+  x0s = rint((x_n-1)/(maxX-minX)*minX);
+  y0s = rint((y_n-1)/(maxY-minY)*minY);
+  x1s = x0s+x_n-1;
+  y1s = y0s+y_n-1;
+
+  xminori = rint((x_n-1)/(maxX-minX)*arg->x_minor);
+  yminori = rint((y_n-1)/(maxY-minY)*arg->y_minor);
+  xmajori = rint(xminori/arg->x_minor*arg->x_major);
+  ymajori = rint(yminori/arg->y_minor*arg->y_major);
+
+  if(xminori<1 || yminori<1 || xmajori<1 || ymajori<1){
+    fprintf(stderr,"Bad xmajor/xminor/ymajor/yminor value.\n");
+    exit(1);
+  }
+
+  if( rint(xmajori*(maxX-minX)/arg->x_major) != x_n-1){
+    fprintf(stderr,"Xmajor or Xminor results in non-integer pixel increment.\n");
+    exit(1);
+  }
+
+  if( rint(ymajori*(maxY-minY)/arg->y_major) != y_n-1){
+    fprintf(stderr,"Ymajor or Yminor results in non-integer pixel increment.\n");
+    exit(1);
+  }
+
   /* determine ~ padding needed */
-  setup_graphs(arg->x0,arg->x1,arg->xmajor,
-               arg->y0,arg->y1,arg->ymajor,
+  setup_graphs(x0s,x1s,xmajori,arg->x_major,
+               y0s,y1s,ymajori,arg->y_major,
                (arg->subtitle1!=0)+(arg->subtitle2!=0)+(arg->subtitle3!=0),
                arg->fontsize);
 
-  /* don't check W; it's always swept in this graph */
-  if(arg->min_est_A != arg->max_est_A ||
-     arg->min_est_P != arg->max_est_P ||
-     arg->min_est_dA != arg->max_est_dA ||
-     arg->min_est_dW != arg->max_est_dW ||
+  if(!(arg->x_dim==DIM_ESTIMATE_A || arg->y_dim==DIM_ESTIMATE_A) &&
+     arg->min_est_A != arg->max_est_A) est_swept=1;
+  if(!(arg->x_dim==DIM_ESTIMATE_P || arg->y_dim==DIM_ESTIMATE_P) &&
+     arg->min_est_P != arg->max_est_P) est_swept=1;
+  if(!(arg->x_dim==DIM_ESTIMATE_W || arg->y_dim==DIM_ESTIMATE_W) &&
+     arg->min_est_W != arg->max_est_W) est_swept=1;
+  if(!(arg->x_dim==DIM_ESTIMATE_dA || arg->y_dim==DIM_ESTIMATE_dA) &&
+     arg->min_est_dA != arg->max_est_dA) est_swept=1;
+  if(!(arg->x_dim==DIM_ESTIMATE_dW || arg->y_dim==DIM_ESTIMATE_dW) &&
+     arg->min_est_dW != arg->max_est_dW) est_swept=1;
+  if(!(arg->x_dim==DIM_ESTIMATE_ddA || arg->y_dim==DIM_ESTIMATE_ddA) &&
      arg->min_est_ddA != arg->max_est_ddA) est_swept=1;
 
-  if(arg->min_chirp_A != arg->max_chirp_A ||
-     arg->min_chirp_P != arg->max_chirp_P ||
-     arg->min_chirp_dA != arg->max_chirp_dA ||
-     arg->min_chirp_dW != arg->max_chirp_dW ||
+  if(!(arg->x_dim==DIM_CHIRP_A || arg->y_dim==DIM_CHIRP_A) &&
+     arg->min_chirp_A != arg->max_chirp_A) chirp_swept=1;
+  if(!(arg->x_dim==DIM_CHIRP_P || arg->y_dim==DIM_CHIRP_P) &&
+     arg->min_chirp_P != arg->max_chirp_P) chirp_swept=1;
+  if(!(arg->x_dim==DIM_CHIRP_W || arg->y_dim==DIM_CHIRP_W) &&
+     arg->min_chirp_W != arg->max_chirp_W) chirp_swept=1;
+  if(!(arg->x_dim==DIM_CHIRP_dA || arg->y_dim==DIM_CHIRP_dA) &&
+     arg->min_chirp_dA != arg->max_chirp_dA) chirp_swept=1;
+  if(!(arg->x_dim==DIM_CHIRP_dW || arg->y_dim==DIM_CHIRP_dW) &&
+     arg->min_chirp_dW != arg->max_chirp_dW) chirp_swept=1;
+  if(!(arg->x_dim==DIM_CHIRP_ddA || arg->y_dim==DIM_CHIRP_ddA) &&
      arg->min_chirp_ddA != arg->max_chirp_ddA) chirp_swept=1;
+
   swept = est_swept | chirp_swept;
 
-  cC = draw_page(!swept?"Convergence":"Worst Case Convergence",
-                 arg->subtitle1,
-                 arg->subtitle2,
-                 arg->subtitle3,
-                 xaxis_label,
-                 yaxis_label,
-                 "Iterations:",
-                 DT_iterations);
-  if(swept)
+  if(arg->graph_convergence_max)
+    cC = draw_page(!swept?"Convergence":"Worst Case Convergence",
+                   arg->subtitle1,
+                   arg->subtitle2,
+                   arg->subtitle3,
+                   arg->xaxis_label,
+                   arg->yaxis_label,
+                   "Iterations:",
+                   DT_iterations);
+  if(swept && arg->graph_convergence_delta)
     cC_d = draw_page("Convergence Delta",
                      arg->subtitle1,
                      arg->subtitle2,
                      arg->subtitle3,
-                     xaxis_label,
-                     yaxis_label,
+                     arg->xaxis_label,
+                     arg->yaxis_label,
                      "Iteration span:",
                      DT_iterations);
 
-
-  cA = draw_page(!swept?"A (Amplitude) Error":"Maximum A (Amplitude) Error",
-                 arg->subtitle1,
-                 arg->subtitle2,
-                 arg->subtitle3,
-                 xaxis_label,
-                 yaxis_label,
-                 "Percentage Error:",
-                 DT_percent);
-  if(swept)
+  if(arg->graph_Aerror_max)
+    cA = draw_page(!swept?"A (Amplitude) Error":"Maximum A (Amplitude) Error",
+                   arg->subtitle1,
+                   arg->subtitle2,
+                   arg->subtitle3,
+                   arg->xaxis_label,
+                   arg->yaxis_label,
+                   "Percentage Error:",
+                   DT_percent);
+  if(swept && arg->graph_Aerror_delta)
     cA_d = draw_page("A (Amplitude) Delta",
                      arg->subtitle1,
                      arg->subtitle2,
                      arg->subtitle3,
-                     xaxis_label,
-                     yaxis_label,
+                     arg->xaxis_label,
+                     arg->yaxis_label,
                      "Delta:",
                      DT_abserror);
 
-
-  cP = draw_page(!swept?"P (Phase) Error":"Maximum P (Phase) Error",
-                 arg->subtitle1,
-                 arg->subtitle2,
-                 arg->subtitle3,
-                 xaxis_label,
-                 yaxis_label,
-                 "Error (rads/block):",
-                 DT_percent);
-  if(swept)
+  if(arg->graph_Perror_max)
+    cP = draw_page(!swept?"P (Phase) Error":"Maximum P (Phase) Error",
+                   arg->subtitle1,
+                   arg->subtitle2,
+                   arg->subtitle3,
+                   arg->xaxis_label,
+                   arg->yaxis_label,
+                   "Error (radians):",
+                   DT_abserror);
+  if(swept && arg->graph_Perror_delta)
     cP_d = draw_page("Phase Delta",
                      arg->subtitle1,
                      arg->subtitle2,
                      arg->subtitle3,
-                     xaxis_label,
-                     yaxis_label,
-                     "Delta (rads/block):",
+                     arg->xaxis_label,
+                     arg->yaxis_label,
+                     "Delta (radians):",
                      DT_abserror);
 
   if(arg->fit_W){
-    cW = draw_page(!swept?"W (Frequency) Error":"Maximum W (Frequency) Error",
-                   arg->subtitle1,
-                   arg->subtitle2,
-                   arg->subtitle3,
-                   xaxis_label,
-                   yaxis_label,
-                   "Error (rads/block):",
-                   DT_abserror);
-    if(swept)
+    if(arg->graph_Werror_max)
+      cW = draw_page(!swept?"W (Frequency) Error":"Maximum W (Frequency) Error",
+                     arg->subtitle1,
+                     arg->subtitle2,
+                     arg->subtitle3,
+                     arg->xaxis_label,
+                     arg->yaxis_label,
+                     "Error (cycles/block):",
+                     DT_abserror);
+    if(swept && arg->graph_Werror_delta)
       cW_d = draw_page("Frequency Delta",
                        arg->subtitle1,
                        arg->subtitle2,
                        arg->subtitle3,
-                       xaxis_label,
-                       yaxis_label,
-                       "Delta (rads/block):",
+                       arg->xaxis_label,
+                       arg->yaxis_label,
+                       "Delta (cycles/block):",
                        DT_abserror);
   }
 
   if(arg->fit_dA){
-    cdA = draw_page(!swept?"dA (Amplitude Modulation) Error":"Maximum dA (Amplitude Modulation) Error",
-                    arg->subtitle1,
-                    arg->subtitle2,
-                    arg->subtitle3,
-                    xaxis_label,
-                    yaxis_label,
-                    "Error:",
-                    DT_abserror);
-    if(swept)
+    if(arg->graph_dAerror_max)
+      cdA = draw_page(!swept?"dA (Amplitude Modulation) Error":"Maximum dA (Amplitude Modulation) Error",
+                      arg->subtitle1,
+                      arg->subtitle2,
+                      arg->subtitle3,
+                      arg->xaxis_label,
+                      arg->yaxis_label,
+                      "Error:",
+                      DT_abserror);
+    if(swept && arg->graph_dAerror_delta)
       cdA_d = draw_page("Amplitude Modulation Delta",
                         arg->subtitle1,
                         arg->subtitle2,
                         arg->subtitle3,
-                        xaxis_label,
-                        yaxis_label,
+                        arg->xaxis_label,
+                        arg->yaxis_label,
                         "Delta:",
                         DT_abserror);
   }
 
   if(arg->fit_dW){
-    cdW = draw_page(!swept?"dW (Chirp Rate) Error":"Maximum dW (Chirp Rate) Error",
-                   arg->subtitle1,
-                   arg->subtitle2,
-                   arg->subtitle3,
-                   xaxis_label,
-                   yaxis_label,
-                   "Error (rads/block):",
-                   DT_abserror);
-    if(swept)
+    if(arg->graph_dWerror_max)
+      cdW = draw_page(!swept?"dW (Chirp Rate) Error":"Maximum dW (Chirp Rate) Error",
+                      arg->subtitle1,
+                      arg->subtitle2,
+                      arg->subtitle3,
+                      arg->xaxis_label,
+                      arg->yaxis_label,
+                      "Error (cycles/block):",
+                      DT_abserror);
+    if(swept && arg->graph_dWerror_delta)
       cdW_d = draw_page("Chirp Rate Delta",
                         arg->subtitle1,
                         arg->subtitle2,
                         arg->subtitle3,
-                        xaxis_label,
-                        yaxis_label,
-                        "Delta (rads/block):",
+                        arg->xaxis_label,
+                        arg->yaxis_label,
+                        "Delta (cycles/block):",
                         DT_abserror);
   }
 
   if(arg->fit_ddA){
-    cddA = draw_page(!swept?"ddA (Amplitude Modulation Squared) Error":
-                     "Maximum ddA (Amplitude Modulation Squared) Error",
-                     arg->subtitle1,
-                     arg->subtitle2,
-                     arg->subtitle3,
-                     xaxis_label,
-                     yaxis_label,
-                     "Error:",
-                     DT_abserror);
-    if(swept)
+    if(arg->graph_ddAerror_max)
+      cddA = draw_page(!swept?"ddA (Amplitude Modulation Squared) Error":
+                       "Maximum ddA (Amplitude Modulation Squared) Error",
+                       arg->subtitle1,
+                       arg->subtitle2,
+                       arg->subtitle3,
+                       arg->xaxis_label,
+                       arg->yaxis_label,
+                       "Error:",
+                       DT_abserror);
+    if(swept && arg->graph_ddAerror_delta)
       cddA_d = draw_page("Amplitude Modulation Squared Delta",
                          arg->subtitle1,
                          arg->subtitle2,
                          arg->subtitle3,
-                         xaxis_label,
-                         yaxis_label,
+                         arg->xaxis_label,
+                         arg->yaxis_label,
                          "Delta:",
                          DT_abserror);
   }
 
-  cRMS = draw_page(swept?"Maximum RMS Fit Error":
-                   "RMS Fit Error",
-                   arg->subtitle1,
-                   arg->subtitle2,
-                   arg->subtitle3,
-                   xaxis_label,
-                   yaxis_label,
-                   "Percentage Error:",
-                   DT_percent);
-  if(swept)
+  if(arg->graph_RMSerror_max)
+    cRMS = draw_page(swept?"Maximum RMS Fit Error":
+                     "RMS Fit Error",
+                     arg->subtitle1,
+                     arg->subtitle2,
+                     arg->subtitle3,
+                     arg->xaxis_label,
+                     arg->yaxis_label,
+                     "Percentage Error:",
+                     DT_percent);
+  if(swept && arg->graph_RMSerror_delta)
     cRMS_d = draw_page("RMS Error Delta",
                        arg->subtitle1,
                        arg->subtitle2,
                        arg->subtitle3,
-                       xaxis_label,
-                       yaxis_label,
+                       arg->xaxis_label,
+                       arg->yaxis_label,
                        "Percentage Delta:",
                        DT_percent);
 
@@ -505,15 +767,13 @@
 
   /* graph computation */
   for(xi=0;xi<x_n;xi++){
-    int x = xi+arg->x0;
-    float w = ((float)x/arg->xmajor)/blocksize*2.*M_PI;
     chirp chirps[y_n];
     chirp estimates[y_n];
     int iter[y_n];
     float rms[y_n];
     int si,sn=(swept && arg->sweep_steps>1 ? arg->sweep_steps : 1);
 
-    fprintf(stderr,"\rW estimate distance vs. W graphs: column %d/%d...",x-arg->x0,x_n-1);
+    fprintf(stderr,"\r%s: column %d/%d...",filebase,xi,x_n-1);
 
     memset(targ,0,sizeof(targ));
 
@@ -547,24 +807,55 @@
 
       /* compute/set chirp and est parameters, potentially compute the
          chirp waveform */
-      for(i=0;i<y_n;i++){
-        int y = arg->y1-i;
-        float we=((float)y/arg->ymajor)/blocksize*2.*M_PI+w;
+      for(yi=0;yi<y_n;yi++){
 
-        set_chirp(chirps+i,arg->sweep_or_rand_p,si,sn,
-                  arg->min_chirp_A,arg->max_chirp_A,
-                  arg->min_chirp_P,arg->max_chirp_P,
-                  w,w,
-                  arg->min_chirp_dA,arg->min_chirp_dA,
-                  arg->min_chirp_dW,arg->min_chirp_dW,
-                  arg->min_chirp_ddA,arg->max_chirp_ddA);
-        set_chirp(estimates+i,arg->sweep_or_rand_p,si,sn,
-                  arg->min_est_A,arg->max_est_A,
-                  arg->min_est_P,arg->max_est_P,
-                  we,we,
-                  arg->min_est_dA,arg->min_est_dA,
-                  arg->min_est_dW,arg->min_est_dW,
-                  arg->min_est_ddA,arg->max_est_ddA);
+        set_chirp(chirps+yi,
+                  xi,x_n,
+                  arg->x_dim&DIM_CHIRP_MASK,
+                  y_n-yi-1,y_n,
+                  arg->y_dim&DIM_CHIRP_MASK,
+                  si,sn,
+                  arg->sweep_or_rand_p,
+                  arg->min_chirp_A,
+                  arg->max_chirp_A,
+                  arg->min_chirp_P*2.*M_PI,
+                  arg->max_chirp_P*2.*M_PI,
+                  arg->min_chirp_W*2.*M_PI/blocksize,
+                  arg->max_chirp_W*2.*M_PI/blocksize,
+                  arg->min_chirp_dA,
+                  arg->max_chirp_dA,
+                  arg->min_chirp_dW*2.*M_PI/blocksize/blocksize,
+                  arg->max_chirp_dW*2.*M_PI/blocksize/blocksize,
+                  arg->min_chirp_ddA,
+                  arg->max_chirp_ddA);
+
+        set_chirp(estimates+yi,
+                  xi,x_n,
+                  arg->x_dim&DIM_ESTIMATE_MASK,
+                  y_n-yi-1,y_n,
+                  arg->y_dim&DIM_ESTIMATE_MASK,
+                  si,sn,
+                  arg->sweep_or_rand_p,
+                  arg->min_est_A,
+                  arg->max_est_A,
+                  arg->min_est_P*2.*M_PI,
+                  arg->max_est_P*2.*M_PI,
+                  arg->min_est_W*2.*M_PI/blocksize,
+                  arg->max_est_W*2.*M_PI/blocksize,
+                  arg->min_est_dA,
+                  arg->max_est_dA,
+                  arg->min_est_dW*2.*M_PI/blocksize/blocksize,
+                  arg->max_est_dW*2.*M_PI/blocksize/blocksize,
+                  arg->min_est_ddA,
+                  arg->max_est_ddA);
+
+        if(arg->rel_est_A) estimates[yi].A += chirps[yi].A;
+        if(arg->rel_est_P) estimates[yi].P += chirps[yi].P;
+        if(arg->rel_est_W) estimates[yi].W += chirps[yi].W;
+        if(arg->rel_est_dA) estimates[yi].dA += chirps[yi].dA;
+        if(arg->rel_est_dW) estimates[yi].dW += chirps[yi].dW;
+        if(arg->rel_est_ddA) estimates[yi].ddA += chirps[yi].ddA;
+
       }
 
       if(!chirp_swept){
@@ -634,42 +925,49 @@
 
     /* Column sweep complete; plot */
     for(yi=0;yi<y_n;yi++){
-      int y=arg->y1-yi;
+      int x=x0s+xi;
+      int y=y1s-yi;
       float a = 1.;
-      if(x%arg->xminor==0 || y%arg->yminor==0) a = .8;
-      if(x%arg->xmajor==0 || y%arg->ymajor==0) a = .3;
+      if(x%xminori==0 || y%yminori==0) a = .8;
+      if(x%xmajori==0 || y%ymajori==0) a = .3;
 
       /* Convergence graph */
-      set_iter_color(cC,ret_maxiter[yi],a);
-      cairo_rectangle(cC,xi+leftpad,yi+toppad,1,1);
-      cairo_fill(cC);
+      if(cC){
+        set_iter_color(cC,ret_maxiter[yi],a);
+        cairo_rectangle(cC,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cC);
+      }
 
       /* Convergence delta graph */
-      if(swept){
+      if(cC_d){
         set_iter_color(cC_d,ret_maxiter[yi]-ret_miniter[yi],a);
         cairo_rectangle(cC_d,xi+leftpad,yi+toppad,1,1);
         cairo_fill(cC_d);
       }
 
       /* A error graph */
-      set_error_color(cA,MAX(fabs(ret_maxA[yi]),fabs(ret_minA[yi])),a);
-      cairo_rectangle(cA,xi+leftpad,yi+toppad,1,1);
-      cairo_fill(cA);
+      if(cA){
+        set_error_color(cA,MAX(fabs(ret_maxA[yi]),fabs(ret_minA[yi])),a);
+        cairo_rectangle(cA,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cA);
+      }
 
       /* A delta graph */
-      if(swept){
+      if(cA_d){
         set_error_color(cA_d,ret_maxA[yi]-ret_minA[yi],a);
         cairo_rectangle(cA_d,xi+leftpad,yi+toppad,1,1);
         cairo_fill(cA_d);
       }
 
       /* P error graph */
-      set_error_color(cP,MAX(fabs(ret_maxP[yi]),fabs(ret_minP[yi])),a);
-      cairo_rectangle(cP,xi+leftpad,yi+toppad,1,1);
-      cairo_fill(cP);
+      if(cP){
+        set_error_color(cP,MAX(fabs(ret_maxP[yi]),fabs(ret_minP[yi])),a);
+        cairo_rectangle(cP,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cP);
+      }
 
       /* P delta graph */
-      if(swept){
+      if(cP_d){
         set_error_color(cP_d,circular_distance(ret_maxP[yi],ret_minP[yi]),a);
         cairo_rectangle(cP_d,xi+leftpad,yi+toppad,1,1);
         cairo_fill(cP_d);
@@ -680,13 +978,13 @@
         set_error_color(cW,MAX(fabs(ret_maxW[yi]),fabs(ret_minW[yi]))/2./M_PI*blocksize,a);
         cairo_rectangle(cW,xi+leftpad,yi+toppad,1,1);
         cairo_fill(cW);
+      }
 
         /* W delta graph */
-        if(swept){
-          set_error_color(cW_d,(ret_maxW[yi]-ret_minW[yi])/2./M_PI*blocksize,a);
-          cairo_rectangle(cW_d,xi+leftpad,yi+toppad,1,1);
-          cairo_fill(cW_d);
-        }
+      if(cW_d){
+        set_error_color(cW_d,(ret_maxW[yi]-ret_minW[yi])/2./M_PI*blocksize,a);
+        cairo_rectangle(cW_d,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cW_d);
       }
 
       if(cdA){
@@ -694,27 +992,26 @@
         set_error_color(cdA,MAX(fabs(ret_maxdA[yi]),fabs(ret_mindA[yi]))*blocksize,a);
         cairo_rectangle(cdA,xi+leftpad,yi+toppad,1,1);
         cairo_fill(cdA);
-
+      }
         /* dA delta graph */
-        if(swept){
-          set_error_color(cdA_d,(ret_maxdA[yi]-ret_mindA[yi])*blocksize,a);
-          cairo_rectangle(cdA_d,xi+leftpad,yi+toppad,1,1);
-          cairo_fill(cdA_d);
-        }
+      if(cdA_d){
+        set_error_color(cdA_d,(ret_maxdA[yi]-ret_mindA[yi])*blocksize,a);
+        cairo_rectangle(cdA_d,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cdA_d);
       }
 
       if(cdW){
         /* dW error graph */
-        set_error_color(cdW,MAX(fabs(ret_maxdW[yi]),fabs(ret_mindW[yi]))/M_PI*blocksize*blocksize,a);
+        set_error_color(cdW,MAX(fabs(ret_maxdW[yi]),fabs(ret_mindW[yi]))/2./M_PI*blocksize*blocksize,a);
         cairo_rectangle(cdW,xi+leftpad,yi+toppad,1,1);
         cairo_fill(cdW);
+      }
 
-        /* dW delta graph */
-        if(swept){
-          set_error_color(cdW_d,(ret_maxdW[yi]-ret_mindW[yi])/M_PI*blocksize*blocksize,a);
-          cairo_rectangle(cdW_d,xi+leftpad,yi+toppad,1,1);
-          cairo_fill(cdW_d);
-        }
+      /* dW delta graph */
+      if(cdW_d){
+        set_error_color(cdW_d,(ret_maxdW[yi]-ret_mindW[yi])/2./M_PI*blocksize*blocksize,a);
+        cairo_rectangle(cdW_d,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cdW_d);
       }
 
       if(cddA){
@@ -722,22 +1019,24 @@
         set_error_color(cddA,MAX(fabs(ret_maxddA[yi]),fabs(ret_minddA[yi]))*blocksize*blocksize,a);
         cairo_rectangle(cddA,xi+leftpad,yi+toppad,1,1);
         cairo_fill(cddA);
+      }
 
-        /* dA delta graph */
-        if(swept){
-          set_error_color(cddA_d,(ret_maxddA[yi]-ret_minddA[yi])*blocksize*blocksize,a);
-          cairo_rectangle(cddA_d,xi+leftpad,yi+toppad,1,1);
-          cairo_fill(cddA_d);
-        }
+      /* dA delta graph */
+      if(cddA_d){
+        set_error_color(cddA_d,(ret_maxddA[yi]-ret_minddA[yi])*blocksize*blocksize,a);
+        cairo_rectangle(cddA_d,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cddA_d);
       }
 
       /* RMS error graph */
-      set_error_color(cRMS,ret_maxRMS[yi],a);
-      cairo_rectangle(cRMS,xi+leftpad,yi+toppad,1,1);
-      cairo_fill(cRMS);
+      if(cRMS){
+        set_error_color(cRMS,ret_maxRMS[yi],a);
+        cairo_rectangle(cRMS,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cRMS);
+      }
 
       /* RMS delta graph */
-      if(swept){
+      if(cRMS_d){
         set_error_color(cRMS_d,ret_maxRMS[yi]-ret_minRMS[yi],a);
         cairo_rectangle(cRMS_d,xi+leftpad,yi+toppad,1,1);
         cairo_fill(cRMS_d);
@@ -796,23 +1095,16 @@
 int main(){
   graph_run arg={
     /* fontsize */      12,
-    /* subtitle1 */     "graphtest1",
-    /* subtitle2 */     "graphtest2",
-    /* subtitle3 */     "graphtest3",
+    /* subtitle1 */     "Nonlinear (W, dW recentered), G/S, symmetric norm, no ddA fit",
+    /* subtitle2 */     "chirp: A=1.0, swept phase, dA=0 | estimate A=P=dA=dW=0, estimate W = chirp W",
+    /* subtitle3 */     "hanning window",
+    /* xaxis label */   "W (cycles/block)",
+    /* yaxis label */   "dW (cycles/block)",
+
     /* blocksize */     256,
     /* threads */       8,
 
-    /* x0 */            0,
-    /* x1 */            1280,
-    /* xmajor */        10,
-    /* xminor */        5,
-
-    /* y0 */            -300,
-    /* y1 */            300,
-    /* ymajor */        30,
-    /* yminor */        15,
-
-    /* window */        window_functions.maxwell1,
+    /* window */        window_functions.hanning,
     /* fit_tol */       .000001,
     /* gauss_seidel */  1,
     /* fit_W */         1,
@@ -825,25 +1117,61 @@
     /* symm_norm */     1,
     /* bound_zero */    0,
 
+    /* x dimension */   DIM_CHIRP_W,
+    /* x steps */       801,
+    /* x major */       1.,
+    /* x minor */       .25,
+    /* y dimension */   DIM_CHIRP_dW,
+    /* y steps */       601,
+    /* y major */       1.,
+    /* y minor */       .5,
     /* sweep_steps */   8,
     /* randomize_p */   0,
 
-    /* est A range */   0.,0.,
-    /* est P range */   0.,0.,
-    /* est W range */   0.,0.,
-    /* est dA range */  0.,0.,
-    /* est dW range */  0.,0.,
-    /* est ddA range */ 0.,0.,
+    /* est A range */     0.,  0.,  0, /* relative flag */
+    /* est P range */     0.,  0.,  0, /* relative flag */
+    /* est W range */     0.,  0.,  1, /* relative flag */
+    /* est dA range */    0.,  0.,  0, /* relative flag */
+    /* est dW range */    0.,  0.,  0, /* relative flag */
+    /* est ddA range */   0.,  0.,  0, /* relative flag */
 
     /* ch A range */    1.,1.,
-    /* ch P range */    0,2*M_PI,
-    /* ch W range */    0.,0.,
+    /* ch P range */    0,1.-1./8,
+    /* ch W range */    0.,8.,
     /* ch dA range */   0.,0.,
-    /* ch dW range */   0.,0.,
+    /* ch dW range */   -6.,6.,
     /* ch ddA range */  0.,0.,
+
+    /* converge max */    1,
+    /* converge del */    1,
+    /* max A error */     1,
+    /* A error delta */   0,
+    /* max P error */     1,
+    /* P error delta */   0,
+    /* max W error */     1,
+    /* W error delta */   0,
+    /* max dA error */    1,
+    /* dA error delta */  0,
+    /* max dW error */    1,
+    /* dW error delta */  0,
+    /* max ddA error */   1,
+    /* ddA error delta */ 0,
+    /* max RMS error */   1,
+    /* RMS error delta */ 0,
+
   };
 
-  w_e("graphs-A",&arg);
+  w_e("dW-vs-W",&arg);
+
+  arg.yaxis_label="initial distance from W (cycles/block)";
+  arg.y_dim = DIM_ESTIMATE_W;
+  arg.min_est_W = -6.;
+  arg.max_est_W =  6.;
+  arg.min_chirp_dW=0.;
+  arg.max_chirp_dW=0.;
+
+  w_e("estW-vs-W",&arg);
+
   return 0;
 }
 



More information about the commits mailing list