[xiph-cvs] cvs commit: postfish declip.c declip.h postfish-gtkrc Makefile input.c lpc.c mainpanel.c minimize.c multibar.c postfish.h reconstruct.c reconstruct.h version.h

Monty xiphmont at xiph.org
Sat Dec 20 03:24:18 PST 2003



xiphmont    03/12/20 06:24:18

  Modified:    .        Makefile input.c lpc.c mainpanel.c minimize.c
                        multibar.c postfish.h reconstruct.c reconstruct.h
                        version.h
  Added:       .        declip.c declip.h postfish-gtkrc
  Log:
  Ongoing filter work; add declipper (not inline yet, but avoid losing work)

Revision  Changes    Path
1.9       +4 -2      postfish/Makefile

Index: Makefile
===================================================================
RCS file: /usr/local/cvsroot/postfish/Makefile,v
retrieving revision 1.8
retrieving revision 1.9
diff -u -r1.8 -r1.9
--- Makefile	18 Oct 2003 09:53:22 -0000	1.8
+++ Makefile	20 Dec 2003 11:24:17 -0000	1.9
@@ -5,8 +5,10 @@
 CC=gcc
 LD=gcc
 
-SRC = main.c mainpanel.c multibar.c readout.c input.c output.c clippanel.c
-OBJ = main.o mainpanel.o multibar.o readout.o input.o output.o clippanel.o
+SRC = main.c mainpanel.c multibar.c readout.c input.c output.c clippanel.c declip.c \
+	reconstruct.c smallft.c
+OBJ = main.o mainpanel.o multibar.o readout.o input.o output.o clippanel.o declip.o \
+	reconstruct.o smallft.o
 GCF = `pkg-config --cflags gtk+-2.0` -DG_DISABLE_DEPRECATED -DGDK_DISABLE_DEPRECATED -DGTK_DISABLE_DEPRECATED -DGDK_PIXBUF_DISABLE_DEPRECATED
 
 all:	

<p><p>1.11      +3 -2      postfish/input.c

Index: input.c
===================================================================
RCS file: /usr/local/cvsroot/postfish/input.c,v
retrieving revision 1.10
retrieving revision 1.11
diff -u -r1.10 -r1.11
--- input.c	17 Oct 2003 08:40:18 -0000	1.10
+++ input.c	20 Dec 2003 11:24:17 -0000	1.11
@@ -31,10 +31,11 @@
 sig_atomic_t loop_active;
 int seekable;
 
-static int input_rate;
+int input_rate;
 int input_ch;
 static int inbytes;
 static int signp;
+int input_size;
 
 typedef struct {
   FILE *f;
@@ -308,7 +309,7 @@
     }
   }
 
-  out.size=2048;
+  input_size=out.size=2048;
   input_ch=out.channels=ch;
   input_rate=out.rate=rate;
   out.data=malloc(sizeof(*out.data)*ch);

<p><p>1.2       +13 -173   postfish/lpc.c

Index: lpc.c
===================================================================
RCS file: /usr/local/cvsroot/postfish/lpc.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -r1.1 -r1.2
--- lpc.c	20 Oct 2003 03:08:43 -0000	1.1
+++ lpc.c	20 Dec 2003 11:24:17 -0000	1.2
@@ -47,6 +47,7 @@
 *********************************************************************/
 
 #define _ISOC99_SOURCE
+#define _GNU_SOURCE
 
 #include <stdlib.h>
 #include <string.h>
@@ -66,23 +67,6 @@
   }
 }
 
-static void autocorr_incomplete(double *data,double *aut,int n,int m){
-  int i,j;
-  
-  /* autocorrelation, p+1 lag coefficients */
-  j=m+1;
-  while(j--){
-    int count=0;
-    double d=0.; /* double needed for accumulator depth */
-    for(i=j;i<n;i++){
-      if(!isnan(data[i]) && !isnan(data[i-j])){
-	d+=data[i]*data[i-j];
-	count++;
-      }
-    }
-    aut[j]=d;
-  }
-}
 
 /* Autocorrelation LPC coeff generation algorithm invented by
    N. Levinson in 1947, modified by J. Durbin in 1959. */
@@ -91,14 +75,14 @@
   double error=aut[0];
   int i,j;
   
