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

xiphmont at svn.xiph.org xiphmont at svn.xiph.org
Wed Apr 27 14:04:50 PDT 2011


Author: xiphmont
Date: 2011-04-27 14:04:50 -0700 (Wed, 27 Apr 2011)
New Revision: 17937

Modified:
   trunk/ghost/monty/chirp/chirp.c
   trunk/ghost/monty/chirp/chirpgraph.c
   trunk/ghost/monty/chirp/chirptest.c
   trunk/ghost/monty/chirp/harness.c
   trunk/ghost/monty/chirp/window.c
   trunk/ghost/monty/chirp/window.h
Log:
minor updated work


Modified: trunk/ghost/monty/chirp/chirp.c
===================================================================
--- trunk/ghost/monty/chirp/chirp.c	2011-04-27 12:38:53 UTC (rev 17936)
+++ trunk/ghost/monty/chirp/chirp.c	2011-04-27 21:04:50 UTC (rev 17937)
@@ -347,7 +347,7 @@
         }
       }
     }
-    if(flag) iter_limit--;
+    iter_limit--;
   }
   return iter_limit;
 }
@@ -397,6 +397,8 @@
   float ei[n], fi[n];
   int i,j;
   int flag=1;
+  float lasterr=0;
+  float thiserr=0;
 
   for (i=0;i<n;i++){
     float tmpa=0;
@@ -488,8 +490,10 @@
     }
   }
 
-  while(flag && iter_limit){
+  while((flag || lasterr>thiserr) && iter_limit>0){
     flag=0;
+    lasterr=thiserr;
+    thiserr=0;
 
     for (i=0;i<n;i++){
 
@@ -543,6 +547,21 @@
           tmpf*ttsin_table[i][j];
       }
 
+
+      /* base convergence on basis projection movement this
+         iteration */
+      {
+        float A2 = ai[i]*ai[i]+bi[i]*bi[i];
+        float move = (tmpa*tmpa + tmpb*tmpb)/(A2 + fit_limit*fit_limit) +
+          (tmpc*tmpc + tmpd*tmpd)/(A2 + fit_limit*fit_limit) +
+          (tmpe*tmpe + tmpf*tmpf)/(A2 + fit_limit*fit_limit);
+        thiserr+=move;
+
+        if(fit_limit>0 && move>fit_limit*fit_limit)flag=1;
+        if(fit_limit<0 && move>1e-14)flag=1;
+      }
+
+
       ai[i] += tmpa;
       bi[i] += tmpb;
       ci[i] += tmpc;
@@ -550,15 +569,18 @@
       ei[i] += tmpe;
       fi[i] += tmpf;
 
-      /* base convergence on basis projection movement this
-         iteration */
-      if( fit_limit<0 ||
-          (tmpa*tmpa + tmpb*tmpb)/(ai[i]*ai[i]+bi[i]*bi[i]+fit_limit*fit_limit) +
-          (tmpc*tmpc + tmpd*tmpd)/(ci[i]*ci[i]+di[i]*di[i]+fit_limit*fit_limit) +
-          (tmpe*tmpe + tmpf*tmpf)/(ei[i]*ei[i]+fi[i]*fi[i]+fit_limit*fit_limit) > fit_limit*fit_limit) flag=1;
+      /* guard overflow; if we're this far out, assume we're never
+         coming back. drop out now. */
+      if((ai[i]*ai[i] + bi[i]*bi[i])>1e10 ||
+         (ci[i]*ci[i] + di[i]*di[i])>1e10 ||
+         (ei[i]*ei[i] + fi[i]*fi[i])>1e10){
+        iter_limit=0;
+        i=n;
+      }
 
+
     }
-    if(flag) iter_limit--;
+    iter_limit--;
   }
 
   for(i=0;i<n;i++){
@@ -672,9 +694,9 @@
   }
 
   if(!nonlinear){
-    if(bound_zero) return -1;
-    if(fit_W_alpha!=1.0) return -1;
-    if(fit_dW_alpha!=1.0) return -1;
+    //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_gs,
                                 fitW,fitdA,fitdW,fitddA,

Modified: trunk/ghost/monty/chirp/chirpgraph.c
===================================================================
--- trunk/ghost/monty/chirp/chirpgraph.c	2011-04-27 12:38:53 UTC (rev 17936)
+++ trunk/ghost/monty/chirp/chirpgraph.c	2011-04-27 21:04:50 UTC (rev 17937)
@@ -211,7 +211,6 @@
     exit(1);
   }
   cairo_t *c = cairo_create(cs);
