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

xiphmont at svn.xiph.org xiphmont at svn.xiph.org
Mon Apr 18 06:30:14 PDT 2011


Author: xiphmont
Date: 2011-04-18 06:30:14 -0700 (Mon, 18 Apr 2011)
New Revision: 17925

Modified:
   trunk/ghost/monty/chirp/chirp.c
   trunk/ghost/monty/chirp/chirpgraph.c
   trunk/ghost/monty/chirp/chirptest.c
Log:
Improve chirp gen precision (input and basis); this will need to be cleaned up a bit

Add '0' to the iteration scale

Correct non-0 starting x-axis and add more spacing code 



Modified: trunk/ghost/monty/chirp/chirp.c
===================================================================
--- trunk/ghost/monty/chirp/chirp.c	2011-04-18 06:47:51 UTC (rev 17924)
+++ trunk/ghost/monty/chirp/chirp.c	2011-04-18 13:30:14 UTC (rev 17925)
@@ -136,7 +136,6 @@
   /* outer fit iteration */
   while(flag && iter_limit>0){
     flag=0;
-    iter_limit--;
 
     /* precompute the portion of the projection/fit estimate shared by
        the zero, first and second order fits.  Subtracts the current
@@ -168,21 +167,25 @@
          of the flag is to emulate/study the behavior of the simplified
          algorithm */
       float cdW = fit_recenter_dW ? c->dW : 0;
+      float Wjj = fmodf((.5f-len/2)*c->W,2.f*M_PI);
 
       for(j=0;j<len;j++){
         float jj = j-len*.5+.5;
         float jj2 = jj*jj;
-        float co,si,c2,s2;
+        float co,si;
+        float c2,s2;
         float yy=r[j];
 
-        sincosf((c->W + cdW*jj)*jj,&si,&co);
+        sincosf(Wjj+cdW*jj*jj,&si,&co);
+        Wjj += c->W;
+        if(Wjj>2.*M_PI)Wjj-=2.*M_PI;
         si*=window[j];
         co*=window[j];
         c2 = co*co*jj;
         s2 = si*si*jj;
 
         /* add the current estimate back to the residue vector */
-        r[j] += (aC*co-aS*si) * (c->A + (c->dA + c->ddA*jj)*jj);
+        r[j] += (aC*co-aS*si) * (c->A + c->dA*jj + c->ddA*jj*jj);
 
         /* zero order projection */
         aP += co*yy;
@@ -317,14 +320,15 @@
       {
         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);
+          double jj = j-len*.5+.5;
+          float a = c->A + c->dA*jj + c->ddA*jj*jj;
+          float v = a*cos(cdW*jj*jj + c->P + c->W*jj);
           r[j] -= v*window[j];
           y[j] += v;
         }
       }
     }
+    if(flag) iter_limit--;
   }
   return iter_limit;
 }
@@ -467,7 +471,7 @@
 
   while(flag && iter_limit){
     flag=0;
-    iter_limit--;
+
     for (i=0;i<n;i++){
 
       float tmpa=0, tmpb=0;
@@ -535,6 +539,7 @@
           (tmpe*tmpe + tmpf*tmpf)/(ei[i]*ei[i]+fi[i]*fi[i]+fit_limit*fit_limit) > fit_limit*fit_limit) flag=1;
 
     }
+    if(flag) iter_limit--;
   }
 
   for(i=0;i<n;i++){

Modified: trunk/ghost/monty/chirp/chirpgraph.c
===================================================================
--- trunk/ghost/monty/chirp/chirpgraph.c	2011-04-18 06:47:51 UTC (rev 17924)
+++ trunk/ghost/monty/chirp/chirpgraph.c	2011-04-18 13:30:14 UTC (rev 17925)
@@ -45,8 +45,10 @@
     cairo_set_source_rgba(cC,0,.6,.6,a);
   }else if (ret==2){  /* 2 dark cyan */
     cairo_set_source_rgba(cC,0,.4,.8,a);
-  }else{ /* 1 blue */
-    cairo_set_source_rgba(cC,.1,.1,1,a);
+  }else if (ret==1){ /* 1 blue */
+    cairo_set_source_rgba(cC,.1,.2,1,a);
+  }else{ /* dark blue */
+    cairo_set_source_rgba(cC,.2,0,1,a);
   }
 }
 
@@ -167,7 +169,7 @@
     float w1,w2,w3;
 
     cairo_text_extents(ct, "100", &extents);
-    w1=extents.width*2*11;
+    w1=extents.width*2*12;
     cairo_text_extents(ct, ".000001", &extents);
     w2=extents.width*1.5*7;
     cairo_text_extents(ct, ".0001%", &extents);
@@ -235,18 +237,24 @@
   }
 
   /* X axis labels */
-  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/xmajor);
+  {
+    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;
+        }
       }
-      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);
     }
   }
 
@@ -325,8 +333,9 @@
       int i;
       float w=cw/11,px=leftpad+x_n;
 
-      for(i=1;i<=100;i++){
+      for(i=0;i<=100;i++){
         switch(i){
+        case 0:
         case 1:
         case 2:
         case 3:

Modified: trunk/ghost/monty/chirp/chirptest.c
===================================================================
--- trunk/ghost/monty/chirp/chirptest.c	2011-04-18 06:47:51 UTC (rev 17924)
+++ trunk/ghost/monty/chirp/chirptest.c	2011-04-18 13:30:14 UTC (rev 17925)
@@ -128,10 +128,10 @@
        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);
+        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;
+        chirp[i] = A*cos(P);
       }
     }
 
@@ -513,7 +513,7 @@
     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);
+    fprintf(stderr,"\rW estimate distance vs. W graphs: column %d/%d...",x-arg->x0,x_n-1);
 
     memset(targ,0,sizeof(targ));
 
@@ -572,10 +572,10 @@
           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);
+          double jj = i-blocksize/2+.5;
+          double A = chirps[0].A + chirps[0].dA*jj + chirps[0].ddA*jj*jj;
+          double P = chirps[0].P + chirps[0].W*jj  + chirps[0].dW *jj*jj;
+          in[i] = A*cos(P);
         }
       }
 
@@ -803,16 +803,16 @@
     /* threads */       8,
 
     /* x0 */            0,
-    /* x1 */            800,
-    /* xmajor */        100,
-    /* xminor */        25,
+    /* x1 */            1280,
+    /* xmajor */        10,
+    /* xminor */        5,
 
     /* y0 */            -300,
     /* y1 */            300,
-    /* ymajor */        100,
-    /* yminor */        25,
+    /* ymajor */        30,
+    /* yminor */        15,
 
-    /* window */        window_functions.hanning,
+    /* window */        window_functions.maxwell1,
     /* fit_tol */       .000001,
     /* gauss_seidel */  1,
     /* fit_W */         1,
@@ -825,8 +825,8 @@
     /* symm_norm */     1,
     /* bound_zero */    0,
 
-    /* sweep_steps */   16,
-    /* randomize_p */   1,
+    /* sweep_steps */   8,
+    /* randomize_p */   0,
 
     /* est A range */   0.,0.,
     /* est P range */   0.,0.,



More information about the commits mailing list