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

xiphmont at svn.xiph.org xiphmont at svn.xiph.org
Sun Apr 17 22:39:44 PDT 2011


Author: xiphmont
Date: 2011-04-17 22:39:44 -0700 (Sun, 17 Apr 2011)
New Revision: 17921

Added:
   trunk/ghost/monty/chirp/chirpgraph.h
   trunk/ghost/monty/chirp/chirptest.c
Modified:
   trunk/ghost/monty/chirp/chirp.c
   trunk/ghost/monty/chirp/chirp.h
   trunk/ghost/monty/chirp/chirpgraph.c
   trunk/ghost/monty/chirp/window.c
Log:
Get some of the new graphing code in SVN; the various code paths are not yet well tested.



Modified: trunk/ghost/monty/chirp/chirp.c
===================================================================
--- trunk/ghost/monty/chirp/chirp.c	2011-04-12 06:45:19 UTC (rev 17920)
+++ trunk/ghost/monty/chirp/chirp.c	2011-04-18 05:39:44 UTC (rev 17921)
@@ -113,16 +113,21 @@
 
                              float fit_limit,
                              int iter_limit,
+                             int fit_gs,
 
                              int fitW,
                              int fitdA,
                              int fitdW,
                              int fitddA,
+                             int fit_recenter_dW,
+                             float fit_W_alpha,
+                             float fit_dW_alpha,
+
                              int symm_norm,
                              float E0,
                              float E1,
                              float E2,
