[xiph-cvs] cvs commit: postfish clippanel.c declip.c reconstruct.c reconstruct.h version.h
Monty
xiphmont at xiph.org
Wed Jan 14 08:03:35 PST 2004
xiphmont 04/01/14 11:03:34
Modified: . clippanel.c declip.c reconstruct.c reconstruct.h
version.h
Log:
Oops, forgot to commit a bunch of work from a week or two ago
Add preconditioner to CG declip filter
Revision Changes Path
1.6 +70 -28 postfish/clippanel.c
Index: clippanel.c
===================================================================
RCS file: /usr/local/cvsroot/postfish/clippanel.c,v
retrieving revision 1.5
retrieving revision 1.6
diff -u -r1.5 -r1.6
--- clippanel.c 26 Dec 2003 09:55:57 -0000 1.5
+++ clippanel.c 14 Jan 2004 16:03:34 -0000 1.6
@@ -41,8 +41,8 @@
GtkWidget *samplereadout;
GtkWidget *msreadout;
GtkWidget *hzreadout;
-GtkWidget *creadout;
-GtkWidget *ireadout;
+GtkWidget *depth_readout;
+GtkWidget *limit_readout;
GtkWidget *mainpanel_inbar;
@@ -85,19 +85,23 @@
declip_setblock(blocksize);
}
-static void converge_slider_change(GtkWidget *w,gpointer in){
+static void depth_slider_change(GtkWidget *w,gpointer in){
char buffer[80];
- double percent=gtk_range_get_value(GTK_RANGE(w));
- double sigfigs=percent*.05+2.8;
- double epsilon=pow(10.,-sigfigs);
+ double dB=gtk_range_get_value(GTK_RANGE(w));
+
+ sprintf(buffer,"%3ddB",(int)dB);
+ readout_set(READOUT(depth_readout),buffer);
+
+ declip_setconvergence(fromdB(-dB));
+}
- sprintf(buffer,"%3.1f",sigfigs);
- readout_set(READOUT(creadout),buffer);
+static void limit_slider_change(GtkWidget *w,gpointer in){
+ char buffer[80];
+ double percent=gtk_range_get_value(GTK_RANGE(w));
- sprintf(buffer,"%3.0f%%",percent);
- readout_set(READOUT(ireadout),buffer);
+ sprintf(buffer,"%3d%%",(int)percent);
+ readout_set(READOUT(limit_readout),buffer);
- declip_setconvergence(epsilon);
declip_setiterations(percent*.01);
}
@@ -114,16 +118,21 @@
"_Declipping filter setup"," [d] ");
GtkWidget *framebox=gtk_hbox_new(1,0);
+ GtkWidget *framebox_right=gtk_vbox_new(0,0);
GtkWidget *blocksize_box=gtk_vbox_new(0,0);
GtkWidget *blocksize_frame=gtk_frame_new (" filter width ");
GtkWidget *converge_frame=gtk_frame_new (" filter convergence ");
+ GtkWidget *limit_frame=gtk_frame_new (" filter CPU throttle ");
GtkWidget *converge_box=gtk_vbox_new(0,0);
+ GtkWidget *limit_box=gtk_vbox_new(0,0);
GtkWidget *channel_table=gtk_table_new(input_ch,4,0);
gtk_widget_set_name(blocksize_box,"choiceframe");
gtk_widget_set_name(converge_box,"choiceframe");
+ gtk_widget_set_name(limit_box,"choiceframe");
gtk_container_set_border_width(GTK_CONTAINER(blocksize_box),2);
gtk_container_set_border_width(GTK_CONTAINER(converge_box),2);
+ gtk_container_set_border_width(GTK_CONTAINER(limit_box),2);
feedback_bars=calloc(input_ch,sizeof(*feedback_bars));
@@ -165,48 +174,81 @@
G_CALLBACK (slider_keymodify), NULL);
g_signal_connect_after (G_OBJECT(slider), "value-changed",
G_CALLBACK(blocksize_slider_change), 0);
- gtk_range_set_value(GTK_RANGE(slider),3.);
+ gtk_range_set_value(GTK_RANGE(slider),2.);
}
gtk_container_add(GTK_CONTAINER(blocksize_frame),blocksize_box);
/* set up convergence config */
{
- GtkWidget *table=gtk_table_new(3,2,0);
+ GtkWidget *table=gtk_table_new(2,2,0);
GtkWidget *sliderbox=gtk_hbox_new(0,0);
GtkWidget *fastlabel=gtk_label_new("fastest");
GtkWidget *qualitylabel=gtk_label_new("best");
- GtkWidget *slider=gtk_hscale_new_with_range(10,200,1);
- GtkWidget *clabel=gtk_label_new("significant figures target");
- GtkWidget *ilabel=gtk_label_new("limit predicted iterations");
- creadout=readout_new("000 ");
- ireadout=readout_new("000%");
+ GtkWidget *slider=gtk_hscale_new_with_range(1,140,1);
+ GtkWidget *label=gtk_label_new("solution depth");
+ depth_readout=readout_new("000dB");
gtk_scale_set_draw_value(GTK_SCALE(slider),FALSE);
- gtk_misc_set_alignment(GTK_MISC(clabel),1,.5);
- gtk_misc_set_alignment(GTK_MISC(ilabel),1,.5);
+ gtk_misc_set_alignment(GTK_MISC(label),1,.5);
gtk_box_pack_start(GTK_BOX(sliderbox),fastlabel,0,0,4);
gtk_box_pack_start(GTK_BOX(sliderbox),slider,1,1,0);
gtk_box_pack_start(GTK_BOX(sliderbox),qualitylabel,0,0,4);
gtk_table_attach(GTK_TABLE(table),sliderbox,0,2,0,1,GTK_FILL|GTK_EXPAND,0,0,8);
- gtk_table_attach(GTK_TABLE(table),clabel,0,1,1,2,GTK_FILL|GTK_EXPAND,GTK_FILL,0,0);
- gtk_table_attach(GTK_TABLE(table),ilabel,0,1,2,3,GTK_FILL|GTK_EXPAND,GTK_FILL,0,0);
+ gtk_table_attach(GTK_TABLE(table),label,0,1,1,2,GTK_FILL|GTK_EXPAND,GTK_FILL,0,0);
+
+ gtk_table_attach(GTK_TABLE(table),depth_readout,1,2,1,2,GTK_FILL,0,5,0);
- gtk_table_attach(GTK_TABLE(table),creadout,1,2,1,2,GTK_FILL,0,5,0);
- gtk_table_attach(GTK_TABLE(table),ireadout,1,2,2,3,GTK_FILL,0,5,0);
gtk_container_add(GTK_CONTAINER(converge_box),table);
g_signal_connect (G_OBJECT (slider), "key-press-event",
G_CALLBACK (slider_keymodify), NULL);
g_signal_connect_after (G_OBJECT(slider), "value-changed",
- G_CALLBACK(converge_slider_change), 0);
- gtk_range_set_value(GTK_RANGE(slider),25.);
+ G_CALLBACK(depth_slider_change), 0);
+ gtk_range_set_value(GTK_RANGE(slider),60.);
+ }
+
+
+ /* set up limit config */
+ {
+ GtkWidget *table=gtk_table_new(2,2,0);
+ GtkWidget *sliderbox=gtk_hbox_new(0,0);
+ GtkWidget *fastlabel=gtk_label_new("fastest");
+ GtkWidget *qualitylabel=gtk_label_new("best");
+ GtkWidget *slider=gtk_hscale_new_with_range(1,100,1);
+ GtkWidget *label=gtk_label_new("hard iteration limit");
+ limit_readout=readout_new("000%");
+
+ gtk_scale_set_draw_value(GTK_SCALE(slider),FALSE);
+ gtk_misc_set_alignment(GTK_MISC(label),1,.5);
+
+ gtk_box_pack_start(GTK_BOX(sliderbox),fastlabel,0,0,4);
+ gtk_box_pack_start(GTK_BOX(sliderbox),slider,1,1,0);
+ gtk_box_pack_start(GTK_BOX(sliderbox),qualitylabel,0,0,4);
+ gtk_table_attach(GTK_TABLE(table),sliderbox,0,2,0,1,GTK_FILL|GTK_EXPAND,0,0,8);
+ gtk_table_attach(GTK_TABLE(table),label,0,1,1,2,GTK_FILL|GTK_EXPAND,GTK_FILL,0,0);
+
+ gtk_table_attach(GTK_TABLE(table),limit_readout,1,2,1,2,GTK_FILL,0,5,0);
+
+ gtk_container_add(GTK_CONTAINER(limit_box),table);
+
+ g_signal_connect (G_OBJECT (slider), "key-press-event",
+ G_CALLBACK (slider_keymodify), NULL);
+ g_signal_connect_after (G_OBJECT(slider), "value-changed",
+ G_CALLBACK(limit_slider_change), 0);
+ gtk_range_set_value(GTK_RANGE(slider),100.);
}
+
+
gtk_container_add(GTK_CONTAINER(converge_frame),converge_box);
+ gtk_container_add(GTK_CONTAINER(limit_frame),limit_box);
+
+ gtk_box_pack_start(GTK_BOX(framebox),blocksize_frame,1,1,4);
+ gtk_box_pack_start(GTK_BOX(framebox),framebox_right,1,1,4);
- gtk_box_pack_start(GTK_BOX(framebox),blocksize_frame,0,1,4);
- gtk_box_pack_start(GTK_BOX(framebox),converge_frame,0,1,4);
+ gtk_box_pack_start(GTK_BOX(framebox_right),converge_frame,1,1,0);
+ gtk_box_pack_start(GTK_BOX(framebox_right),limit_frame,1,1,0);
gtk_box_pack_start(GTK_BOX(panel->subpanel_box),framebox,0,1,4);
gtk_box_pack_start(GTK_BOX(panel->subpanel_box),channel_table,0,1,4);
<p><p>1.4 +9 -13 postfish/declip.c
Index: declip.c
===================================================================
RCS file: /usr/local/cvsroot/postfish/declip.c,v
retrieving revision 1.3
retrieving revision 1.4
diff -u -r1.3 -r1.4
--- declip.c 26 Dec 2003 09:55:57 -0000 1.3
+++ declip.c 14 Jan 2004 16:03:34 -0000 1.4
@@ -217,14 +217,14 @@
double epsilon, double iteration,
int *runningtotal, int *runningcount){
double freq[blocksize];
- int flag[blocksize];
+ double flag[blocksize];
int iterbound,i,j,count=0;
- for(i=0;i<blocksize/8;i++)flag[i]=0;
+ for(i=0;i<blocksize/8;i++)flag[i]=0.;
for(;i<blocksize*7/8;i++){
- flag[i]=0;
+ flag[i]=0.;
if(data[i]>=trigger || data[i]<=-trigger){
- flag[i]=1;
+ flag[i]=1.;
count++;
}
}
@@ -234,21 +234,16 @@
if(declip_active){
- for(;i<blocksize;i++)flag[i]=0;
+ for(;i<blocksize;i++)flag[i]=0.;
for(i=0;i<blocksize;i++)data[i]*=window[i];
memcpy(freq,data,sizeof(freq));
drft_forward(&fft,freq);
sliding_bark_average(freq,freq,blocksize,width);
- if(iteration<1.){
- iterbound=count*iteration;
- }else{
- iterbound=count+blocksize*(iteration-1.);
- }
- if(iterbound<20)iterbound=20;
+ iterbound=count*iteration;
+ if(iterbound<10)iterbound=10;
+ if(count)reconstruct(&fft,data,freq,flag,epsilon,iterbound,blocksize);
- if(count)reconstruct(&fft,data,freq,flag,epsilon*count,iterbound,blocksize);
-
if(out)
for(i=0;i<blocksize/2;i++)
out[i]=lap[i]+data[i]*window[i];
@@ -360,6 +355,7 @@
}
memcpy(work,temp+j,sizeof(*work)*blocksize/2);
memcpy(work+blocksize/2,in->data[i],sizeof(*work)*blocksize/2);
+
declip(work,lap[i],out.data[i]+j,blocksize,
local_trigger[i],local_convergence,local_iterations,
&total,count+i);
<p><p>1.4 +65 -47 postfish/reconstruct.c
Index: reconstruct.c
===================================================================
RCS file: /usr/local/cvsroot/postfish/reconstruct.c,v
retrieving revision 1.3
retrieving revision 1.4
diff -u -r1.3 -r1.4
--- reconstruct.c 20 Dec 2003 11:24:17 -0000 1.3
+++ reconstruct.c 14 Jan 2004 16:03:34 -0000 1.4
@@ -28,83 +28,101 @@
Agonizing Pain' for the additional understanding needed to make the
n^3 -> n^2 log n jump possible. Google for it, you'll find it. */
+#include <string.h>
#include "smallft.h"
#include "reconstruct.h"
-static double inner_product(double *a, double *b, int n){
+static void AtWA(drft_lookup *fft, double *x, double *w,int n){
int i;
- double acc=0.;
- for(i=0;i<n;i++)acc+=a[i]*b[i];
- return acc;
+ drft_forward(fft,x);
+ for(i=0;i<n;i++)x[i]*=w[i];
+ drft_backward(fft,x); /* this is almost the same as A'; see the
+ correction factor rolled into w at the
+ beginning of reconstruct() */
}
-static void compute_AtAx(drft_lookup *fft,
- double *x,double *w,int *flag,int mask,int n,
- double *out){
+/* This is not the inverse of XA'WAX; the algebra isn't valid (due to
+ the singularity of the selection matrix X) and as such is useless
+ for direct solution. However, it does _approximate_ the inverse
+ and as such makes an excellent system preconditioner. */
+static void precondition(drft_lookup *fft, double *x, double *w,int n){
int i;
- if(mask){
- for(i=0;i<n;i++)
- if(!flag[i])
- out[i]=0;
- else
- out[i]=x[i];
- }else
- for(i=0;i<n;i++)
- if(flag[i])
- out[i]=0;
- else
- out[i]=x[i];
-
- drft_forward(fft,out);
- for(i=0;i<n;i++)out[i]*=w[i];
- drft_backward(fft,out);
-
- for(i=0;i<n;i++)
- if(!flag[i])out[i]=0;
-
+ /* no need to remove scaling of result; the relative stretching of
+ the solution space is what's important */
+
+ drft_forward(fft,x); /* almost the same as A^-1'; see the correction
+ factor rolled into w at the beginning of
+ reconstruct() */
+ for(i=0;i<n;i++)x[i]/=w[i];
+ drft_backward(fft,x);
}
-static void compute_Atb_minus_AtAx(drft_lookup *fft,
- double *x,double *w,double *Atb,int *flag,
- int n,double *out){
+static double inner_product(double *a, double *b, int n){
int i;
- compute_AtAx(fft,x,w,flag,1,n,out);
- for(i=0;i<n;i++)out[i]=Atb[i]-out[i];
+ double acc=0.;
+ for(i=0;i<n;i++)acc+=a[i]*b[i];
+ return acc;
}
+#include <stdio.h>
void reconstruct(drft_lookup *fft,
- double *x, double *w, int *flag, double e,int max,int n){
+ double *x, double *w,
+ double *flag, double e,int max,int n){
int i,j;
double Atb[n];
double r[n];
double d[n];
double q[n];
- double phi_new,phi_old,phi_0;
+ double s[n];
+ double phi_new,phi_old,res_0,res_new;
double alpha,beta;
+ /* hack; roll a correction factor for A'/A-1 into w */
+ for(j=1;j<n-1;j++)w[j]*=.5;
+
/* compute initial Atb */
- compute_AtAx(fft,x,w,flag,0,n,Atb);
- for(j=0;j<n;j++)Atb[j]= -Atb[j];
+ for(j=0;j<n;j++)Atb[j]=x[j]*(flag[j]-1.);
+ AtWA(fft,Atb,w,n);
- compute_Atb_minus_AtAx(fft,x,w,Atb,flag,n,r);
- memcpy(d,r,sizeof(d));
- phi_0=phi_new=inner_product(r,r,n);
+ /* compute initial residue */
+ for(j=0;j<n;j++)r[j]=x[j]*flag[j];
+ AtWA(fft,r,w,n);
+ for(j=0;j<n;j++)d[j]=r[j]=(Atb[j]-r[j])*flag[j];
+
+ /* initial preconditioning */
+ precondition(fft,d,w,n);
+ for(j=0;j<n;j++)q[j]=d[j]*=flag[j];
- for(i=0;i<max && phi_new>e*e*phi_0;i++){
- compute_AtAx(fft,d,w,flag,1,n,q);
+ phi_new=inner_product(r,d,n);
+ res_new=res_0=inner_product(Atb,Atb,n);
+
+ for(i=0;i<max && sqrt(res_new)/sqrt(res_0)>e;i++){
+ AtWA(fft,q,w,n);
alpha=phi_new/inner_product(d,q,n);
for(j=0;j<n;j++)x[j]+=alpha*d[j];
- if((i & 0x3f)==0x3f)
- compute_Atb_minus_AtAx(fft,x,w,Atb,flag,n,r);
- else
- for(j=0;j<n;j++)r[j]-=alpha*q[j];
+ if((i & 0x3f)==0x3f){
+ for(j=0;j<n;j++)r[j]=x[j]*flag[j];
+ AtWA(fft,r,w,n);
+ for(j=0;j<n;j++)r[j]=(Atb[j]-r[j])*flag[j];
+ }else
+ for(j=0;j<n;j++)r[j]-=alpha*q[j]*flag[j];
+ /* apply preconditioner */
+ for(j=0;j<n;j++)s[j]=r[j]*flag[j];
+ precondition(fft,s,w,n);
+ for(j=0;j<n;j++)s[j]*=flag[j];
+
phi_old=phi_new;
- phi_new=inner_product(r,r,n);
+ phi_new=inner_product(r,s,n);
+ res_new=inner_product(r,r,n);
beta=phi_new/phi_old;
- for(j=0;j<n;j++) d[j]=r[j]+beta*d[j];
+ for(j=0;j<n;j++) q[j]=d[j]=s[j]+beta*d[j];
+
}
+
+ fprintf(stderr,"converged in %d with res=%f\n",i,sqrt(res_new)/sqrt(res_0));
+
}
<p><p>1.3 +1 -1 postfish/reconstruct.h
Index: reconstruct.h
===================================================================
RCS file: /usr/local/cvsroot/postfish/reconstruct.h,v
retrieving revision 1.2
retrieving revision 1.3
diff -u -r1.2 -r1.3
--- reconstruct.h 20 Dec 2003 11:24:17 -0000 1.2
+++ reconstruct.h 14 Jan 2004 16:03:34 -0000 1.3
@@ -21,5 +21,5 @@
*
*/
-extern void reconstruct(drft_lookup *fft,double *x, double *w, int *flag,
+extern void reconstruct(drft_lookup *fft,double *x, double *w, double *flag,
double e,int max,int n);
<p><p>1.24 +2 -2 postfish/version.h
Index: version.h
===================================================================
RCS file: /usr/local/cvsroot/postfish/version.h,v
retrieving revision 1.23
retrieving revision 1.24
diff -u -r1.23 -r1.24
--- version.h 26 Dec 2003 09:55:57 -0000 1.23
+++ version.h 14 Jan 2004 16:03:34 -0000 1.24
@@ -1,2 +1,2 @@
-#define VERSION "$Id: version.h,v 1.23 2003/12/26 09:55:57 xiphmont Exp $ "
-/* DO NOT EDIT: Automated versioning hack [Fri Dec 26 04:54:29 EST 2003] */
+#define VERSION "$Id: version.h,v 1.24 2004/01/14 16:03:34 xiphmont Exp $ "
+/* DO NOT EDIT: Automated versioning hack [Mon Jan 12 14:46:49 EST 2004] */
<p><p>--- >8 ----
List archives: http://www.xiph.org/archives/
Ogg project homepage: http://www.xiph.org/ogg/
To unsubscribe from this list, send a message to 'cvs-request at xiph.org'
containing only the word 'unsubscribe' in the body. No subject is needed.
Unsubscribe messages sent to the list will be ignored/filtered.
More information about the commits
mailing list