[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