-  cairo_save(c);
 
   /* clear page to white */
   cairo_set_source_rgb(c,1,1,1);
@@ -220,6 +219,7 @@
   cairo_set_font_size(c, fontsize);
 
   /* set graph area to transparent */
+  cairo_save(c);
   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);

Modified: trunk/ghost/monty/chirp/chirptest.c
===================================================================
--- trunk/ghost/monty/chirp/chirptest.c	2011-04-27 12:38:53 UTC (rev 17936)
+++ trunk/ghost/monty/chirp/chirptest.c	2011-04-27 21:04:50 UTC (rev 17937)
@@ -234,7 +234,7 @@
       rms_acc += ee*ee;
     }
     arg->rms_error[y] = sqrt(rms_acc)/sqrt(e_acc);
-    arg->iterations[y] = arg->max_iterations-ret;
+    arg->iterations[y] = arg->max_iterations-ret-1;
     feclearexcept(FE_ALL_EXCEPT);
     feenableexcept(except);
 
@@ -1094,35 +1094,35 @@
 
 int main(){
   graph_run arg={
-    /* fontsize */      12,
-    /* subtitle1 */     "Nonlinear (W, dW recentered), G/S, symmetric norm, no ddA fit",
+    /* fontsize */      18,
+    /* subtitle1 */     "Linear estimation, no ddA fit",
     /* subtitle2 */     "chirp: A=1.0, dA=0., swept phase | estimate A=P=dA=dW=0, estimate W=chirp W",
-    /* subtitle3 */     "hanning window",
+    /* subtitle3 */     "sine window",
     /* xaxis label */   "W (cycles/block)",
     /* yaxis label */   "dW (cycles/block)",
 
     /* blocksize */     256,
     /* threads */       8,
 
-    /* window */        window_functions.hanning,
+    /* window */        window_functions.sine,
     /* fit_tol */       .000001,
-    /* gauss_seidel */  1,
+    /* gauss_seidel */  0,
     /* fit_W */         1,
     /* fit_dA */        1,
     /* fit_dW */        1,
     /* fit_ddA */       0,
-    /* nonlinear */     2,
+    /* nonlinear */     0,
     /* W_alpha */       1.,
     /* dW_alpha */      1.75,
-    /* symm_norm */     1,
+    /* symm_norm */     0,
     /* bound_zero */    0,
 
     /* x dimension */   DIM_CHIRP_W,
-    /* x steps */       801,
+    /* x steps */       1001,
     /* x major */       1.,
     /* x minor */       .25,
     /* y dimension */   DIM_CHIRP_dW,
-    /* y steps */       501,
+    /* y steps */       601,
     /* y major */       1.,
     /* y minor */       .5,
     /* sweep_steps */   32,
@@ -1137,9 +1137,9 @@
 
     /* ch A range */    1.,1.,
     /* ch P range */    0.,1.-1./32.,
-    /* ch W range */    0.,8.,
+    /* ch W range */    0.,10.,
     /* ch dA range */   0.,0.,
-    /* ch dW range */   -5.,5.,
+    /* ch dW range */   -1.2,1.2,
     /* ch ddA range */  0.,0.,
 
     /* converge max */    1,
@@ -1161,17 +1161,28 @@
 
   };
 
-  w_e("dW-vs-W-rand",&arg);
+  w_e("linear-dW-vs-W",&arg);
+  arg.nonlinear=1;
+  w_e("partial-nonlinear-dW-vs-W",&arg);
+  arg.nonlinear=2;
+  w_e("full-nonlinear-dW-vs-W",&arg);
 
+  arg.nonlinear=0;
   arg.yaxis_label="initial distance from W (cycles/block)";
   arg.y_dim = DIM_ESTIMATE_W;
-  arg.min_est_W = -5.;
-  arg.max_est_W =  5.;
+  arg.min_est_W = -1.2;
+  arg.max_est_W =  1.2;
   arg.min_chirp_dW=-1.;
   arg.max_chirp_dW=1.;
 
-  w_e("estW-vs-W-rand",&arg);
+  w_e("linear-estW-vs-W",&arg);
+  arg.nonlinear=1;
+  w_e("partial-nonlinear-estW-vs-W",&arg);
+  arg.nonlinear=2;
+  w_e("full-nonlinear-estW-vs-W",&arg);
 
+  arg.nonlinear=0;
+
   return 0;
 }
 