-  if(error==0){
-    memset(lpc,0,m*sizeof(*lpc));
-    return 0.;
-  }
-
   for(i=0;i<m;i++){
     double r= -aut[i+1];
 
+    if(error==0){
+      memset(lpc,0,m*sizeof(*lpc));
+      return 0.;
+    }
+
     /* Sum up this iteration's reflection coefficient; note that in
        Vorbis we don't save it.  If anyone wants to recycle this code
        and needs reflection coefficients, save the results of 'r' from
@@ -135,12 +119,6 @@
   return levinson_durbin(aut,lpc,m);
 }
 
-double lpc_from_incomplete_data(double *data,double *lpc,int n,int m){
-  double *aut=alloca(sizeof(*aut)*(m+1));
-  autocorr_incomplete(data,aut,n,m);
-  return levinson_durbin(aut,lpc,m);
-}
-
 void lpc_extrapolate(double *lpc,double *prime,int m,
                      double *data,int n){
   
@@ -165,15 +143,14 @@
     p=m;
     for(j=0;j<m;j++)
       y-=work[o++]*lpc[--p];
-    
-    data[i]=work[o]=y+data[i];
+
+    data[i]=work[o]=y;
   }
 }
 
-void lpc_subtract(double *lpc,double *prime,int m,
-		  double *data,int n){
-  
-  long i,j,o,p;
+void lpc_extrapolateB(double *lpc,double *prime,int m,
+		     double *data,int n){
+    long i,j,o,p;
   double y;
   double *work=alloca(sizeof(*work)*(m+n));
   
@@ -182,7 +159,7 @@
       work[i]=0.;
   else
     for(i=0;i<m;i++)
-      work[i]=prime[i];
+      work[i]=prime[m-i-1];
   
   for(i=0;i<n;i++){
     y=0;
@@ -190,148 +167,11 @@
     p=m;
     for(j=0;j<m;j++)
       y-=work[o++]*lpc[--p];
-    work[o]=data[i];
-    data[i]-=y;
-  }
-}
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <time.h>
-#include <math.h>
-
-#define todB(x)   ((x)==0?-140.f:log((x)*(x))*4.34294480f)
-#define toBARK(n)   (13.1f*atan(.00074f*(n))+2.24f*atan((n)*(n)*1.85e-8f)+1e-4f*(n))
-
-void _analysis(char *base,int i,double *v,int n,int bark,int dB){
-  int j;
-  FILE *of;
-  char buffer[80];
 
-  sprintf(buffer,"%s_%d.m",base,i);
-  of=fopen(buffer,"w");
-  
-  if(!of)perror("failed to open data dump file");
-  
-  for(j=0;j<n;j++){
-    if(bark){
-      float b=toBARK((4000.f*j/n)+.25);
-      fprintf(of,"%f ",b);
-    }else
-      fprintf(of,"%f ",(double)j);
-    
-    if(dB){
-      fprintf(of,"%f\n",todB(v[j]));
-    }else{
-      fprintf(of,"%f\n",v[j]);
-    }
+    data[n-i-1]=work[o]=y;
   }
-  fclose(of);
 }
 
-void lpc_to_curve(double *curve,double *lpc,int n,int m,double amp){
-  int i;
-  drft_lookup fft;
-  double *work=alloca(sizeof(*work)*n*2);
-  memset(work,0,sizeof(double)*n*2);
-
-  for(i=0;i<m;i++){
-    work[i*2+1]=lpc[i]/(4*amp);
-    work[i*2+2]=-lpc[i]/(4*amp);
-  }
-  
-  drft_init(&fft,n*2);
-  drft_backward(&fft,work); /* reappropriated ;-) */
-  drft_clear(&fft);
-
-  {
-    int n2=n*2;
-    double unit=1./amp;
-    curve[0]=(1./(work[0]*2+unit));
-    for(i=1;i<n;i++){
-      double real=(work[i]+work[n2-i]);
-      double imag=(work[i]-work[n2-i]);
-      double a = real + unit;
-      curve[i] = 1.0 / hypot(a, imag);
-    }
-  }
-}  
-
-#define N    1024
-#define M    32
-#define READ (N+M)
-
-signed char readbuffer[READ*4+44];
-
-main(){
-  static int seq;
-  int eos=0,ret;
-  int i, founddata;
-
-  readbuffer[0] = '\0';
-  for (i=0, founddata=0; i<30 && ! feof(stdin) && ! ferror(stdin); i++)
-  {
-    fread(readbuffer,1,2,stdin);
-
-    if ( ! strncmp((char*)readbuffer, "da", 2) )
-    {
-      founddata = 1;
-      fread(readbuffer,1,6,stdin);
-      break;
-    }
-  }
-
-  while(!eos){
-    double      samples[READ];
-    double      curve[N];
-    double      lpc[M];
-
-    long i;
-    long bytes=fread(readbuffer,1,READ*4,stdin); /* stereo hardwired here */
-    
-    if(bytes<READ*4)break;
-      
-    /* uninterleave samples */
-    for(i=0;i<bytes/4;i++)
-      samples[i]=
-	((readbuffer[i*4+1]<<8)|
-	 (0x00ff&(int)readbuffer[i*4]))/32768. +
-	((readbuffer[i*4+3]<<8)|
-	 (0x00ff&(int)readbuffer[i*4+2]))/32768.f;
-    
-    
-    _analysis("pcm",seq,samples,READ,0,0);
-    lpc_from_data(samples,lpc,READ,M);
-    _analysis("lpc",seq,lpc,M,0,0);
-    lpc_to_curve(curve,lpc,N,M,1.0);
-    _analysis("H",seq,curve,N,1,1);
-
-
-    /* kill off samples */
-    {
-      double max=0;
-      for(i=0;i<READ;i++)
-	if(samples[i]>max)max=samples[i];
-      for(i=0;i<READ;i++){
-	if(samples[i]>max*.50)samples[i]=max*.50;
-	if(samples[i]<-max*.50)samples[i]=-max*.50;
-      }
-    }
-    _analysis("pcmi",seq,samples,READ,0,0);
-    lpc_from_incomplete_data(samples,lpc,READ,M);
-    _analysis("lpci",seq,lpc,M,0,0);
-    lpc_to_curve(curve,lpc,N,M,1.0);
-    _analysis("Hi",seq,curve,N,1,1);
-    lpc_subtract(lpc,samples,M,samples+M,N);
-    _analysis("resi",seq,samples,READ,0,0);
-
 
 
-    seq++;
-  }
-  
-
-  return 0;
-}
 

<p><p>1.23      +2 -2      postfish/mainpanel.c

Index: mainpanel.c
===================================================================
RCS file: /usr/local/cvsroot/postfish/mainpanel.c,v
retrieving revision 1.22
retrieving revision 1.23
diff -u -r1.22 -r1.23
--- mainpanel.c	18 Oct 2003 11:10:09 -0000	1.22
+++ mainpanel.c	20 Dec 2003 11:24:17 -0000	1.23
@@ -900,9 +900,9 @@
   int i;
 
   memset(&p,0,sizeof(p));
