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