Modified: trunk/ghost/monty/chirp/harness.c
===================================================================
--- trunk/ghost/monty/chirp/harness.c	2011-04-27 12:38:53 UTC (rev 17936)
+++ trunk/ghost/monty/chirp/harness.c	2011-04-27 21:04:50 UTC (rev 17937)
@@ -10,7 +10,7 @@
 #include "chirp.h"
 #include "lpc.h"
 
-static int BLOCKSIZE=2048;
+static int BLOCKSIZE=1024;
 static int MAX_CPC=200; /* maximum of 200 chirp candidates per channel */
 
 static float MINdB=-100;
@@ -18,18 +18,18 @@
 
 #define AREAS 2
 static float AREA_SIZE[AREAS]={1.,1.};
-static int   AREA_HZOOM[AREAS]={1,1};
-static int   AREA_VZOOM[AREAS]={16,1};
+static int   AREA_HZOOM[AREAS]={1,4};
+static int   AREA_VZOOM[AREAS]={4,4};
 static int   AREA_FRAMEOUT_TRIGGER=1;
 
-static int   WAVEFORM_AREA=-1;
-static int   SPECTROGRAM_AREA=0;
+static int   WAVEFORM_AREA=0;
+static int   SPECTROGRAM_AREA=-1;
 static int   CHIRPOGRAM_AREA=-1;
 
 static int   SPECTRUM_AREA=1;
 static int   SPECTRUM_GRIDFONTSIZE=8;
 
-static int   ENVELOPE_AREA=-11;
+static int   ENVELOPE_AREA=-1;
 
 static int   TRACK_AREA=1;
 
@@ -347,9 +347,12 @@
   grid_lines(c,SPECTRUM_AREA);
   cairo_set_line_width(c,.7);
 
-  for(i=0;i<pcm->ch;i++)
+  for(i=0;i<pcm->ch;i++){
+    int j;
+    for(j=0;j<BLOCKSIZE/2+1;j++)
+      current_dB[i][j] = current_dB[i][j];
     frequency_vector_dB(current_dB[i],i,BLOCKSIZE/2+1,SPECTRUM_AREA,c,0);
-
+  }
   cairo_restore(c);
 }
 
@@ -387,6 +390,9 @@
 void waveform(pcm_t *pcm,float *window,float **current,cairo_t *c,off_t pos){
   int i,j;
   float work[BLOCKSIZE];
+  drft_lookup fft;
+  drft_init(&fft, BLOCKSIZE);
+
   cairo_save(c);
   clip_area(c,WAVEFORM_AREA);
   clear_area(c,WAVEFORM_AREA);
@@ -395,10 +401,12 @@
   for(i=0;i<pcm->ch;i++){
     for(j=0;j<BLOCKSIZE;j++)
       work[j]= -current[i][j]*window[j];
+
     waveform_vector(work,i,BLOCKSIZE,WAVEFORM_AREA,c,0);
-    waveform_vector(window,i,BLOCKSIZE,WAVEFORM_AREA,c,0);
+    //waveform_vector(window,i,BLOCKSIZE,WAVEFORM_AREA,c,0);
   }
 
+  drft_clear(&fft);
   cairo_restore(c);
 }
 
@@ -702,8 +710,8 @@
     for(j=0;j<chirps_used[i];j++)
       fprintf(stderr,"in>>%f:%f:%f::%f:%f ",chirps[i][j].A,chirps[i][j].W*BLOCKSIZE/2./M_PI,chirps[i][j].P,chirps[i][j].dA*BLOCKSIZE,chirps[i][j].dW*BLOCKSIZE*BLOCKSIZE/2/M_PI);
 