-                             int fit_compound,
+
                              int bound_edges){
   int i,j;
   int flag=1;
@@ -130,8 +135,8 @@
 
   /* outer fit iteration */
   while(flag && iter_limit>0){
+    flag=0;
     iter_limit--;
-    flag=0;
 
     /* precompute the portion of the projection/fit estimate shared by
        the zero, first and second order fits.  Subtracts the current
@@ -157,13 +162,20 @@
       float aC = cos(c->P);
       float aS = sin(c->P);
 
+      /* Not recentering dW allows potential simplification of the
+         nonlinear solver. This code is designed to recenter dW as
+         easily as not, so the following looks a bit silly.  The point
+         of the flag is to emulate/study the behavior of the simplified
+         algorithm */
+      float cdW = fit_recenter_dW ? c->dW : 0;
+
       for(j=0;j<len;j++){
         float jj = j-len*.5+.5;
         float jj2 = jj*jj;
         float co,si,c2,s2;
         float yy=r[j];
 
-        sincosf((c->W + c->dW*jj)*jj,&si,&co);
+        sincosf((c->W + cdW*jj)*jj,&si,&co);
         si*=window[j];
         co*=window[j];
         c2 = co*co*jj;
@@ -182,7 +194,7 @@
         eP += co*yy*jj2;
         fP += si*yy*jj2;
 
-        if(fit_compound){
+        if(fit_gs){
           /* subtract zero order estimate from first */
           cP2 += c2;
           dP2 += s2;
@@ -227,6 +239,17 @@
       if(!fitdW && !fitddA)
         eP=fP=0;
 
+      /* base convergence on relative projection movement this
+         iteration */
+      {
+        float move = (aP*aP + bP*bP)/(c->A*c->A + fit_limit*fit_limit) +
+          (cP*cP + dP*dP)/(c->A*c->A + fit_limit*fit_limit) +
+          (eP*eP + fP*fP)/(c->A*c->A + fit_limit*fit_limit);
+
+        if(fit_limit>0 && move>fit_limit*fit_limit)flag=1;
+        if(fit_limit<0 && move>1e-14)flag=1;
+      }
+
       /* we're fitting to the remaining error; add the fit to date
          back in to relate our newest incremental results to the
          global fit so far.  Note that this does not include W or dW,
@@ -235,30 +258,13 @@
       {
         float A = toAi(c->A, c->P);
         float B = toBi(c->A, c->P);
-        float C = dAtoCi(A,B,c->dA) + cP;
-        float D = dAtoDi(A,B,c->dA) + dP;
-        float E = ddAtoEi(A,B,c->ddA) + eP;
-        float F = ddAtoFi(A,B,c->ddA) + fP;
-        A += aP;
-        B += bP;
+        aP += A;
+        bP += B;
+        cP += dAtoCi(A,B,c->dA);
+        dP += dAtoDi(A,B,c->dA);
+        eP += ddAtoEi(A,B,c->ddA);
+        fP += ddAtoFi(A,B,c->ddA);
 
-        /* base convergence on relative projection movement this
-           iteration */
-
-        float move = (aP*aP + bP*bP)/(A*A + B*B + fit_limit*fit_limit) +
-          (cP*cP + dP*dP)/(A*A + B*B + fit_limit*fit_limit) +
-          (eP*eP + fP*fP)/(A*A + B*B + fit_limit*fit_limit);
-
-        if(fit_limit>0 && move>fit_limit*fit_limit)flag=1;
-        if(fit_limit<0 && move>1e-14)flag=1;
-
-        aP = A;
-        bP = B;
-        cP = C;
-        dP = D;
-        eP = E;
-        fP = F;
-
         /* guard overflow; if we're this far out, assume we're never
            coming back. drop out now. */
         if((aP*aP + bP*bP)>1e20 ||
@@ -272,9 +278,12 @@
       /* save new estimate */
       c->A = toA(aP,bP);
       c->P = toP(aP,bP);
-      c->W += (fitW ? toW(aP,bP,cP,dP) : 0);
+      c->W += fit_W_alpha*(fitW ? toW(aP,bP,cP,dP) : 0);
       c->dA = (fitdA ? todA(aP,bP,cP,dP) : 0);
-      c->dW += 1.75*(fitdW ? todW(aP,bP,eP,fP) : 0);
+      if(fit_recenter_dW)
+        c->dW += fit_dW_alpha*(fitdW ? todW(aP,bP,eP,fP) : 0);
+      else
+        c->dW = (fitdW ? todW(aP,bP,eP,fP) : 0);
       c->ddA = (fitddA ? toddA(aP,bP,eP,fP) : 0);
 
       if(bound_edges){
@@ -282,7 +291,7 @@
         if(c->W<0){
           c->W = 0; /* clamp frequency to 0 (DC) */
           c->dW = 0; /* if W is 0, the chirp rate must also be 0 to
-                           avoid negative frequencies */
+                                 avoid negative frequencies */
         }
         if(c->W>M_PI){
           c->W = M_PI; /* clamp frequency to Nyquist */
@@ -299,22 +308,24 @@
           c->dW = c->W/len;
         /* ...or exceed Nyquist */
         if(c->W + c->dW*len > M_PI)
-          c->dW = (M_PI - c->W)/len;
+          c->dW = M_PI/len - c->W/len;
         if(c->W - c->dW*len > M_PI)
-          c->dW = (M_PI + c->W)/len;
+          c->dW = c->W/len - M_PI/len;
       }
 
       /* update the reconstruction/residue vectors with new fit */
-      for(j=0;j<len;j++){
-        float jj = j-len*.5+.5;
-        float a = c->A + (c->dA + c->ddA*jj)*jj;
-        float v = a*cosf((c->W + c->dW*jj)*jj + c->P);
-        r[j] -= v*window[j];
-        y[j] += v;
+      {
+        float cdW = fit_recenter_dW ? c->dW : 0;
+        for(j=0;j<len;j++){
+          float jj = j-len*.5+.5;
+          float a = c->A + (c->dA + c->ddA*jj)*jj;
+          float v = a*cosf((c->W + cdW*jj)*jj + c->P);
+          r[j] -= v*window[j];
+          y[j] += v;
+        }
       }
     }
   }
-
   return iter_limit;
 }
 
@@ -334,6 +345,7 @@
 
                           float fit_limit,
                           int iter_limit,
+                          int fit_gs,
 
                           int fitW,
                           int fitdA,
@@ -342,8 +354,7 @@
                           int symm_norm,
                           float E0,
                           float E1,
-                          float E2,
-                          int fit_compound){
+                          float E2){
 
   float *cos_table[n];
   float *sin_table[n];
@@ -410,7 +421,7 @@
       tmpe += ttcos_table[i][j]*ttcos_table[i][j]*window[j]*window[j];
       tmpf += ttsin_table[i][j]*ttsin_table[i][j]*window[j]*window[j];
 
-      /* compound fit terms */
+      /* gs fit terms */
       tmpc2 += cos_table[i][j]*tcos_table[i][j]*window[j]*window[j];
       tmpd2 += sin_table[i][j]*tsin_table[i][j]*window[j]*window[j];
       tmpe3 += tcos_table[i][j]*ttcos_table[i][j]*window[j]*window[j];
@@ -432,7 +443,7 @@
       ttsinE[i] = (tmpf>0.f ? 1./tmpf : 0);
     }
 
-    /* set compound fit terms */
+    /* set gs fit terms */
     tcosC2[i] = tmpc2;
     tsinC2[i] = tmpd2;
     ttcosC2[i] = tmpc;
@@ -455,8 +466,8 @@
   }
 
   while(flag && iter_limit){
+    flag=0;
     iter_limit--;
-    flag=0;
     for (i=0;i<n;i++){
 
       float tmpa=0, tmpb=0;
@@ -477,7 +488,7 @@
       tmpa*=cosE[i];
       tmpb*=sinE[i];
 
-      if(fit_compound){
+      if(fit_gs){
         tmpc -= tmpa*tcosC2[i];
         tmpd -= tmpb*tsinC2[i];
       }
@@ -485,7 +496,7 @@
       tmpc*=tcosE[i];
       tmpd*=tsinE[i];
 
-      if(fit_compound){
+      if(fit_gs){
         tmpe -= tmpa*ttcosC2[i] + tmpc*ttcosC3[i];
         tmpf -= tmpb*ttsinC2[i] + tmpd*ttsinC3[i];
       }
@@ -562,7 +573,7 @@
    fitdW : flag indicating that fitting should include the dW parameter.
    fitddA : flag indicating that fitting should include the ddA parameter.
    symm_norm
-   int fit_compound
+   int fit_gs
    int bound_zero
 
    Input estimates affect convergence region and speed and fit
@@ -584,16 +595,19 @@
                     int len,
                     chirp *c,
                     int n,
+
+                    float fit_limit,
                     int iter_limit,
-                    float fit_limit,
+                    int fit_gs,
 
-                    int linear,
                     int fitW,
                     int fitdA,
                     int fitdW,
                     int fitddA,
+                    int nonlinear,
+                    float fit_W_alpha,
+                    float fit_dW_alpha,
                     int symm_norm,
-                    int fit_compound,
                     int bound_zero
                     ){
 
@@ -633,19 +647,21 @@
     }
   }
 
-  if(linear){
+  if(!nonlinear){
     if(bound_zero) return -1;
+    if(fit_W_alpha!=1.0) return -1;
+    if(fit_dW_alpha!=1.0) return -1;
     iter_limit = linear_iterate(x,y,window,len,c,n,
-                                fit_limit,iter_limit,
+                                fit_limit,iter_limit,fit_gs,
                                 fitW,fitdA,fitdW,fitddA,
-                                symm_norm,E0,E1,E2,
-                                fit_compound);
+                                symm_norm,E0,E1,E2);
   }else{
     iter_limit = nonlinear_iterate(x,y,window,len,c,n,
-                                   fit_limit,iter_limit,
+                                   fit_limit,iter_limit,fit_gs,
                                    fitW,fitdA,fitdW,fitddA,
+                                   nonlinear-1,fit_W_alpha,fit_dW_alpha,
                                    symm_norm,E0,E1,E2,
-                                   fit_compound,bound_zero);
+                                   bound_zero);
   }
 
   /* Sort by ascending frequency */

Modified: trunk/ghost/monty/chirp/chirp.h
===================================================================
--- trunk/ghost/monty/chirp/chirp.h	2011-04-12 06:45:19 UTC (rev 17920)
+++ trunk/ghost/monty/chirp/chirp.h	2011-04-18 05:39:44 UTC (rev 17921)
@@ -25,23 +25,29 @@
   int label;/* used for tracking by outside code */
 } chirp;
 
-extern int estimate_chirps(const float *x,
-                           float *y,
-                           const float *window,
-                           int len,
-                           chirp *c,
-                           int n,
-                           int iter_limit,
-                           float fit_limit,
+extern int
+estimate_chirps(const float *x, /* unwindowed input to fit */
+                float *y,       /* unwindowed outptu reconstruction */
+                const float *window, /* window to apply to input/bases */
+                int len,        /* block length */
+                chirp *c,       /* list of chirp estimates/outputs */
+                int n,          /* number of chirps */
+                float fit_limit,/* minimum basis movement to continue iteration */
+                int iter_limit, /* maximum number of iterations */
+                int fit_gs,     /* Use Gauss-Seidel partial updates */
 
-                           int linear,
-                           int fitW,
-                           int fitdA,
-                           int fitdW,
-                           int fitddA,
-                           int symm_norm,
-                           int fit_compound,
-                           int bound_zero);
+                int fitW,       /* fit the W parameter */
+                int fitdA,      /* fit the dA parameter */
+                int fitdW,      /* fit the dW parameter */
+                int fitddA,     /* fit the ddA parameter */
+                int nonlinear,  /* perform a linear fit (0),
+                                   nonlinear fit recentering W only (1)
+                                   nonlinear fit recentering W and dW (2) */
+                float fit_W_alpha, /* W alpha multiplier for nonlinear fit */
+                float fit_dW_alpha,/* dW alpha multiplier for nonlinear fit */
+                int symm_norm,     /* Use symmetric normalization optimization */
+                int bound_zero);   /* prevent W or dW from fitting to negative or
+                                      greater-then-Nyquist frequencies */
 
 extern void advance_chirps(chirp *c, int n, int len);
 

Modified: trunk/ghost/monty/chirp/chirpgraph.c
===================================================================
--- trunk/ghost/monty/chirp/chirpgraph.c	2011-04-12 06:45:19 UTC (rev 17920)
+++ trunk/ghost/monty/chirp/chirpgraph.c	2011-04-18 05:39:44 UTC (rev 17921)
@@ -18,6 +18,7 @@
 #define _GNU_SOURCE
 #include <math.h>
 #include "chirp.h"
+#include "chirpgraph.h"
 #include "scales.h"
 #include <stdio.h>
 #include <stdlib.h>
@@ -25,61 +26,6 @@
 #include <cairo/cairo.h>
 #include <pthread.h>
 
-/* configuration globals. */
-int BLOCKSIZE=2048;
-int cores=8;
-int xstepd=100;
-int ystepd=100;
-int xstepg=25;
-int ystepg=25;
-int x_n=801;
-int y_n=300;
-int fontsize=12;
-
-int randomize_est_A=0;
-int randomize_est_P=0;
-int randomize_est_dA=0;
-int randomize_est_dW=0;
-int randomize_est_ddA=0;
-
-int randomize_chirp_P=0;
-int randomize_chirp_dA=0;
-int randomize_chirp_dW=0;
-int randomize_chirp_ddA=0;
-
-float initial_est_A = 0.;
-float initial_est_P = 0.;
-float initial_est_dA = 0.;
-float initial_est_dW = 0.;
-float initial_est_ddA = 0.;
-
-float initial_chirp_A = 1.;
-float initial_chirp_P = M_PI/4;
-float initial_chirp_dA = 0.;
-float initial_chirp_dW = 0.;
-float initial_chirp_ddA = 0.;
-
-int fit_linear=0;
-int fit_W=1;
-int fit_dA=1;
-int fit_dW=1;
-int fit_ddA=1;
-int fit_symm_norm=1;
-int fit_gauss_seidel=1;
-int fit_bound_zero=0;
-float fit_tolerance=.000001;
-int fit_hanning=1;
-int iterations=-1;
-
-void hanning(float *x, int n){
-  float scale = 2*M_PI/n;
-  int i;
-  for(i=0;i<n;i++){
-    float i5 = i+.5;
-    x[i] = .5-.5*cos(scale*i5);
-  }
-}
-
 /********************** colors for graphs *********************************/
 
 void set_iter_color(cairo_t *cC, int ret, float a){
@@ -91,7 +37,7 @@
     cairo_set_source_rgba(cC,1-(ret-30)*.05,0,0,a);
   }else if (ret>10){ /* 11(yellow) - 30(red) */
     cairo_set_source_rgba(cC,1,1-(ret-10)*.05,0,a);
-  }else if (ret>4){ /*  5 (green) - 10(yellow) */
+  }else if (ret>4){ /*  5 (green) - 10 (yellow) */
     cairo_set_source_rgba(cC,(ret-5)*.15+.25,1,0,a);
   }else if (ret==4){ /*  4 brighter cyan */
     cairo_set_source_rgba(cC,0,.8,.4,a);
@@ -130,121 +76,95 @@
 
 /********* draw everything in the graph except the graph data itself *******/
 
-#define DT_iterations 0
-#define DT_abserror 1
-#define DT_percent 2
+static float fontsize=12;
+static int x0s,x1s,xmajor;
+static int y0s,y1s,ymajor;
 
 /* computed configuration globals */
 int leftpad=0;
-int rightpad=0;
+static int rightpad=0;
 int toppad=0;
-int bottompad=0;
-float legendy=0;
-float legendh=0;
-int pic_w=0;
-int pic_h=0;
-float titleheight=0.;
+static int bottompad=0;
+static float legendy=0;
+static float legendh=0;
+static int pic_w=0;
+static int pic_h=0;
+static float titleheight=0.;
+static int x_n=0;
+static int y_n=0;
 
-/* determines padding, etc.  Call several times to determine maximums. */
-void setup_graph(char *title1,
-                 char *title2,
-                 char *title3,
-                 char *xaxis_label,
-                 char *yaxis_label,
-                 char *legend_label,
-                 int datatype){
+/* determines padding, etc.  Will not expand the frame to prevent
+   overruns of user-set titles/legends */
+void setup_graphs(int start_x_step,
+                  int end_x_step, /* inclusive; not one past */
+                  int x_major_d,
 
+                  int start_y_step,
+                  int end_y_step, /* inclusive; not one past */
+                  int y_major_d,
+
+                  int subtitles,
+                  float fs){
+
   /* determine ~ padding needed */
   cairo_surface_t *cs=cairo_image_surface_create(CAIRO_FORMAT_RGB24,100,100);
   cairo_t *ct=cairo_create(cs);
-  int minx=-leftpad,maxx=x_n+rightpad;
-  int miny=-toppad,maxy=2*y_n+1+bottompad;
-  int textpad=fontsize;
   cairo_text_extents_t extents;
-  int x,y;
+  cairo_font_extents_t fextents;
+  int y;
 
+  x_n = end_x_step-start_x_step+1;
+  y_n = end_y_step-start_y_step+1;
+  fontsize=fs;
+  x0s=start_x_step;
+  x1s=end_x_step;
+  xmajor=x_major_d;
+  y0s=start_y_step;
+  y1s=end_y_step;
+  ymajor=y_major_d;
+
   /* y labels */
   cairo_set_font_size(ct, fontsize);
-  for(y=0;y<=y_n;y+=ystepd){
-    char buf[80];
-    snprintf(buf,80,"%.0f",(float)y/ystepd);
-    cairo_text_extents(ct, buf, &extents);
-    if(-textpad-extents.height*.5<miny)miny=-textpad-extents.height*.5;
-    if(-textpad-extents.width-extents.height*1.5<minx)minx=-textpad-extents.width-extents.height*1.5;
-
-    if(y>0 && y<y_n){
-      snprintf(buf,80,"%.0f",(float)-y/ystepd);
+  for(y=y0s+fontsize;y<=y1s;y++){
+    if(y%ymajor==0){
+      char buf[80];
+      snprintf(buf,80,"%.0f",(float)y/ymajor);
       cairo_text_extents(ct, buf, &extents);
-      if(-textpad-extents.height*.5<miny)miny=-textpad-extents.height*.5;
-      if(-textpad-extents.width<minx)minx=-textpad-extents.width;
+      if(extents.width + extents.height*3.5>leftpad)leftpad=extents.width + extents.height*3.5;
     }
   }
 
   /* x labels */
-  for(x=0;x<=x_n;x+=xstepd){
-    char buf[80];
-    snprintf(buf,80,"%.0f",(float)x/xstepd);
-    cairo_text_extents(ct, buf, &extents);
+  cairo_font_extents(ct, &fextents);
+  if(fextents.height*3>bottompad)bottompad=fextents.height*3;
 
-    if(y_n*2+1+extents.height*3>maxy)maxy=y_n*2+1+extents.height*3;
-    if(x-(extents.width/2)-textpad<minx)minx=x-(extents.width/2)-textpad;
-    if(x+(extents.width/2)+textpad>maxx)maxx=x+(extents.width/2)+textpad;
-  }
-
   /* center horizontally */
-  if(maxx-x_n < -minx)maxx = x_n-minx;
-  if(maxx-x_n > -minx)minx = x_n-maxx;
+  if(leftpad<rightpad)leftpad=rightpad;
+  if(rightpad<leftpad)rightpad=leftpad;
 
   /* top legend */
   {
     float sofar=0;
-    if(title1){
-      cairo_save(ct);
-      cairo_select_font_face (ct, "",
-                              CAIRO_FONT_SLANT_NORMAL,
-                              CAIRO_FONT_WEIGHT_BOLD);
-      cairo_set_font_size(ct,fontsize*1.25);
-      cairo_text_extents(ct, title1, &extents);
+    cairo_save(ct);
+    cairo_select_font_face (ct, "",
+                            CAIRO_FONT_SLANT_NORMAL,
+                            CAIRO_FONT_WEIGHT_BOLD);
+    cairo_set_font_size(ct,fontsize*1.25);
+    cairo_font_extents(ct, &fextents);
+    sofar+=fextents.height;
+    cairo_restore(ct);
 
-      if(extents.width+textpad > maxx-minx){
-        int moar = extents.width+textpad-maxx+minx;
-        minx-=moar/2;
-        maxx+=moar/2;
-      }
-      sofar+=extents.height;
-      cairo_restore(ct);
-    }
-    if(title2){
-      cairo_text_extents(ct, title2, &extents);
-
-      if(extents.width+textpad > maxx-minx){
-        int moar = extents.width+textpad-maxx+minx;
-        minx-=moar/2;
-        maxx+=moar/2;
-      }
-      sofar+=extents.height;
-    }
-    if(title3){
-      cairo_text_extents(ct, title3, &extents);
-
-      if(extents.width+textpad > maxx-minx){
-        int moar = extents.width+textpad-maxx+minx;
-        minx-=moar/2;
-        maxx+=moar/2;
-      }
-      sofar+=extents.height;
-    }
+    cairo_font_extents(ct, &fextents);
+    while(subtitles--)
+      sofar+=fextents.height;
     if(sofar>titleheight)titleheight=sofar;
-    if(-miny<titleheight+textpad*2)miny=-titleheight-textpad*2;
-
+    if(toppad<titleheight+fontsize*2)toppad=titleheight+fontsize*2;
   }
 
   /* color legend */
   {
     float width=0.;
     float w1,w2,w3;
-    cairo_text_extents(ct, legend_label, &extents);
-    width=extents.width;
 
     cairo_text_extents(ct, "100", &extents);
     w1=extents.width*2*11;
@@ -255,89 +175,79 @@
 
     width+=MAX(w1,MAX(w2,w3));
 
-    legendy = y_n*2+1+extents.height*4;
+    legendy = y_n+extents.height*4;
     legendh = extents.height*1.5;
 
-    if(width + textpad > x_n-minx){
-      int moar = width+textpad-x_n;
-      minx-=moar;
-      maxx+=moar;
-    }
-    if(legendy+legendh*4>maxy)
-      maxy=legendy+legendh*2.5;
+    if(legendy+legendh*4>bottompad)
+      bottompad=(legendy-y_n)+legendh*2.5;
   }
 
-  toppad = -miny;
-  leftpad = -minx;
-  rightpad = maxx-x_n;
-  bottompad = maxy-(y_n*2+1);
-  if(toppad<leftpad*.75)toppad = leftpad*.75;
-  if(leftpad<toppad)leftpad=toppad;
-  if(rightpad<toppad)rightpad=toppad;
-
-  pic_w = maxx-minx;
-  pic_h = maxy-miny;
+  pic_w = x_n + leftpad + rightpad;
+  pic_h = y_n + toppad + bottompad;
 }
 
 /* draws the page surrounding the graph data itself */
-void draw_page(cairo_t *c,
-               char *title1,
-               char *title2,
-               char *title3,
-               char *xaxis_label,
-               char *yaxis_label,
-               char *legend_label,
-               int datatype){
+cairo_t *draw_page(char *title,
+                   char *subtitle1,
+                   char *subtitle2,
+                   char *subtitle3,
+                   char *xaxis_label,
+                   char *yaxis_label,
+                   char *legend_label,
+                   int datatype){
 
   int i;
   cairo_text_extents_t extents;
+  cairo_font_extents_t fextents;
+  cairo_surface_t *cs=cairo_image_surface_create(CAIRO_FORMAT_ARGB32,pic_w, pic_h);
+  if(!cs || cairo_surface_status(cs)!=CAIRO_STATUS_SUCCESS){
+    fprintf(stderr,"Could not set up Cairo surface.\n\n");
+    exit(1);
+  }
+  cairo_t *c = cairo_create(cs);
   cairo_save(c);
 
   /* clear page to white */
   cairo_set_source_rgb(c,1,1,1);
   cairo_rectangle(c,0,0,pic_w,pic_h);
   cairo_fill(c);
+  cairo_set_font_size(c, fontsize);
 
   /* set graph area to transparent */
   cairo_set_source_rgba(c,0,0,0,0);
   cairo_set_operator(c,CAIRO_OPERATOR_SOURCE);
-  cairo_rectangle(c,leftpad,toppad,x_n,y_n*2+1);
+  cairo_rectangle(c,leftpad,toppad,x_n,y_n);
   cairo_fill(c);
   cairo_restore(c);
 
   /* Y axis numeric labels */
   cairo_set_source_rgb(c,0,0,0);
-  for(i=0;i<=y_n;i+=ystepd){
-    int y = toppad+y_n;
-    char buf[80];
+  for(i=y0s+fontsize;i<=y1s;i++){
+    if(i%ymajor==0){
+      int y = toppad+y_n-i+y0s;
+      char buf[80];
 
-    snprintf(buf,80,"%.0f",(float)i/ystepd);
-    cairo_text_extents(c, buf, &extents);
-    cairo_move_to(c,leftpad - fontsize*.5 - extents.width,y-i+extents.height*.5);
-    cairo_show_text(c,buf);
-
-    if(i>0 && i<y_n){
-
-      snprintf(buf,80,"%.0f",(float)-i/ystepd);
+      snprintf(buf,80,"%.0f",(float)i/ymajor);
       cairo_text_extents(c, buf, &extents);
-      cairo_move_to(c,leftpad - fontsize*.5 - extents.width,y+i+extents.height*.5);
+      cairo_move_to(c,leftpad - fontsize*.5 - extents.width,y+extents.height*.5);
       cairo_show_text(c,buf);
-
     }
   }
 
   /* X axis labels */
-  for(i=0;i<=x_n;i+=xstepd){
-    char buf[80];
+  for(i=x0s;i<=x1s;i++){
+    if(i%xmajor==0){
+      char buf[80];
 
-    if(i==0){
-      snprintf(buf,80,"DC");
-    }else{
-      snprintf(buf,80,"%.0f",(float)i/xstepd);
+      if(i==0){
+        snprintf(buf,80,"DC");
+      }else{
+        snprintf(buf,80,"%.0f",(float)i/xmajor);
+      }
+      cairo_text_extents(c, buf, &extents);
+      cairo_move_to(c,leftpad + i - extents.width/2,y_n+toppad+extents.height+fontsize*.5);
+      cairo_show_text(c,buf);
     }
-    cairo_text_extents(c, buf, &extents);
-    cairo_move_to(c,leftpad + i - extents.width/2,y_n*2+1+toppad+extents.height+fontsize*.5);
-    cairo_show_text(c,buf);
   }
 
   /* Y axis label */
@@ -346,7 +256,7 @@
     cairo_matrix_t b = {0.,-1., 1.,0., 0.,0.}; // account for border!
     cairo_matrix_t d;
     cairo_text_extents(c, yaxis_label, &extents);
-    cairo_move_to(c,leftpad-extents.height*2,y_n+toppad+extents.width*.5);
+    cairo_move_to(c,extents.height+fontsize,y_n/2+toppad+extents.width*.5);
 
     cairo_save(c);
     cairo_get_matrix(c,&a);
@@ -359,37 +269,45 @@
   /* X axis caption */
   {
     cairo_text_extents(c, xaxis_label, &extents);
-    cairo_move_to(c,pic_w/2-extents.width/2,y_n*2+1+toppad+extents.height*2+fontsize*.5);
+    cairo_move_to(c,pic_w/2-extents.width/2,y_n+toppad+extents.height*2+fontsize*.5);
     cairo_show_text(c,xaxis_label);
   }
 
   /* top title(s) */
   {
     float y = (toppad-titleheight);
-    if(title1){
+    if(title){
       cairo_save(c);
       cairo_select_font_face (c, "",
                               CAIRO_FONT_SLANT_NORMAL,
                               CAIRO_FONT_WEIGHT_BOLD);
       cairo_set_font_size(c,fontsize*1.25);
-      cairo_text_extents(c, title1, &extents);
+      cairo_font_extents(c, &fextents);
+      cairo_text_extents(c, title, &extents);
       cairo_move_to(c,pic_w/2-extents.width/2,y);
-      cairo_show_text(c,title1);
-      y+=extents.height;
+      cairo_show_text(c,title);
+      y+=fextents.height;
       cairo_restore(c);
     }
-    if(title2){
-      cairo_text_extents(c, title2, &extents);
+    cairo_font_extents(c, &fextents);
+    if(subtitle1){
+      cairo_text_extents(c, subtitle1, &extents);
       cairo_move_to(c,pic_w/2-extents.width/2,y);
-      cairo_show_text(c,title2);
-      y+=extents.height;
+      cairo_show_text(c,subtitle1);
+      y+=fextents.height;
     }
-    if(title3){
-      cairo_text_extents(c, title3, &extents);
+    if(subtitle2){
+      cairo_text_extents(c, subtitle2, &extents);
       cairo_move_to(c,pic_w/2-extents.width/2,y);
-      cairo_show_text(c,title3);
-      y+=extents.height;
+      cairo_show_text(c,subtitle2);
+      y+=fextents.height;
     }
+    if(subtitle3){
+      cairo_text_extents(c, subtitle3, &extents);
+      cairo_move_to(c,pic_w/2-extents.width/2,y);
+      cairo_show_text(c,subtitle3);
+      y+=fextents.height;
+    }
   }
 
   /* color legend */
@@ -500,331 +418,23 @@
       cairo_show_text(c,legend_label);
     }
   }
+
+  return c;
 }
 
-void to_png(cairo_surface_t *c,char *base, char *name){
+void to_png(cairo_t *c,char *base, char *name){
   if(c){
+    cairo_surface_t *cs=cairo_get_target(c);
     char buf[320];
-    snprintf(buf,320,"%s-%s.png",name,base);
-    cairo_surface_write_to_png(c,buf);
+    snprintf(buf,320,"%s-%s.png",base,name);
+    cairo_surface_write_to_png(cs,buf);
   }
 }
 
-/*********************** Plot w estimate error against w **********************/
-
-typedef struct {
-  float *in;
-  float *window;
-
-  chirp *c;
-  int *ret;
-  int max_y;
-} vearg;
-
-
-pthread_mutex_t ymutex = PTHREAD_MUTEX_INITIALIZER;
-int next_y=0;
-
-#define _GNU_SOURCE
-#include <fenv.h>
-
-void *w_e_column(void *in){
-  float rec[BLOCKSIZE];
-  vearg *arg = (vearg *)in;
-  int y;
-  chirp save;
-
-  while(1){
-    pthread_mutex_lock(&ymutex);
-    y=next_y;
-    if(y>=arg->max_y){
-      pthread_mutex_unlock(&ymutex);
-      return NULL;
-    }
-    next_y++;
-    pthread_mutex_unlock(&ymutex);
-
-    int except=feenableexcept(FE_ALL_EXCEPT);
-    fedisableexcept(FE_INEXACT);
-    fedisableexcept(FE_UNDERFLOW);
-    save=arg->c[y];
-
-    /* iterate only until convergence */
-    arg->ret[y]=
-      estimate_chirps(arg->in,rec,arg->window,BLOCKSIZE,
-                      arg->c+y,1,iterations,fit_tolerance,
-                      fit_linear,
-                      fit_W,
-                      fit_dA,
-                      fit_dW,
-                      fit_ddA,
-                      fit_symm_norm,
-                      fit_gauss_seidel,
-                      fit_bound_zero);
-
-    /* continue iterating to get error numbers for a fixed large
-       number of iterations.  The linear estimator must be restarted
-       from the beginning, else 'continuing' causes it to recenter the
-       basis-- which renders it nonlinear 
-    if(fit_linear) arg->c[y]=save;
-    int ret=estimate_chirps(arg->in,rec,arg->window,BLOCKSIZE,
-                    arg->c+y,1,fit_linear?iterations:arg->ret[y],-1.,
-                    fit_linear,
-                    fit_W,
-                    fit_dA,
-                    fit_dW,
-                    fit_ddA,
-                    fit_symm_norm,
-                    fit_gauss_seidel,
-                    fit_bound_zero);*/
-    arg->ret[y] = iterations-arg->ret[y];
-    feclearexcept(FE_ALL_EXCEPT);
-    feenableexcept(except);
-
+void destroy_page(cairo_t *c){
+  if(c){
+    cairo_surface_t *cs=cairo_get_target(c);
+    cairo_destroy(c);
+    cairo_surface_destroy(cs);
   }
-  return NULL;
 }
-
-void w_e(){
-  float window[BLOCKSIZE];
-  float in[BLOCKSIZE];
-  int i,x,y;
-
-  cairo_surface_t *converge=NULL;
-  cairo_surface_t *Aerror=NULL;
-  cairo_surface_t *Perror=NULL;
-  cairo_surface_t *Werror=NULL;
-  cairo_surface_t *dAerror=NULL;
-  cairo_surface_t *dWerror=NULL;
-  cairo_surface_t *ddAerror=NULL;
-
-  cairo_t *cC=NULL;
-  cairo_t *cA=NULL;
-  cairo_t *cP=NULL;
-  cairo_t *cW=NULL;
-  cairo_t *cdA=NULL;
-  cairo_t *cdW=NULL;
-  cairo_t *cddA=NULL;
-
-  pthread_t threads[cores];
-  vearg arg[cores];
-
-  char *filebase="W-vs-Westimate";
-  char *yaxis_label = "initial distance from W (cycles/block)";
-  char *xaxis_label = "W (cycles/block)";
-  char *title2=NULL;
-  char *title3=NULL;
-
-  pic_w = x_n;
-  pic_h = y_n*2+1;
-  if(iterations<0)iterations=100;
-
-  /* determine ~ padding needed */
-
-  setup_graph("Convergence",title2,title3,xaxis_label,yaxis_label,
-              "Iterations:",DT_iterations);
-  setup_graph("A (Amplitude) Error",title2,title3,xaxis_label,yaxis_label,
-              "Percentage Error:",DT_percent);
-  setup_graph("P (Phase) Error",title2,title3,xaxis_label,yaxis_label,
-              "Error (rads/block)",DT_abserror);
-  if(fit_W)
-    setup_graph("W (Frequency) Error",title2,title3,xaxis_label,yaxis_label,
-                "Error (cycles/block)",DT_abserror);
-  if(fit_dA)
-    setup_graph("dA (Amplitude Modulation) Error",title2,title3,xaxis_label,yaxis_label,
-                "Percentage Error:",DT_percent);
-  if(fit_dW)
-    setup_graph("dW (Chirp Rate) Error",title2,title3,xaxis_label,yaxis_label,
-                "Error (cycles/block)",DT_abserror);
-  if(fit_ddA)
-    setup_graph("ddA (Amplitude Modulation Squared) Error",title2,title3,xaxis_label,yaxis_label,
-                "Percentage Error:",DT_percent);
-
-  /* Make cairo drawables */
-  converge=cairo_image_surface_create(CAIRO_FORMAT_ARGB32,pic_w, pic_h);
-  if(!converge || cairo_surface_status(converge)!=CAIRO_STATUS_SUCCESS){
-    fprintf(stderr,"Could not set up Cairo surface.\n\n");
-    exit(1);
-  }
-  cC=cairo_create(converge);
-  draw_page(cC,"Convergence",title2,title3,xaxis_label,yaxis_label,
-            "Iterations:",DT_iterations);
-
-  Aerror=cairo_image_surface_create(CAIRO_FORMAT_ARGB32,pic_w, pic_h);
-  if(!Aerror || cairo_surface_status(Aerror)!=CAIRO_STATUS_SUCCESS){
-    fprintf(stderr,"Could not set up Cairo surface.\n\n");
-    exit(1);
-  }
-  cA=cairo_create(Aerror);
-  draw_page(cA,"A (Amplitude) Error",title2,title3,xaxis_label,yaxis_label,
-            "Percentage Error:",DT_percent);
-
-  Perror=cairo_image_surface_create(CAIRO_FORMAT_ARGB32,pic_w, pic_h);
-  if(!Perror || cairo_surface_status(Perror)!=CAIRO_STATUS_SUCCESS){
-    fprintf(stderr,"Could not set up Cairo surface.\n\n");
-    exit(1);
-  }
-  cP=cairo_create(Perror);
-  draw_page(cP,"P (Phase) Error",title2,title3,xaxis_label,yaxis_label,
-            "Error (rads/block)",DT_abserror);
-
-  if(fit_W){
-    Werror=cairo_image_surface_create(CAIRO_FORMAT_ARGB32,pic_w, pic_h);
-    if(!Werror || cairo_surface_status(Werror)!=CAIRO_STATUS_SUCCESS){
-      fprintf(stderr,"Could not set up Cairo surface.\n\n");
-      exit(1);
-    }
-    cW=cairo_create(Werror);
-    draw_page(cW,"W (Frequency) Error",title2,title3,xaxis_label,yaxis_label,
-              "Error (cycles/block)",DT_abserror);
-  }
-
-  if(fit_dA){
-    dAerror=cairo_image_surface_create(CAIRO_FORMAT_ARGB32,pic_w, pic_h);
-    if(!dAerror || cairo_surface_status(dAerror)!=CAIRO_STATUS_SUCCESS){
-      fprintf(stderr,"Could not set up Cairo surface.\n\n");
-      exit(1);
-    }
-    cdA=cairo_create(dAerror);
-    draw_page(cdA,"dA (Amplitude Modulation) Error",
-              title2,title3,xaxis_label,yaxis_label,
-              "Percentage Error:",DT_percent);
-  }
-
-  if(fit_dW){
-    dWerror=cairo_image_surface_create(CAIRO_FORMAT_ARGB32,pic_w, pic_h);
-    if(!dWerror || cairo_surface_status(dWerror)!=CAIRO_STATUS_SUCCESS){
-      fprintf(stderr,"Could not set up Cairo surface.\n\n");
-      exit(1);
-    }
-    cdW=cairo_create(dWerror);
-    draw_page(cdW,"dW (Chirp Rate) Error",title2,title3,xaxis_label,yaxis_label,
-              "Error (cycles/block)",DT_abserror);
-  }
-
-  if(fit_ddA){
-    ddAerror=cairo_image_surface_create(CAIRO_FORMAT_ARGB32,pic_w, pic_h);
-    if(!ddAerror || cairo_surface_status(ddAerror)!=CAIRO_STATUS_SUCCESS){
-      fprintf(stderr,"Could not set up Cairo surface.\n\n");
-      exit(1);
-    }
-    cddA=cairo_create(ddAerror);
-    draw_page(cddA,"ddA (Amplitude Modulation Squared) Error",
-              title2,title3,xaxis_label,yaxis_label,
-              "Percentage Error:",DT_percent);
-  }
-
-  for(i=0;i<BLOCKSIZE;i++)
-    window[i]=1.;
-  if(fit_hanning)
-    hanning(window,BLOCKSIZE);
-
-  /* graph computation */
-  for(x=0;x<x_n;x++){
-    float w = (float)x/xstepd;
-    chirp ch[y_n*2+1];
-    int ret[y_n*2+1];
-
-    float rand_P_ch = 2*M_PI*(float)drand48()-M_PI;
-
-    for(i=0;i<BLOCKSIZE;i++){
-      float jj = i-BLOCKSIZE*.5+.5;
-      in[i]=initial_chirp_A * cos(w*jj*2.*M_PI/BLOCKSIZE+initial_chirp_P);
-    }
-
-    fprintf(stderr,"\rW estimate distance vs. W graphs: %d%%...",100*x/x_n);
-
-    for(y=y_n;y>=-y_n;y--){
-      float rand_A_est = (float)drand48();
-      float rand_P_est = 2*M_PI*(float)drand48()-M_PI;
-      int yi=y_n-y;
-      float we=(float)y/ystepd+w;
-      ch[yi]=(chirp){0,
-                     (we)*2.*M_PI/BLOCKSIZE,
-                     0,
-                     initial_est_dA,
-                     initial_est_dW,
-                     initial_est_ddA,
-                     yi};
-    }
-    next_y=0;
-    for(y=0;y<cores;y++){
-      arg[y].in = in;
-      arg[y].window=window;
-      arg[y].c=ch;
-      arg[y].ret=ret;
-      arg[y].max_y=y_n*2+1;
-      pthread_create(threads+y,NULL,w_e_column,arg+y);
-    }
-    for(y=0;y<cores;y++)
-      pthread_join(threads[y],NULL);
-
-    for(y=-y_n;y<=y_n;y++){
-      int yi=y+y_n;
-      float a = (x%xstepg==0 || y%ystepg==0 ?
-                 (x%xstepd==0 || y%ystepd==0 ? .3 : .8) : 1.);
-
-      set_iter_color(cC,ret[yi],a);
-      cairo_rectangle(cC,x+leftpad,yi+toppad,1,1);
-      cairo_fill(cC);
-
-      set_error_color(cA,fabs(initial_chirp_A-ch[yi].A),a);
-      cairo_rectangle(cA,x+leftpad,yi+toppad,1,1);
-      cairo_fill(cA);
-
-      set_error_color(cP,fabs(initial_chirp_P-ch[yi].P),a);
-      cairo_rectangle(cP,x+leftpad,yi+toppad,1,1);
-      cairo_fill(cP);
-
-      if(cW){
-        set_error_color(cW,fabs(w-(ch[yi].W/2./M_PI*BLOCKSIZE)),a);
-        cairo_rectangle(cW,x+leftpad,yi+toppad,1,1);
-        cairo_fill(cW);
-      }
-
-      if(cdA){
-        set_error_color(cdA,ch[yi].dA*BLOCKSIZE,a);
-        cairo_rectangle(cdA,x+leftpad,yi+toppad,1,1);
-        cairo_fill(cdA);
-      }
-
-      if(cdW){
-        set_error_color(cdW,ch[yi].dW/M_PI*BLOCKSIZE*BLOCKSIZE,a);
-        cairo_rectangle(cdW,x+leftpad,yi+toppad,1,1);
-        cairo_fill(cdW);
-      }
-
-      if(cddA){
-        set_error_color(cddA,ch[yi].ddA/BLOCKSIZE*BLOCKSIZE,a);
-        cairo_rectangle(cddA,x+leftpad,yi+toppad,1,1);
-        cairo_fill(cddA);
-      }
-    }
-
-    if((x&15)==0){
-      to_png(converge,filebase,"converge");
-      to_png(Aerror,filebase,"Aerror");
-      to_png(Perror,filebase,"Perror");
-      to_png(Werror,filebase,"Werror");
-      to_png(dAerror,filebase,"dAerror");
-      to_png(dWerror,filebase,"dWerror");
-      to_png(ddAerror,filebase,"ddAerror");
-    }
-
-  }
-
-  to_png(converge,filebase,"converge");
-  to_png(Aerror,filebase,"Aerror");
-  to_png(Perror,filebase,"Perror");
-  to_png(Werror,filebase,"Werror");
-  to_png(dAerror,filebase,"dAerror");
-  to_png(dWerror,filebase,"dWerror");
-  to_png(ddAerror,filebase,"ddAerror");
-
-}
-
-int main(){
-  w_e();
-  return 0;
-}
-

Added: trunk/ghost/monty/chirp/chirpgraph.h
===================================================================
--- trunk/ghost/monty/chirp/chirpgraph.h	                        (rev 0)
+++ trunk/ghost/monty/chirp/chirpgraph.h	2011-04-18 05:39:44 UTC (rev 17921)
@@ -0,0 +1,46 @@
+/********************************************************************
+ *                                                                  *
+ * THIS FILE IS PART OF THE OggGhost SOFTWARE CODEC SOURCE CODE.    *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
+ *                                                                  *
+ * THE OggGhost SOURCE CODE IS (C) COPYRIGHT 2007-2011              *
+ * by the Xiph.Org Foundation http://www.xiph.org/                  *
+ *                                                                  *
+ ********************************************************************
+
+ function: research-grade chirp extraction code
+ last mod: $Id$
+
+ ********************************************************************/
+
+#include <cairo/cairo.h>
+#define DT_iterations 0
+#define DT_abserror 1
+#define DT_percent 2
+extern int toppad;
+extern int leftpad;
+
+extern void set_error_color(cairo_t *c, float err,float a);
+extern void set_iter_color(cairo_t *cC, int ret, float a);
+extern void to_png(cairo_t *c,char *base, char *name);
+extern void destroy_page(cairo_t *c);
+extern cairo_t *draw_page(char *title,
+                          char *subtitle1,
+                          char *subtitle2,
+                          char *subtitle3,
+                          char *xaxis_label,
+                          char *yaxis_label,
+                          char *legend_label,
+                          int datatype);
+extern void setup_graphs(int start_x_step,
+                         int end_x_step, /* inclusive; not one past */
+                         int x_major_d,
+
+                         int start_y_step,
+                         int end_y_step, /* inclusive; not one past */
+                         int y_major_d,
+
+                         int subtitles,
+                         float fs);

Added: trunk/ghost/monty/chirp/chirptest.c
===================================================================
--- trunk/ghost/monty/chirp/chirptest.c	                        (rev 0)
+++ trunk/ghost/monty/chirp/chirptest.c	2011-04-18 05:39:44 UTC (rev 17921)
@@ -0,0 +1,821 @@
+/********************************************************************
+ *                                                                  *
+ * THIS FILE IS PART OF THE OggGhost SOFTWARE CODEC SOURCE CODE.    *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
+ *                                                                  *
+ * THE OggGhost SOURCE CODE IS (C) COPYRIGHT 2007-2011              *
+ * by the Xiph.Org Foundation http://www.xiph.org/                  *
+ *                                                                  *
+ ********************************************************************
+
+ function: graphing code for chirp tests
+ last mod: $Id$
+
+ ********************************************************************/
+
+#define _GNU_SOURCE
+#include <math.h>
+#include "chirp.h"
+#include "chirpgraph.h"
+#include "scales.h"
+#include "window.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <cairo/cairo.h>
+#include <pthread.h>
+
+float circular_distance(float A,float B){
+  float ret = A-B;
+  while(ret<-M_PI)ret+=2*M_PI;
+  while(ret>=M_PI)ret-=2*M_PI;
+  return ret;
+}
+
+/* ranges are inclusive */
+void set_chirp(chirp *c,int rand_p, int i, int n,
+               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;
+  }
+}
+
+/*********************** Plot w estimate error against w **********************/
+
+typedef struct {
+  float *in;
+  float *window;
+  int blocksize;
+  int max_iterations;
+  float fit_tolerance;
+
+  int fit_gauss_seidel;
+  int fit_W;
+  int fit_dA;
+  int fit_dW;
+  int fit_ddA;
+  int fit_nonlinear; /* 0==linear, 1==W-recentered, 2==W,dW-recentered */
+  float fit_W_alpha;
+  float fit_dW_alpha;
+  int fit_symm_norm;
+  int fit_bound_zero;
+
+  chirp *chirp;
+  chirp *estimate;
+  float *rms_error;
+  int *iterations;
+
+} colarg;
+
+
+pthread_mutex_t ymutex = PTHREAD_MUTEX_INITIALIZER;
+int next_y=0;
+int max_y=0;
+
+#define _GNU_SOURCE
+#include <fenv.h>
+
+void *compute_column(void *in){
+  colarg *arg = (colarg *)in;
+  int blocksize=arg->blocksize;
+  float rec[blocksize];
+  int y,i,ret;
+  int except;
+  int localinit = !arg->in;
+  float *chirp = localinit ? malloc(sizeof(*chirp)*blocksize) : arg->in;
+
+  while(1){
+    float s=0.;
+
+    pthread_mutex_lock(&ymutex);
+    y=next_y;
+    if(y>=max_y){
+      pthread_mutex_unlock(&ymutex);
+      return NULL;
+    }
+    next_y++;
+    pthread_mutex_unlock(&ymutex);
+
+    /* if the input is uninitialized, it's because we're sweeping or
+       randomizing chirp components across the column; generate the
+       input here in the thread */
+    if(localinit){
+      for(i=0;i<blocksize;i++){
+        float jj = i - blocksize/2 + .5;
+        float A = arg->chirp[y].A + (arg->chirp[y].dA + arg->chirp[y].ddA*jj)*jj;
+        float P = arg->chirp[y].P + (arg->chirp[y].W  + arg->chirp[y].dW *jj)*jj;
+        chirp[i] = A*cosf(P);
+      }
+    }
+
+    except=feenableexcept(FE_ALL_EXCEPT);
+    fedisableexcept(FE_INEXACT);
+    fedisableexcept(FE_UNDERFLOW);
+
+    ret=estimate_chirps(chirp,rec,arg->window,blocksize,
+                        arg->estimate+y,1,arg->fit_tolerance,arg->max_iterations,
+                        arg->fit_gauss_seidel,
+                        arg->fit_W,
+                        arg->fit_dA,
+                        arg->fit_dW,
+                        arg->fit_ddA,
+                        arg->fit_nonlinear,
+                        arg->fit_W_alpha,
+                        arg->fit_dW_alpha,
+                        arg->fit_symm_norm,
+                        arg->fit_bound_zero);
+
+    for(i=0;i<blocksize;i++){
+      float v = (chirp[i]-rec[i])*arg->window[i];
+      s += v*v;
+    }
+    arg->rms_error[y] = sqrt(s/blocksize);
+    arg->iterations[y] = arg->max_iterations-ret;
+    feclearexcept(FE_ALL_EXCEPT);
+    feenableexcept(except);
+
+  }
+
+  if(localinit)free(chirp);
+  return NULL;
+}
+
+typedef struct {
+  float fontsize;
+  char *subtitle1;
+  char *subtitle2;
+  char *subtitle3;
+
+  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;
+
+  int fit_gauss_seidel;
+  int fit_W;
+  int fit_dA;
+  int fit_dW;
+  int fit_ddA;
+  int fit_nonlinear;
+  float fit_W_alpha;
+  float fit_dW_alpha;
+  int fit_symm_norm;
+  int fit_bound_zero;
+
+  /* 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;
+
+  float min_est_P;
+  float max_est_P;
+
+  float min_est_W;
+  float max_est_W;
+
+  float min_est_dA;
+  float max_est_dA;
+
+  float min_est_dW;
+  float max_est_dW;
+
+  float min_est_ddA;
+  float max_est_ddA;
+
+  float min_chirp_A;
+  float max_chirp_A;
+
+  float min_chirp_P;
+  float max_chirp_P;
+
+  float min_chirp_W;
+  float max_chirp_W;
+
+  float min_chirp_dA;
+  float max_chirp_dA;
+
+  float min_chirp_dW;
+  float max_chirp_dW;
+
+  float min_chirp_ddA;
+  float max_chirp_ddA;
+
+}  graph_run;
+
+/* performs a W initial estimate error vs chirp W plot.  Ignores the
+   est and chirp arguments for W; these are pulled from the x and y setup */
+
+void w_e(char *filebase,graph_run *arg){
+  int threads=arg->threads;
+  int blocksize = arg->blocksize;
+  float window[blocksize];
+  float in[blocksize];
+  int i,xi,yi;
+
+  int x_n = arg->x1-arg->x0+1;
+  int y_n = arg->y1-arg->y0+1;
+
+  /* graphs:
+
+     convergence
+     Aerror
+     Perror
+     Werror
+     dAerror
+     dWerror
+     ddAerror
+     rms fit error
+
+     generate for:
+     worst case across sweep (0-7)
+     delta across sweep (8-15)
+  */
+
+  cairo_t *cC=NULL;
+  cairo_t *cA=NULL;
+  cairo_t *cP=NULL;
+  cairo_t *cW=NULL;
+  cairo_t *cdA=NULL;
+  cairo_t *cdW=NULL;
+  cairo_t *cddA=NULL;
+  cairo_t *cRMS=NULL;
+
+  cairo_t *cC_d=NULL;
+  cairo_t *cA_d=NULL;
+  cairo_t *cP_d=NULL;
+  cairo_t *cW_d=NULL;
+  cairo_t *cdA_d=NULL;
+  cairo_t *cdW_d=NULL;
+  cairo_t *cddA_d=NULL;
+  cairo_t *cRMS_d=NULL;
+
+  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;
+
+  float ret_minA[y_n];
+  float ret_minP[y_n];
+  float ret_minW[y_n];
+  float ret_mindA[y_n];
+  float ret_mindW[y_n];
+  float ret_minddA[y_n];
+  float ret_minRMS[y_n];
+  float ret_miniter[y_n];
+
+  float ret_maxA[y_n];
+  float ret_maxP[y_n];
+  float ret_maxW[y_n];
+  float ret_maxdA[y_n];
+  float ret_maxdW[y_n];
+  float ret_maxddA[y_n];
+  float ret_maxRMS[y_n];
+  float ret_maxiter[y_n];
+
+  /* determine ~ padding needed */
+  setup_graphs(arg->x0,arg->x1,arg->xmajor,
+               arg->y0,arg->y1,arg->ymajor,
+               (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 ||
+     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 ||
+     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)
+    cC_d = draw_page("Convergence Delta",
+                     arg->subtitle1,
+                     arg->subtitle2,
+                     arg->subtitle3,
+                     xaxis_label,
+                     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)
+    cA_d = draw_page("A (Amplitude) Delta",
+                     arg->subtitle1,
+                     arg->subtitle2,
+                     arg->subtitle3,
+                     xaxis_label,
+                     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)
+    cP_d = draw_page("Phase Delta",
+                     arg->subtitle1,
+                     arg->subtitle2,
+                     arg->subtitle3,
+                     xaxis_label,
+                     yaxis_label,
+                     "Delta (rads/block):",
+                     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)
+      cW_d = draw_page("Frequency Delta",
+                       arg->subtitle1,
+                       arg->subtitle2,
+                       arg->subtitle3,
+                       xaxis_label,
+                       yaxis_label,
+                       "Delta (rads/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)
+      cdA_d = draw_page("Amplitude Modulation Delta",
+                        arg->subtitle1,
+                        arg->subtitle2,
+                        arg->subtitle3,
+                        xaxis_label,
+                        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)
+      cdW_d = draw_page("Chirp Rate Delta",
+                        arg->subtitle1,
+                        arg->subtitle2,
+                        arg->subtitle3,
+                        xaxis_label,
+                        yaxis_label,
+                        "Delta (rads/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)
+      cddA_d = draw_page("Amplitude Modulation Squared Delta",
+                         arg->subtitle1,
+                         arg->subtitle2,
+                         arg->subtitle3,
+                         xaxis_label,
+                         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)
+    cRMS_d = draw_page("RMS Error Delta",
+                       arg->subtitle1,
+                       arg->subtitle2,
+                       arg->subtitle3,
+                       xaxis_label,
+                       yaxis_label,
+                       "Percentage Delta:",
+                       DT_percent);
+
+  if(arg->window)
+    arg->window(window,blocksize);
+  else
+    for(i=0;i<blocksize;i++)
+      window[i]=1.;
+
+  /* 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,x_n-1);
+
+    memset(targ,0,sizeof(targ));
+
+    for(i=0;i<threads;i++){
+      if(!chirp_swept)targ[i].in=in;
+      targ[i].window=window;
+      targ[i].blocksize=blocksize;
+      targ[i].max_iterations=100;
+      targ[i].fit_tolerance=arg->fit_tolerance;
+      targ[i].fit_gauss_seidel=arg->fit_gauss_seidel;
+      targ[i].fit_W=arg->fit_W;
+      targ[i].fit_dA=arg->fit_dA;
+      targ[i].fit_dW=arg->fit_dW;
+      targ[i].fit_ddA=arg->fit_ddA;
+      targ[i].fit_nonlinear=arg->fit_nonlinear;
+      targ[i].fit_W_alpha=arg->fit_W_alpha;
+      targ[i].fit_dW_alpha=arg->fit_dW_alpha;
+      targ[i].fit_symm_norm=arg->fit_symm_norm;
+      targ[i].fit_bound_zero=arg->fit_bound_zero;
+      targ[i].chirp=chirps;
+      //targ[i].in=NULL;
+      targ[i].estimate=estimates;
+      targ[i].rms_error=rms;
+      targ[i].iterations=iter;
+    }
+
+    /* if we're sweeping a parameter, we're going to iterate here for a bit. */
+    for(si=0;si<sn;si++){
+      max_y=y_n;
+      next_y=0;
+
+      /* 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;
+
+        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);
+      }
+
+      if(!chirp_swept){
+        for(i=0;i<threads;i++)
+          targ[i].in=in;
+
+        for(i=0;i<blocksize;i++){
+          float jj = i-blocksize/2+.5;
+          float A = chirps[0].A + (chirps[0].dA + chirps[0].ddA*jj)*jj;
+          float P = chirps[0].P + (chirps[0].W  + chirps[0].dW *jj)*jj;
+          in[i] = A*cosf(P);
+        }
+      }
+
+      /* compute column */
+      for(i=0;i<threads;i++)
+        pthread_create(threadlist+i,NULL,compute_column,targ+i);
+      for(i=0;i<threads;i++)
+        pthread_join(threadlist[i],NULL);
+
+      /* accumulate results for this pass */
+      if(si==0){
+        for(i=0;i<y_n;i++){
+          ret_minA[i]=ret_maxA[i]=chirps[i].A - estimates[i].A;
+          ret_minP[i]=ret_maxP[i]=circular_distance(chirps[i].P,estimates[i].P);
+          ret_minW[i]=ret_maxW[i]=chirps[i].W - estimates[i].W;
+          ret_mindA[i]=ret_maxdA[i]=chirps[i].dA - estimates[i].dA;
+          ret_mindW[i]=ret_maxdW[i]=chirps[i].dW - estimates[i].dW;
+          ret_minddA[i]=ret_maxddA[i]=chirps[i].ddA - estimates[i].ddA;
+          ret_minRMS[i]=ret_maxRMS[i]=rms[i];
+          ret_miniter[i]=ret_maxiter[i]=iter[i];
+        }
+      }else{
+        for(i=0;i<y_n;i++){
+          float v = chirps[i].A - estimates[i].A;
+          if(ret_minA[i]>v)ret_minA[i]=v;
+          if(ret_maxA[i]<v)ret_maxA[i]=v;
+
+          v = chirps[i].P - estimates[i].P;
+          if(ret_minP[i]>v)ret_minP[i]=v;
+          if(ret_maxP[i]<v)ret_maxP[i]=v;
+
+          v = chirps[i].W - estimates[i].W;
+          if(ret_minW[i]>v)ret_minW[i]=v;
+          if(ret_maxW[i]<v)ret_maxW[i]=v;
+
+          v = chirps[i].dA - estimates[i].dA;
+          if(ret_mindA[i]>v)ret_mindA[i]=v;
+          if(ret_maxdA[i]<v)ret_maxdA[i]=v;
+
+          v = chirps[i].dW - estimates[i].dW;
+          if(ret_mindW[i]>v)ret_mindW[i]=v;
+          if(ret_maxdW[i]<v)ret_maxdW[i]=v;
+
+          v = chirps[i].ddA - estimates[i].ddA;
+          if(ret_minddA[i]>v)ret_minddA[i]=v;
+          if(ret_maxddA[i]<v)ret_maxddA[i]=v;
+
+          if(ret_minRMS[i]>rms[i])ret_minRMS[i]=rms[i];
+          if(ret_maxRMS[i]<rms[i])ret_maxRMS[i]=rms[i];
+
+          if(ret_miniter[i]>iter[i])ret_miniter[i]=iter[i];
+          if(ret_maxiter[i]<iter[i])ret_maxiter[i]=iter[i];
+        }
+      }
+    }
+
+    /* Column sweep complete; plot */
+    for(yi=0;yi<y_n;yi++){
+      int y=arg->y1-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;
+
+      /* Convergence graph */
+      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){
+        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);
+
+      /* A delta graph */
+      if(swept){
+        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);
+
+      /* P delta graph */
+      if(swept){
+        set_error_color(cP_d,ret_maxP[yi]-ret_minP[yi],a);
+        cairo_rectangle(cP_d,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cP_d);
+      }
+
+      if(cW){
+        /* W error graph */
+        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(cdA){
+        /* dA error graph */
+        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(cdW){
+        /* dW error graph */
+        set_error_color(cdW,MAX(fabs(ret_maxdW[yi]),fabs(ret_mindW[yi]))/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);
+        }
+      }
+
+      if(cddA){
+        /* ddA error graph */
+        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);
+        }
+      }
+
+      /* RMS error graph */
+      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){
+        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);
+      }
+    }
+
+    if((x&15)==0){
+      to_png(cC,filebase,"converge");
+      to_png(cA,filebase,"Aerror");
+      to_png(cP,filebase,"Perror");
+      to_png(cW,filebase,"Werror");
+      to_png(cdA,filebase,"dAerror");
+      to_png(cdW,filebase,"dWerror");
+      to_png(cddA,filebase,"ddAerror");
+      to_png(cRMS,filebase,"RMSerror");
+
+      to_png(cC_d,filebase,"convdelta");
+      to_png(cA_d,filebase,"Adelta");
+      to_png(cP_d,filebase,"Pdelta");
+      to_png(cW_d,filebase,"Wdelta");
+      to_png(cdA_d,filebase,"dAdelta");
+      to_png(cdW_d,filebase,"dWdelta");
+      to_png(cddA_d,filebase,"ddAdelta");
+      to_png(cRMS_d,filebase,"RMSdelta");
+    }
+
+  }
+
+  fprintf(stderr," done\n");
+}
+
+int main(){
+  graph_run arg={
+    /* fontsize */      12,
+    /* subtitle1 */     "graphtest1",
+    /* subtitle2 */     "graphtest2",
+    /* subtitle3 */     "graphtest3",
+    /* blocksize */     2048,
+    /* threads */       8,
+
+    /* x0 */            0,
+    /* x1 */            800,
+    /* xmajor */        100,
+    /* xminor */        25,
+
+    /* y0 */            -300,
+    /* y1 */            300,
+    /* ymajor */        100,
+    /* yminor */        25,
+
+    /* window */        window_functions.hanning,
+    /* fit_tol */       .000001,
+    /* gauss_seidel */  1,
+    /* fit_W */         1,
+    /* fit_dA */        1,
+    /* fit_dW */        1,
+    /* fit_ddA */       0,
+    /* nonlinear */     2,
+    /* W_alpha */       1.,
+    /* dW_alpha */      1.75,
+    /* symm_norm */     1,
+    /* bound_zero */    0,
+
+    /* sweep_steps */   0,
+    /* 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.,
+
+    /* ch A range */    1.,1.,
+    /* ch P range */    M_PI/2.,M_PI/2.,
+    /* ch W range */    0.,0.,
+    /* ch dA range */   0.,0.,
+    /* ch dW range */   0.,0.,
+    /* ch ddA range */  0.,0.,
+  };
+
+  w_e("graphs-A",&arg);
+  return 0;
+}
+

Modified: trunk/ghost/monty/chirp/window.c
===================================================================
--- trunk/ghost/monty/chirp/window.c	2011-04-12 06:45:19 UTC (rev 17920)
+++ trunk/ghost/monty/chirp/window.c	2011-04-18 05:39:44 UTC (rev 17921)
@@ -1,38 +1,113 @@
+/********************************************************************
+ *                                                                  *
+ * THIS FILE IS PART OF THE OggGhost SOFTWARE CODEC SOURCE CODE.    *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
+ *                                                                  *
+ * THE OggGhost SOURCE CODE IS (C) COPYRIGHT 2007-2011              *
+ * by the Xiph.Org Foundation http://www.xiph.org/                  *
+ *                                                                  *
+ ********************************************************************
+
+ function: window functions for research code
+ last mod: $Id$
+
+ ********************************************************************/
+
 #define _GNU_SOURCE
 #include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <cairo/cairo.h>
-#include <squishyio/squishyio.h>
-#include "sinusoids.h"
-#include "smallft.h"
-#include "scales.h"
+#include "window.h"
 
 /* A few windows */
 
+/* Not-a-window */
+static void rectangle(float *x, int n){
+  int i;
+  for(i=0;i<n;i++)
+    x[i]=1.;
+}
+
+/* Minimum 4-term Blackman Harris; highest sidelobe = -92dB */
 #define A0 .35875f
 #define A1 .48829f
 #define A2 .14128f
 #define A3 .01168f
 
-void blackmann_harris(float *out, float *in, int n){
+static void blackmann_harris(float *x, int n){
   int i;
   float scale = 2*M_PI/n;
 
   for(i=0;i<n;i++){
     float i5 = i+.5;
     float w = A0 - A1*cos(scale*i5) + A2*cos(scale*i5*2) - A3*cos(scale*i5*3);
-    out[i] = in[i]*w;
+    x[i] = w;
   }
 }
 
-static void hanning(float *out, float *in, int n){
+/* Good 'ol Hanning window (narrow mainlobe, fast falloff, high sidelobes) */
+static void hanning(float *x, int n){
   int i;
   float scale = 2*M_PI/n;
 
   for(i=0;i<n;i++){
     float i5 = i+.5;
-    out[i] = in[i]*(.5-.5*cos(scale*i5));
+    x[i] = (.5-.5*cos(scale*i5));
   }
 }
+
+/* Triangular * gaussian (sidelobeless window) */
+#define TGA 1.e-4
+#define TGB 21.6
+static void tgauss_deep(float *x, int n){
+  int i;
+  for(i=0;i<n;i++){
+    float f = (i-n/2.)/(n/2.);
+    x[i] = exp(-TGB*pow(f,2)) * pow(1.-fabs(f),TGA);
+  }
+}
+
+static void vorbis(float *d, int n){
+  int i;
+  for(i=0;i<n;i++)
+    d[i] = sin(0.5 * M_PI * sin((i+.5)/n * M_PI)*sin((i+.5)/n * M_PI));
+}
+
+static float beta(int n, float alpha){
+  return cosh (acosh(pow(10,alpha))/(n-1));
+}
+
+static double T(double n, double x){
+  if(fabs(x)<=1){
+    return cos(n*acos(x));
+  }else{
+    return cosh(n*acosh(x));
+  }
+}
+
+/* Dolph-Chebyshev window (a=6., all sidelobes < -120dB) */
+static void dolphcheb(float *d, int n){
+  int i,k;
+  float a = 6.;
+  int M=n/2;
+  int N=M*2;
+  double b = beta(N,a);
+
+  for(i=0;i<n;i++){
+    double sum=0;
+    for(k=0;k<M;k++)
+      sum += (k&1?-1:1)*T(N,b*cos(M_PI*k/N)) * cos (2*i*k*M_PI/N);
+    sum /= T(N,b);
+    sum-=.5;
+    d[i]=sum;
+  }
+}
+
+window_bundle window_functions = {
+  rectangle,
+  hanning,
+  vorbis,
+  blackmann_harris,
+  tgauss_deep,
+  dolphcheb
+};



More information about the commits mailing list