-  gtk_rc_add_default_file("/etc/postfish/postfishrc");
+  gtk_rc_add_default_file("/etc/postfish/postfish-gtkrc");
   if(homedir){
-    char *rcfile="/.postfishrc";
+    char *rcfile="/.postfish-gtkrc";
     char *homerc=calloc(1,strlen(homedir)+strlen(rcfile)+1);
     strcat(homerc,homedir);
     strcat(homerc,rcfile);

<p><p>1.2       +253 -106  postfish/minimize.c

Index: minimize.c
===================================================================
RCS file: /usr/local/cvsroot/postfish/minimize.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -r1.1 -r1.2
--- minimize.c	7 Nov 2003 17:10:40 -0000	1.1
+++ minimize.c	20 Dec 2003 11:24:17 -0000	1.2
@@ -20,123 +20,270 @@
  *
  * 
  */
+#define _GNU_SOURCE
+#include <stdio.h>
+#include "smallft.h"
+
+#include "test.h"
+
+#include <stdlib.h>
+#include <math.h>
+#define toBark(n)   (13.1f*atan(.00074f*(n))+2.24f*atan((n)*(n)*1.85e-8f)+1e-4f*(n))
+#define fromBark(z) (102.f*(z)-2.f*pow(z,2.f)+.4f*pow(z,3.f)+pow(1.46f,z)-1.f)
+
+static void sliding_bark_average(double *f,double *w, int n,double width){
+  int i=0,j;
+  int lopad=1-rint(fromBark(toBark(0.)-width)*n/44100.);
+  int hipad=rint(fromBark(toBark(22050.)+width)*n/44100.)+lopad;
+  double acc=0.,del=0.;
+  double sec[hipad+1];
+
+  memset(sec,0,sizeof(sec));
+
+  {
+    double bark=toBark(0.);
+    int hi=rint(fromBark(bark-width)*n/44100.)-1;
+    int lo=rint(fromBark(bark+width)*n/44100.)+1;
+    double del=fabs(f[0])/(lo-hi);
+
+    double hidel=del/(-hi);
+    double lodel=del/(lo);
+
+    sec[hi+lopad]+=hidel;
+    sec[lopad]-=hidel;
+    sec[lopad]-=lodel;
+    sec[lo+lopad]+=lodel;
+    
+  }
+
+  for(i=1;i<n/2;i++){
+    double bark=toBark(44100.*i/n);
+    int hi=rint(fromBark(bark-width)*n/44100.)-1;
+    int lo=rint(fromBark(bark+width)*n/44100.)+1;
+    double del=hypot(f[(i<<1)-1],f[i<<1])/(lo-hi);
+
+    double hidel=del/((i-hi));
+    double lodel=del/((lo-i));
+
+    sec[hi+lopad]+=hidel;
+    sec[i+lopad]-=hidel;
+    sec[i+lopad]-=lodel;
+    sec[lo+lopad]+=lodel;
 
-/* Implment a linear system least-squares minimizer.  Typical
-   implementations are inefficient due to our need for a weighting
-   vector that changes from frame to frame, eliminating our ability to
-   presolve the first- and second-order differential matricies.  This
-   alternate to Newton's Method allows fast worst-case O(N^2) direct
-   solution of the least-squares fit. */
-
-/* if no 'x' vector, assume zeroes.  If no slope vector, assume unit slope */
-static double line_minimize(double **A, double *x,double *slope,
-			    double *c,double *w,
-			    int *flag,int n, double *grad){
-  double c2[n];
-  double A2[n];
-  double f1=0.,f2=0.,delta;
-  int i,j;
-
-  for(i=0;i<n;i++){
-    c2[i]=c[i];
-    A2[i]=0;
-    if(x && flag[i]){
-      double *Ai=A[i];
-      if(slope){
-	for(j=0;j<n;j++)
-	  if(flag[j]){
-	    A2[i]+=Ai[j]*s[j];
-	    c2[i]+=Ai[j]*x[j];
-	  }
-      }else
-	for(j=0;j<n;j++)
-	  if(flag[j]){
-	    A2[i]+=Ai[j];
-	    c2[i]+=Ai[j]*x[j];
-	  }
-    }
   }
-  
-  /* calculate first derivative */
-  for(i=0;i<n;i++)if(flag[i]) f1+= 2. * A2[i] * c2[i] * w[i];
 
-  /* calculate second derivative */
-  for(i=0;i<n;i++)if(flag[i]) f2+= 2. * A2[i] * A2[i] * w[i];
+  for(i=0;i<lopad;i++){
+    del+=sec[i];
+    acc+=del;
+  }
 
-  /* thus our minimum is at... */
-  delta = -f1/f2;
+  w[0]=acc;
+  del+=sec[lopad];
+  acc+=del;
+
+  for(i=1;i<n/2;i++){
+    w[(i<<1)-1]=w[i<<1]=acc;
+    del+=sec[i+lopad];
+    acc+=del;
 
-  /* do we need a new grad? */
-  if(grad){
-    for(i=0;i<n;i++)c2[i]+=A2[i]*delta;
-    for(j=0;j<n;j++){
-      grad[j]=0.;
-      if(flag[j])
-	for(i=0;i<n;i++)
-	  grad[j]+= 2. * c2[i] * A[i][j] * w[i];
-    }
   }
+  w[n-1]=w[n-2];
+}
+
+#define todB(x)   ((x)==0?-140.f:log((x)*(x))*4.34294480f)
+#define toBARK(n)   (13.1f*atan(.00074f*(n))+2.24f*atan((n)*(n)*1.85e-8f)+1e-4f*(n))
 
-  return delta;
+void _analysis(char *base,int i,double *v,int n,int bark,int dB){
+  int j;
+  FILE *of;
+  char buffer[80];
+
+  sprintf(buffer,"%s_%d.m",base,i);
+  of=fopen(buffer,"w");
+  
+  if(!of)perror("failed to open data dump file");
+  
+  for(j=0;j<n;j++){
+    if(bark){
+      float b=toBARK((4000.f*j/n)+.25);
+      fprintf(of,"%f ",b);
+    }else
+      fprintf(of,"%f ",(double)j);
+    
+    if(dB){
+      if(j==0||j==n-1)
+	fprintf(of,"%f\n",todB(v[j]));
+      else{
+	fprintf(of,"%f\n",todB(hypot(v[j],v[j+1])));
+	j++;
+      }
+    }else{
+      fprintf(of,"%f\n",v[j]);
+    }
+  }
+  fclose(of);
 }
 
-double Nspace_minimum(double **A, double *x, double *w, int *flag,
-		      double *outgrad,int n){
-  int i,j;
-  double constants[n];
-  double delta;
-  double delta2;
-  double work[n];
-  double slope[n];
-  double grad[n];
-
-  /* everything we're concerned with is in quadratic form, so Newton's
-     method will solve minimization directly.  We don't do the
-     minimization in n-space as computation of the second derivative
-     of our full cost function is N^3 (due to the weighting vector).
-     However, we can perform the same operation with three 1-d
-     minimizations in n-space */
-
-  /* begin by assuming the only free variables we have are the x's
-     marked by flag.  The fixed x's immediately collapse out of the
-     original cost function and into constants. */
-  for(i=0;i<n;i++){
-    constants[i]=0;
-    if(flag[i]){
-      double *Ai=A[i];
-      for(j=0;j<n;j++)
-	if(!flag[j])constants[i]+=Ai[j]*x[j];
+#define N    128
+
+extern void minimize(drft_lookup *fft,double **A, double *w, double *x, 
+		    int *flag, 
+		    int piecep, int n1, int n2, int n);
+
+
+int main(){
+  drft_lookup fft;
+  
+  double **A;
+  double x[N*2];
+  double s[N*2];
+  double f[N*2];
+  double w[N*2];
+  double wi[N*2];
+  int flag[N*2];
+
+  double lap[N];
+  int i,j,freex=0,eos=0;
+  static int seq=0;
+  signed char readbuffer[N*4+44];
+  long founddata;
+
+  memset(lap,0,sizeof(lap));
+  memset(x,0,sizeof(x));
+  /* set up the A matrix needed to compute gradient */
+  drft_init(&fft,N*2);
+  A=malloc(N*2*sizeof(*A));
+  for(i=0;i<N*2;i++)
+    A[i]=calloc(N*2,sizeof(**A));
+  
+  for(i=0;i<N*2;i++){
+    A[i][i]=1.;
+    drft_forward(&fft,A[i]);
+  }
+
+  for(j=0;j<N*2;j++) wi[j]=sin( (double)j/N*M_PIl*.5 );
+  for(j=0;j<N*2;j++) wi[j]*=wi[j];
+  for(j=0;j<N*2;j++) wi[j]=sin(wi[j]*M_PIl*.5 );
+
+  readbuffer[0] = '\0';
+  for (i=0, founddata=0; i<30 && ! feof(stdin) && ! ferror(stdin); i++){
+    fread(readbuffer,1,2,stdin);
+    fwrite(readbuffer,1,2,stdout);
+    
+    if ( ! strncmp((char*)readbuffer, "da", 2) ){
+      founddata = 1;
+      fread(readbuffer,1,6,stdin);
+      fwrite(readbuffer,1,6,stdout);
+      break;
     }
   }
 
-  /* initially use an arbitrary slope to define a line passing
-     through the origin of the n-space defined by the free variables;
-     a unit slope is convenient. */
+
+  while(!eos){
+    long bytes=fread(readbuffer,1,N*4,stdin); /* stereo hardwired here */
+    if(bytes<N*4)break;
+
+    freex=0;
+
+    /* uninterleave samples */
+    for(i=0;i<bytes/4;i++){
+      x[i+N]=
+        ((readbuffer[i*4+1]<<8)|
+         (0x00ff&(int)readbuffer[i*4]))/65536. +
+        ((readbuffer[i*4+3]<<8)|
+         (0x00ff&(int)readbuffer[i*4+2]))/65536.;
+      flag[i+N]=0;
+    }
+
+    for(i=0;i<N;i++){
+      flag[i+N]=0;
+      if(x[i+N]>=.22){
+	flag[i+N]=1;
+	//x[i+N]=.2;
+	freex++;
+      }else if(x[i+N]<=-.2){
+	flag[i+N]=-1;
+	//x[i+N]=-.2;
+	freex++;
+      }
+    }
+    
+    memcpy(f,x,sizeof(f));
   
-  /* Fix the free variables to zero; rewrite the cost function as a
-     function of a single variable 'delta' that scales the unit slope.
-     Minimize the rewritten cost function for delta. Also get the grad
-     at the minimum as the new slope (which falls out at low cost as a
-     normal vector, allowing easy construction of a parallel line for
-     the second step) */
-  delta=line_minimize(A,NULL,NULL,constants,w,flag,n,grad);
-
-  /* offset the origin by the gradient, and set up a new line
-     minimization */
-  delta2=delta-line_minimize(A,grad,NULL,constants,w,flag,n,NULL);
-   
-  /* construct a line passing through the two minima */
-  /* faster not to if()-out the unused elements */
-  for(i=0;i<n;i++){
-    work[i]=delta;
-    slope[i]=grad[i]+delta2;
-  }
-
-  /* one last minimization.  We might need the gradient */
-  delta2=line_minimize(A,work,slope,constants,w,flag,n,outgrad);
-
-  /* compute new x's */
-  for(i=0;i<n;i++)
-    if(flag[i])
-      x[i]=work[i]+slope[i]*delta2;
+    /* weight the bugger */
+    for(i=0;i<N*2;i++)f[i]*=wi[i];
+    memcpy(s,f,sizeof(s));
+    
+    if(freex){
+      fprintf(stderr,"%d: Free variables: %d\n",seq,freex);
+
+      drft_forward(&fft,f);
+      _analysis("pcm",seq,s,N*2,0,0);
+      _analysis("f",seq,f,N*2,0,1);
+      sliding_bark_average(f,w,N*2,.5);
+      //memcpy(w,f,sizeof(w));
+
+      //for(i=1;i<N*2;i+=2)w[i]=w[i+1]=hypot(f[i],f[i+1]);
+
+
+      _analysis("w",seq,w,N*2,0,1);
+
+      for(i=0;i<N*2;i++)w[i]=1./w[i];
+
+      if(flag[0] && flag[N*2-1]){
+	flag[0]=0;
+	s[0]=0;
+      }
+      minimize(&fft,A,w,s,flag,0,0,N*2,N*2);
+
+      memcpy(f,s,sizeof(f));
+      _analysis("x",seq,f,N*2,0,0);
+      drft_forward(&fft,f);
+      _analysis("xf",seq++,f,N*2,0,1);
+      
+#if 0
+      memcpy(f,s,sizeof(s));
+      memcpy(s,x,sizeof(s));
+      for(i=0;i<N*2;i++)s[i]*=wi[i];
+      
+      drft_forward(&fft,f);
+      sliding_bark_average(f,w,N*2,.5);
+      _analysis("w2",seq,w,N*2,0,1);
+      for(i=0;i<N*2;i++)w[i]=1./w[i];
+      minimize(&fft,A,w,s,flag,0,0,N*2,N*2);
+      
+      memcpy(f,s,sizeof(f));
+      _analysis("x2",seq,f,N*2,0,0);
+      drft_forward(&fft,f);
+      _analysis("xf2",seq++,f,N*2,0,1);
+#endif
+    }
+
+    /* lap */
+    for(i=0;i<N;i++)lap[i]+=s[i]*wi[i];
+
+    for(i=0;i<bytes/4;i++){
+      int val=rint(lap[i]*32768.);
+      if(val>32767)val=32767;
+      if(val<-32768)val=-32768;
+      readbuffer[i*4]=val&0xff;
+      readbuffer[i*4+1]=(val>>8)&0xff;
+      readbuffer[i*4+2]=val&0xff;
+      readbuffer[i*4+3]=(val>>8)&0xff;
+    }
+    fwrite(readbuffer,1,bytes,stdout);
+
+    for(i=0;i<N;i++)lap[i]=s[i+N]*wi[i+N];
+    
+    memmove(x,x+N,sizeof(x)-N*sizeof(*x));
+    memmove(flag,flag+N,sizeof(flag)-N*sizeof(*flag));
+  }
+  return 0;
 }
+
+
+
+
+
+
+

<p><p>1.12      +2 -2      postfish/multibar.c

Index: multibar.c
===================================================================
RCS file: /usr/local/cvsroot/postfish/multibar.c,v
retrieving revision 1.11
retrieving revision 1.12
diff -u -r1.11 -r1.12
--- multibar.c	18 Oct 2003 07:29:47 -0000	1.11
+++ multibar.c	20 Dec 2003 11:24:17 -0000	1.12
@@ -272,8 +272,8 @@
         gdk_draw_line(m->backing,widget->style->fg_gc[0],
                       x-1,y-3,x+1,y-3);
       
-	gdk_draw_line(m->backing,widget->style->fg_gc[0],
-		      x,0,x,y-1);
+	gdk_draw_line(m->backing,widget->style->fg_gc[1],
+		      x,1,x,y-2);
       }
     }
   }

<p><p>1.9       +3 -2      postfish/postfish.h

Index: postfish.h
===================================================================
RCS file: /usr/local/cvsroot/postfish/postfish.h,v
retrieving revision 1.8
retrieving revision 1.9
diff -u -r1.8 -r1.9
--- postfish.h	18 Oct 2003 11:10:09 -0000	1.8
+++ postfish.h	20 Dec 2003 11:24:17 -0000	1.9
@@ -30,8 +30,7 @@
 #define _POSTFISH_H_
 
 #define _GNU_SOURCE
-#define _LARGEFILE_SOURCE 
-#define _LARGEFILE64_SOURCE
+#define _ISOC99_SOURCE
 #define _FILE_OFFSET_BITS 64
 #define _REENTRANT 1
 #include <stdlib.h>
@@ -54,6 +53,8 @@
 #define todB(x)   ((x)==0?-400.f:log((x)*(x))*4.34294480f)
 #define fromdB(x) (exp((x)*.11512925f))  
 #define toOC(n)     (log(n)*1.442695f-5.965784f)
+#define toBark(n)   (13.1f*atan(.00074f*(n))+2.24f*atan((n)*(n)*1.85e-8f)+1e-4f*(n))
+#define fromBark(z) (102.f*(z)-2.f*pow(z,2.f)+.4f*pow(z,3.f)+pow(1.46f,z)-1.f)
 
 typedef struct time_linkage {
   int size;

<p><p>1.3       +16 -29    postfish/reconstruct.c

Index: reconstruct.c
===================================================================
RCS file: /usr/local/cvsroot/postfish/reconstruct.c,v
retrieving revision 1.2
retrieving revision 1.3
diff -u -r1.2 -r1.3
--- reconstruct.c	13 Dec 2003 12:26:55 -0000	1.2
+++ reconstruct.c	20 Dec 2003 11:24:17 -0000	1.3
@@ -31,28 +31,16 @@
 #include "smallft.h"
 #include "reconstruct.h"
 
-/* this setup isn't thread safe, but it's fine for postfish which will
-   only access it from a single thread */
-
-drft_lookup fft;
-void set_up_filter(int n){
-  static int cached_n=0;
-  if(n!=cached_n){
-    if(n)drft_clear(&fft);
-    drft_init(&fft,n);
-    cached_n=n;
-  }
-}
-
-double inner_product(double *a, double *b, int n){
+static double inner_product(double *a, double *b, int n){
   int i;
   double acc=0.;
   for(i=0;i<n;i++)acc+=a[i]*b[i];
   return acc;
 }
 
-void compute_AtAx(double *x,double *w,int *flag,int mask,int n,
-		    double *out){
+static void compute_AtAx(drft_lookup *fft,
+			 double *x,double *w,int *flag,int mask,int n,
+			 double *out){
   int i;
 
   if(mask){
@@ -68,24 +56,25 @@
       else
         out[i]=x[i];
   
-  drft_forward(&fft,out);
+  drft_forward(fft,out);
   for(i=0;i<n;i++)out[i]*=w[i];
-  drft_backward(&fft,out);
+  drft_backward(fft,out);
 
   for(i=0;i<n;i++)
     if(!flag[i])out[i]=0;
   
 }
 
-void compute_Atb_minus_AtAx(double *x,double *w,double *Atb,int *flag,int n,
-			    double *out){
+static void compute_Atb_minus_AtAx(drft_lookup *fft,
+				   double *x,double *w,double *Atb,int *flag,
+				   int n,double *out){
   int i;
-  compute_AtAx(x,w,flag,1,n,out);
+  compute_AtAx(fft,x,w,flag,1,n,out);
   for(i=0;i<n;i++)out[i]=Atb[i]-out[i];
 }
 
-
-void reconstruct(double *x, double *w, int *flag, double e,int max,int n){
+void reconstruct(drft_lookup *fft,
+		 double *x, double *w, int *flag, double e,int max,int n){
   int i,j;
   double Atb[n];
   double r[n];
@@ -94,23 +83,21 @@
   double phi_new,phi_old,phi_0;
   double alpha,beta;
 
-  set_up_filter(n);
-
   /* compute initial Atb */
-  compute_AtAx(x,w,flag,0,n,Atb);
+  compute_AtAx(fft,x,w,flag,0,n,Atb);
   for(j=0;j<n;j++)Atb[j]= -Atb[j];
 
-  compute_Atb_minus_AtAx(x,w,Atb,flag,n,r);
+  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);
 
   for(i=0;i<max && phi_new>e*e*phi_0;i++){
-    compute_AtAx(d,w,flag,1,n,q);
+    compute_AtAx(fft,d,w,flag,1,n,q);
     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(x,w,Atb,flag,n,r);
+      compute_Atb_minus_AtAx(fft,x,w,Atb,flag,n,r);
     else
       for(j=0;j<n;j++)r[j]-=alpha*q[j];
     

<p><p>1.2       +1 -1      postfish/reconstruct.h

Index: reconstruct.h
===================================================================
RCS file: /usr/local/cvsroot/postfish/reconstruct.h,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -r1.1 -r1.2
--- reconstruct.h	13 Dec 2003 12:02:11 -0000	1.1
+++ reconstruct.h	20 Dec 2003 11:24:17 -0000	1.2
@@ -21,5 +21,5 @@
  * 
  */
 
-extern void reconstruct(double *x, double *w, int *flag, 
+extern void reconstruct(drft_lookup *fft,double *x, double *w, int *flag, 
                         double e,int max,int n);

<p><p>1.21      +2 -2      postfish/version.h

Index: version.h
===================================================================
RCS file: /usr/local/cvsroot/postfish/version.h,v
retrieving revision 1.20
retrieving revision 1.21
diff -u -r1.20 -r1.21
--- version.h	18 Oct 2003 11:12:43 -0000	1.20
+++ version.h	20 Dec 2003 11:24:17 -0000	1.21
@@ -1,2 +1,2 @@
-#define VERSION "$Id: version.h,v 1.20 2003/10/18 11:12:43 xiphmont Exp $ "
-/* DO NOT EDIT: Automated versioning hack [Sat Oct 18 07:11:13 EDT 2003] */
+#define VERSION "$Id: version.h,v 1.21 2003/12/20 11:24:17 xiphmont Exp $ "
+/* DO NOT EDIT: Automated versioning hack [Sat Dec 20 05:44:23 EST 2003] */

<p><p>1.1                  postfish/declip.c

Index: declip.c
===================================================================
/*
 *
 *  postfish
 *    
 *      Copyright (C) 2002-2004 Monty and Xiph.Org
 *
 *  Postfish is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2, or (at your option)
 *  any later version.
 *   
 *  Postfish is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *   
 *  You should have received a copy of the GNU General Public License
 *  along with Postfish; see the file COPYING.  If not, write to the
 *  Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 * 
 */

#include "postfish.h"
#include <math.h>
#include <sys/types.h>
#include "smallft.h"
#include "reconstruct.h"
#include <stdio.h>

extern int input_rate;
extern int input_ch;
extern int input_size;

void _analysis(char *base,int i,double *v,int n,int bark,int dB){
  int j;
  FILE *of;
  char buffer[80];

  sprintf(buffer,"%s_%d.m",base,i);
  of=fopen(buffer,"w");
  
  if(!of)perror("failed to open data dump file");
  
  for(j=0;j<n;j++){
    if(bark){
      float b=toBark((4000.f*j/n)+.25);
      fprintf(of,"%f ",b);
    }else
      fprintf(of,"%f ",(double)j);
    
    if(dB){
      if(j==0||j==n-1)
        fprintf(of,"%f\n",todB(v[j]));
      else{
        fprintf(of,"%f\n",todB(hypot(v[j],v[j+1])));
        j++;
      }
    }else{
      fprintf(of,"%f\n",v[j]);
    }
  }
  fclose(of);
}

/* accessed only in playback thread/setup */
static drft_lookup fft;
static int blocksize=0;
static int lopad=0,hipad=0;
static u_int32_t *widthlookup=0;
static double *window=0;
static double width=.5;
static double **lap=0;
static double **cache;
static int cache_samples;
static int fillstate=0; /* 0: uninitialized
                           1: normal
                           2: eof processed */

tatic time_linkage out;

/* accessed across threads */
static sig_atomic_t declip_active=0;
static sig_atomic_t declip_converge=2; /* 0=over, 1=full, 2=half, 3=partial, 4=approx */
static sig_atomic_t *chtrigger=0;
static sig_atomic_t *chactive=0;
static sig_atomic_t pending_blocksize=0;

/* called only by initial setup */
int declip_load(void){
  int i;
  chtrigger=malloc(input_ch*sizeof(*chtrigger));
  chactive=malloc(input_ch*sizeof(*chactive));
  for(i=0;i<input_ch;i++){
    chtrigger[i]=0x80000000UL;
    chactive[i]=1;
  }
  
  out.size=input_size;
  out.channels=input_ch;
  out.rate=input_rate;
  out.data=malloc(input_ch*sizeof(*out.data));
  for(i=0;i<input_ch;i++)
    out.data[i]=malloc(input_size*sizeof(**out.data));

  fillstate=0;
  cache=malloc(input_ch*sizeof(*cache));
  for(i=0;i<input_ch;i++)
    cache[i]=malloc(input_size*sizeof(**cache));
  lap=malloc(input_ch*sizeof(*lap));
  for(i=0;i<input_ch;i++)
    lap[i]=malloc(input_size*sizeof(**lap));

  return(0);
}

int declip_setblock(int n){
  if(n<32)return -1;
  if(n>input_size*2)return -1;
  pending_blocksize=n;
  return 0;
}

int declip_setactive(int activep,int ch){
  if(ch<0 || ch>=input_ch)return -1;
  chactive[ch]=activep;
  return 0;
}

int declip_settrigger(double trigger,int ch){
  if(ch<0 || ch>=input_ch)return -1;
  chtrigger[ch]=rint(trigger*(1.*0x80000000UL));
  return 0;
}

/* called only in playback thread */
int declip_reset(void){
  /* reset cached pipe state */
  fillstate=0;
  return 0;
}

tatic void sliding_bark_average(double *f,double *w, int n,double width){
  int i=0,j;
  double acc=0.,del=0.;
  double sec[hipad+1];

  memset(sec,0,sizeof(sec));

  {
    double bark=toBark(0.);
    int hi=widthlookup[0]>>16;
    int lo=widthlookup[0]&(0xffff);
    double del=fabs(f[0])/(lo-hi);

    double hidel=del/(-hi+lopad);
    double lodel=del/(lo-lopad);

    sec[hi]+=hidel;
    sec[lopad]-=hidel;
    sec[lopad]-=lodel;
    sec[lo]+=lodel;
    
  }

  for(i=1;i<n/2;i++){

    double bark=toBark(44100.*i/n);
    int hi=widthlookup[i]>>16;
    int lo=widthlookup[i]&(0xffff);
    double del=hypot(f[(i<<1)-1],f[i<<1])/(lo-hi);

    double hidel=del/((i-hi+lopad));
    double lodel=del/((lo-i-lopad));

    sec[hi]+=hidel;
    sec[i+lopad]-=hidel;
    sec[i+lopad]-=lodel;
    sec[lo]+=lodel;

  }

  for(i=0;i<lopad;i++){
    del+=sec[i];
    acc+=del;
  }

  w[0]=1./(acc*acc);
  del+=sec[lopad];
  acc+=del;

  for(i=1;i<n/2;i++){
    w[(i<<1)-1]=w[i<<1]=1./(acc*acc);
    del+=sec[i+lopad];
    acc+=del;

  }
  w[n-1]=w[n-2];
}

tatic void declip(double *data,double *lap,double *out,
                   int blocksize,unsigned long trigger){
  double freq[blocksize];
  int    flag[blocksize];
  double triggerlevel=trigger*(1./0x80000000UL);
  double epsilon=1e-12;
  int    iterbound=blocksize,i,j,count=0;
  
  for(i=0;i<blocksize/8;i++)flag[i]=0;
  for(;i<blocksize*7/8;i++){
    flag[i]=0;
    if(data[i]>=trigger || data[i]<=-trigger){
      flag[i]=1;
      count++;
    }
  }
  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);
  
  switch(declip_converge){
  case 0:
    epsilon=1e-12;
    iterbound=blocksize*2;
    break;
  case 1:
    epsilon=1e-8;
    iterbound=count;
    break;
  case 2:
    epsilon=1e-6;
    iterbound=count/2;
    break;
  case 3:
    epsilon=1e-5;
    iterbound=count/4;
    break;
  case 4:
    epsilon=1e-3;
    iterbound=count/8;
    break;
  }
  if(iterbound<20)iterbound=20;

  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];

  for(i=blocksize/2,j=0;i<blocksize;i++)
    lap[j]=data[i]*window[i];
}

/* called only by playback thread */
time_linkage *declip_read(time_linkage *in){
  int i;
  double work[blocksize];

  if(pending_blocksize!=blocksize){
    if(blocksize){
      free(widthlookup);
      free(window);
      drft_clear(&fft);
    }
    blocksize=pending_blocksize;

    widthlookup=malloc((blocksize>>1)*sizeof(*widthlookup));
    for(i=0;i<blocksize/2;i++){
      double bark=toBark(input_rate*i/blocksize);
      int hi=rint(fromBark(bark-width)*blocksize/input_rate)-1+lopad;
      int lo=rint(fromBark(bark+width)*blocksize/input_rate)+1+lopad;
      
      if(hi<0 || lo<0 || hi>65535 || lo<65535) return 0;
      widthlookup[i]=(hi<<16)+lo;
    }
    lopad=1-rint(fromBark(toBark(0.)-width)*blocksize/input_rate);
    hipad=rint(fromBark(toBark(input_rate*.5)+width)*blocksize/input_rate)+lopad;
    
    drft_init(&fft,blocksize);

    window=malloc(blocksize*sizeof(*window));
    for(i=0;i<blocksize/8;i++) window[i]=0.;
    for(;i<blocksize*3/8;i++) window[i]=sin( (double)(i-blocksize/8)/blocksize*M_PIl );
    for(;i<blocksize*5/8;i++) window[i]=1.;
    for(;i<blocksize*7/8;i++) window[i]=sin( (double)(blocksize*7/8-i)/blocksize*M_PIl );
    for(;i<blocksize;i++) window[i]=0.;    
    for(i=0;i<blocksize;i++) window[i]*=window[i];
    for(i=0;i<blocksize;i++) window[i]=sin(window[i]*M_PIl);

  }

  switch(fillstate){
  case 0: /* prime the lapping and cache */
    for(i=0;i<input_ch;i++){
      int j;
      double *temp=in->data[i];
      if(chactive[i] && declip_active){
        memset(work,0,sizeof(*work)*blocksize/2);
        memcpy(work+blocksize/2,temp,sizeof(*work)*blocksize/2);
        declip(work,lap[i],0,blocksize,chtrigger[i]);
      }else{
        for(j=0;j<blocksize/2;j++)
          lap[i][j]=window[j+blocksize/2]*temp[j];
      }
      memset(cache[i],0,sizeof(**cache)*input_size);
      in->data[i]=cache[i];
      cache[i]=temp;
    }
    cache_samples=in->samples;
    fillstate=1;
    if(in->samples==in->size)return 0;
    in->samples=0;
    /* fall through */
  case 1: /* nominal processing */
    for(i=0;i<input_ch;i++){
      double *temp=cache[i];
        int j;
      if(chactive[i] && declip_active){
        for(j=0;j+blocksize<out.size;j+=blocksize/2){
          memcpy(work,temp+j,sizeof(*work)*blocksize);
          declip(work,lap[i],out.data[i]+j,blocksize,chtrigger[i]);
        }
        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,chtrigger[i]);
      }else{
        memcpy(out.data[i],temp,out.size*sizeof(*temp));
        for(j=0;j<blocksize/2;j++)
          lap[i][j]=window[j+blocksize/2]*temp[j];
      }
      
      in->data[i]=cache[i];
      cache[i]=temp;
    }
    out.samples=cache_samples;
    cache_samples=in->samples;
    if(out.samples<out.size)fillstate=2;
    break;
  case 2: /* we've pushed out EOF already */
    return 0;
  }
  return &out;
}

<p><p>1.1                  postfish/declip.h

Index: declip.h
===================================================================
/*
 *
 *  postfish
 *    
 *      Copyright (C) 2002-2004 Monty and Xiph.Org
 *
 *  Postfish is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2, or (at your option)
 *  any later version.
 *   
 *  Postfish is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *   
 *  You should have received a copy of the GNU General Public License
 *  along with Postfish; see the file COPYING.  If not, write to the
 *  Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 * 
 */

extern int declip_load(void);
extern int declip_setblock(int n);
extern int declip_setactive(int activep,int ch);
extern int declip_settrigger(double trigger,int ch);
extern int declip_reset(void);
extern time_linkage *declip_read(time_linkage *in);

<p><p>1.1                  postfish/postfish-gtkrc

Index: postfish-gtkrc
===================================================================
style "button-poppy" {
        bg[NORMAL]="#80a0ff" 
        bg[ACTIVE]="#c0f0ff" 
        bg[PRELIGHT]="#c0f0ff" 

        font_name = "sans 8"	
        GtkButton::focus-padding = 0
        GtkButton::focus-line-width = 1
        GtkButton::interior-focus = 0
}

tyle "check-poppy" {
        bg[NORMAL]="#80a0ff" 

        font_name = "sans 8"	
        GtkButton::focus-padding = 0
        GtkButton::focus-line-width = 1
        GtkButton::interior-focus = 0
}

tyle "slider" {
        bg[NORMAL]="#80a0ff" 
        GtkWidget::focus-padding = 0
        GtkWidget::focus-line-width = 1
        GtkWidget::interior-focus = 0
}

tyle "multibar" {
        fg[NORMAL]="#404040" 
        fg[ACTIVE]="#ff8080" 
        text[NORMAL]="#c0c0d0" 
        text[ACTIVE]="#ffb0b0" 
        font_name = "sans 8"	
}

tyle "readout" {
        base[NORMAL]="#ffffff" 
        base[ACTIVE]="#ffffff" 
        bg[NORMAL]="#ffffff" 
        bg[ACTIVE]="#ffffff" 

        font_name = "fixed"	
        GtkMisc::xpad = 10
        GtkMisc::xalign = 1.0
        GtkWidget::width-request=100
        GtkLabel::justify = right
}

tyle "darkpanel" {
        bg[NORMAL]="#b0b0b0" 
}

tyle "quitbutton" {
        bg[NORMAL]="#d0d0d0"
        bg[PRELIGHT]="#ffc0c0"
        bg[ACTIVE]="#ffc0c0"
        font_name = "sans 8"	
        GtkButton::focus-padding = 0
        GtkButton::focus-line-width = 1
        GtkButton::interior-focus = 0
}

tyle "left" {
        text[NORMAL] = "#606060"
        text[ACTIVE] = "#606060"
        text[SELECTED] = "#606060"
        text[PRELIGHT] = "#606060"
        fg[ACTIVE] = "#606060"
        bg[NORMAL]="#80a0ff" 
}
style "right" {
        text[NORMAL] = "#cc0000"
        text[ACTIVE] = "#cc0000"
        text[SELECTED] = "#cc0000"
        text[PRELIGHT] = "#cc0000"
        bg[NORMAL]="#80a0ff" 
}
style "mid" {
        text[NORMAL] = "#0000fc"
        text[ACTIVE] = "#0000fc"
        text[SELECTED] = "#0000fc"
        text[PRELIGHT] = "#0000fc"
        bg[NORMAL]="#80a0ff" 
}
style "side" {
        text[NORMAL] = "#00B200"
        text[ACTIVE] = "#00B200"
        text[SELECTED] = "#00B200"
        text[PRELIGHT] = "#00B200"
        bg[NORMAL]="#80a0ff" 
}

widget "*.color0" style "left"
widget "*.color1" style "right"
widget "*.color2" style "mid"
widget "*.color3" style "side"

widget "*.winpanel" style "darkpanel"

widget "*.Readout*" style "readout"
widget "*.GtkEntry" style "readout"
widget "*.GtkHScale" style "slider"
widget "*.GtkToggleButton.*" style "button-poppy"
widget "*.GtkToggleButton" style "button-poppy"
widget "*.GtkButton.*" style "button-poppy"
widget "*.reset*" style "quitbutton"
widget "*.GtkButton" style "button-poppy"
widget "*.GtkCheckButton" style "check-poppy"
widget "*.quitbutton" style "quitbutton"
widget "*.quitbutton.GtkLabel" style "quitbutton"
widget "*.Multibar*" style "multibar"

<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