-    ret=estimate_chirps(current_block[i],y[i],r[i],window,BLOCKSIZE,
-                    chirps[i],chirps_used[i],50,.01);
+    //ret=estimate_chirps(current_block[i],y[i],r[i],window,BLOCKSIZE,
+    //                  chirps[i],chirps_used[i],50,.01);
     for(j=0;j<chirps_used[i];j++)
       fprintf(stderr,"out(%d)>>%f:%f:%f::%f:%f   ",ret,chirps[i][j].A,chirps[i][j].W*BLOCKSIZE/2./M_PI,chirps[i][j].P,chirps[i][j].dA*BLOCKSIZE,chirps[i][j].dW*BLOCKSIZE*BLOCKSIZE/2./M_PI);
     dump_waveform("x",fitno,current_block[i],BLOCKSIZE);
@@ -730,8 +738,8 @@
     for(j=0;j<chirps_used[i];j++)
       fprintf(stderr,"in>>%f:%f:%f::%f:%f ",chirps[i][j].A,chirps[i][j].W*BLOCKSIZE/2./M_PI,chirps[i][j].P,chirps[i][j].dA*BLOCKSIZE,chirps[i][j].dW*BLOCKSIZE*BLOCKSIZE/2/M_PI);
 
-    ret=estimate_chirps(current_block[i],y[i],r[i],(count>10?rec:window),BLOCKSIZE,
-                        chirps[i],chirps_used[i],50,.01);
+    //ret=estimate_chirps(current_block[i],y[i],r[i],(count>10?rec:window),BLOCKSIZE,
+    //                  chirps[i],chirps_used[i],50,.01);
 
     for(j=0;j<chirps_used[i];j++)
       fprintf(stderr,"out(%d)>>%f:%f:%f::%f:%f   ",ret,chirps[i][j].A,chirps[i][j].W*BLOCKSIZE/2./M_PI,chirps[i][j].P,chirps[i][j].dA*BLOCKSIZE,chirps[i][j].dW*BLOCKSIZE*BLOCKSIZE/2./M_PI);
@@ -855,7 +863,7 @@
   off_t areas_pos[AREAS]={0,0};
 
   pic_w = BLOCKSIZE/2;
-  pic_h = pic_w*9/16;
+  pic_h = pic_w*9/16/2;
 
   drft_init(&fft, BLOCKSIZE);
 
@@ -917,7 +925,7 @@
     chirps[i]=calloc(MAX_CPC,sizeof(**chirps));
 
   for(j=0;j<BLOCKSIZE;j++)window[j]=1.;
-  hanning(window,window,BLOCKSIZE);
+  //hanning(window,window,BLOCKSIZE);
 
   while(1){
 

Modified: trunk/ghost/monty/chirp/window.c
===================================================================
--- trunk/ghost/monty/chirp/window.c	2011-04-27 12:38:53 UTC (rev 17936)
+++ trunk/ghost/monty/chirp/window.c	2011-04-27 21:04:50 UTC (rev 17937)
@@ -28,6 +28,17 @@
     x[i]=1.;
 }
 
+/* sine window */
+static void sine(float *x, int n){
+  int i;
+  float scale = M_PI/n;
+
+  for(i=0;i<n;i++){
+    float i5 = i+.5;
+    x[i] = sin(scale*i5);
+  }
+}
+
 /* Minimum 4-term Blackman Harris; highest sidelobe = -92dB */
 #define A0 .35875f
 #define A1 .48829f
@@ -120,6 +131,7 @@
 
 window_bundle window_functions = {
   rectangle,
+  sine,
   hanning,
   vorbis,
   blackmann_harris,

Modified: trunk/ghost/monty/chirp/window.h
===================================================================
--- trunk/ghost/monty/chirp/window.h	2011-04-27 12:38:53 UTC (rev 17936)
+++ trunk/ghost/monty/chirp/window.h	2011-04-27 21:04:50 UTC (rev 17937)
@@ -17,6 +17,7 @@
 
 typedef struct {
   void (*rectangle)(float *,int);
+  void (*sine)(float *,int);
   void (*hanning)(float *,int);
   void (*vorbis)(float *,int);
   void (*blackman_harris)(float *,int);



More information about the commits mailing list