Commit a99f7874 authored by Monty Montgomery's avatar Monty Montgomery

Oops, forgot to commit a bunch of work from a week or two ago

Add preconditioner to CG declip filter


git-svn-id: https://svn.xiph.org/trunk/postfish@5759 0101bb08-14d6-0310-b084-bc0e0c8e3800
parent 39c11f95
......@@ -41,8 +41,8 @@ GtkWidget **feedback_bars;
GtkWidget *samplereadout;
GtkWidget *msreadout;
GtkWidget *hzreadout;
GtkWidget *creadout;
GtkWidget *ireadout;
GtkWidget *depth_readout;
GtkWidget *limit_readout;
GtkWidget *mainpanel_inbar;
......@@ -85,19 +85,23 @@ static void blocksize_slider_change(GtkWidget *w,gpointer in){
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 @@ void clippanel_create(postfish_mainpanel *mp,
"_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 @@ void clippanel_create(postfish_mainpanel *mp,
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);
......
......@@ -217,14 +217,14 @@ static void declip(double *data,double *lap,double *out,
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 @@ static void declip(double *data,double *lap,double *out,
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 @@ time_linkage *declip_read(time_linkage *in){
}
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);
......
......@@ -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 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];
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);
phi_new=inner_product(r,d,n);
res_new=res_0=inner_product(Atb,Atb,n);
for(i=0;i<max && phi_new>e*e*phi_0;i++){
compute_AtAx(fft,d,w,flag,1,n,q);
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));
}
......@@ -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);
#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] */
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment