Update to 16x32 pre-/post- filter.
authorTimothy B. Terriberry <tterribe@xiph.org>
Fri, 27 Jul 2012 21:23:46 +0000 (14:23 -0700)
committerTimothy B. Terriberry <tterribe@xiph.org>
Fri, 27 Jul 2012 21:23:46 +0000 (14:23 -0700)
This converts it to the Type IV dyadic rational filter with no ramp
 constraint that negge's optimization routines found.
See http://wiki.xiph.org/TDLT

src/filter.c

index db47125..8c6d131 100644 (file)
    acceptable solution found so far.
   The maximum denominator for all coefficients was allowed to be 64.*/
 
-#undef TEST
-
 #if defined(TEST)
 extern int minv[32];
 extern int maxv[32];
@@ -520,14 +518,6 @@ void od_post_filter8(od_coeff _x[8],const od_coeff _y[8]){
 
 void od_pre_filter16(od_coeff _y[16],const od_coeff _x[16]){
    int t[16];
-   /*
-"Found: "{9.788267439364699, {3/2, 39/32, 77/64, 39/32, 19/16, 87/64, 35/32,`
-    1}, {29/32, -7/16, 13/16, -27/64, 45/64, -3/8, 37/64, -1/4, 7/16, -5/64,`
-    7/64, 1/8, 0, 1/16}}
-Found: {9.791121179618907, {3/2, 39/32, 35/32, 37/32, 77/64, 87/64, 35/32,`
-    1}, {29/32, -7/16, 13/16, -13/32, 47/64, -3/8, 19/32, -9/32, 7/16, -5/64,`
-    7/64, 1/8, 0, 1/16}}
-   */
    /*+1/-1 butterflies (required for FIR, PR, LP).*/
    t[15]=_x[0]-_x[15];
    t[14]=_x[1]-_x[14];
@@ -548,45 +538,45 @@ Found: {9.791121179618907, {3/2, 39/32, 35/32, 37/32, 77/64, 87/64, 35/32,`
    /*U filter (arbitrary invertible, omitted).*/
    /*V filter (arbitrary invertible, can be optimized for speed or accuracy).*/
    /*Scaling factors: the biorthogonal part.*/
-   t[8]=(t[8]<<5)/23;     /*1.41*/
+   t[8]=t[8]*45>>5;       /*1.40568*/
    t[8]+=-t[8]>>15&1;
    /*if(t[8]>0)t[8]++;*/
-   t[9]=(t[9]<<5)/29;     /*1.09*/
+   t[9]=t[9]*37>>5;       /*1.16942*/
    t[9]+=-t[9]>>15&1;
    /*if(t[9]>0)t[9]++;*/
-   t[10]=(t[10]<<4)/15;   /*1.06*/
+   t[10]=t[10]*73>>6;     /*1.12797*/
    t[10]+=-t[10]>>15&1;
    /*if(t[10]>0)t[10]++;*/
-   t[11]=(t[11]<<4)/15;   /*1.06*/
+   t[11]=t[11]*71>>6;     /*1.11222*/
    t[11]+=-t[11]>>15&1;
    /*if(t[11]>0)t[11]++;*/
-   t[12]=(t[12]<<4)/15;   /*1.06*/
+   t[12]=t[12]*67>>6;     /*1.09528*/
    t[12]+=-t[12]>>15&1;
    /*if(t[12]>0)t[12]++;*/
-   t[13]=(t[13]<<4)/15;   /*1.06*/
+   t[13]=t[13]*67>>6;     /*1.09264*/
    t[13]+=-t[13]>>15&1;
    /*if(t[13]>0)t[13]++;*/
-   t[14]=(t[14]<<4)/15;   /*1.08*/
+   t[14]=t[14]*67>>6;     /*1.09645*/
    t[14]+=-t[14]>>15&1;
    /*if(t[14]>0)t[14]++;*/
-   t[15]=(t[15]<<5)/29;   /*1.11*/
+   t[15]=t[15]*9>>3;      /*1.11448*/
    t[15]+=-t[15]>>15&1;
    /*if(t[15]>0)t[15]++;*/
    /*Rotations:*/
-   t[15]-=t[14]*3+16>>5;  /*-0.09*/
-   t[14]+=t[15]*13+32>>6; /*0.21*/
-   t[14]-=t[13]*15+32>>6; /*-0.24*/
-   t[13]+=t[14]*3+4>>3;   /*0.37*/
-   t[13]-=t[12]*23+32>>6; /*-0.36*/
-   t[12]+=t[13]+1>>1;     /*0.50*/
-   t[12]-=t[11]*15+16>>5; /*-0.46*/
-   t[11]+=t[12]*19+16>>5; /*0.60*/
-   t[11]-=t[10]*17+16>>5; /*-0.53*/
-   t[10]+=t[11]*43+32>>6; /*0.67*/
-   t[10]-=t[9]*9+8>>4;    /*-0.56*/
-   t[9]+=t[10]*47+32>>6;  /*0.74*/
-   t[9]-=t[8]*15+16>>5;   /*-0.47*/
-   t[8]+=t[9]*55+32>>6;   /*0.86*/
+   t[9]-=t[8]*3+4>>3;     /*-0.43396*/
+   t[10]-=t[9]*23+32>>6;  /*-0.43365*/
+   t[11]-=t[10]*17+32>>6; /*-0.41797*/
+   t[12]-=t[11]*3+8>>4;   /*-0.36992*/
+   t[13]-=t[12]*7+16>>5;  /*-0.29805*/
+   t[14]-=t[13]*13+32>>6; /*-0.19938*/
+   t[15]-=t[14]*7+32>>6;  /*-0.05825*/
+   t[14]+=t[15]*11+32>>6; /* 0.23047*/
+   t[13]+=t[14]+2>>2;     /* 0.391415*/
+   t[12]+=t[13]*9+16>>5;  /* 0.521553*/
+   t[11]+=t[12]*11+16>>5; /* 0.627901*/
+   t[10]+=t[11]*31+32>>6; /* 0.726061*/
+   t[9]+=t[10]*5+4>>3;    /* 0.818859*/
+   t[8]+=t[9]*25+16>>5;   /* 0.912145*/
    /*More +1/-1 butterflies (required for FIR, PR, LP).*/
    t[0]+=t[15]>>1;
    _y[0]=(od_coeff)t[0];
@@ -632,28 +622,28 @@ void od_post_filter16(od_coeff _x[16],const od_coeff _y[16]){
    t[2]=_y[2]-(t[13]>>1);
    t[1]=_y[1]-(t[14]>>1);
    t[0]=_y[0]-(t[15]>>1);
-   t[8]-=t[9]*55+32>>6;
-   t[9]+=t[8]*15+16>>5;
-   t[9]-=t[10]*47+32>>6;
-   t[10]+=t[9]*9+8>>4;
-   t[10]-=t[11]*43+32>>6;
-   t[11]+=t[10]*17+16>>5;
-   t[11]-=t[12]*19+16>>5;
-   t[12]+=t[11]*15+16>>5;
-   t[12]-=t[13]+1>>1;
-   t[13]+=t[12]*23+32>>6;
-   t[13]-=t[14]*3+4>>3;
-   t[14]+=t[13]*15+32>>6;
-   t[14]-=t[15]*13+32>>6;
-   t[15]+=t[14]*3+16>>5;
-   t[15]=t[15]*29>>5;
-   t[14]=t[14]*15>>4;
-   t[13]=t[13]*15>>4;
-   t[12]=t[12]*15>>4;
-   t[11]=t[11]*15>>4;
-   t[10]=t[10]*15>>4;
-   t[9]=t[9]*29>>5;
-   t[8]=t[8]*23>>5;
+   t[8]-=t[9]*25+16>>5;
+   t[9]-=t[10]*5+4>>3;
+   t[10]-=t[11]*31+32>>6;
+   t[11]-=t[12]*11+16>>5;
+   t[12]-=t[13]*9+16>>5;
+   t[13]-=t[14]+2>>2;
+   t[14]-=t[15]*11+32>>6;
+   t[15]+=t[14]*7+32>>6;
+   t[14]+=t[13]*13+32>>6;
+   t[13]+=t[12]*7+16>>5;
+   t[12]+=t[11]*3+8>>4;
+   t[11]+=t[10]*17+32>>6;
+   t[10]+=t[9]*23+32>>6;
+   t[9]+=t[8]*3+4>>3;
+   t[15]=(t[15]<<3)/9;
+   t[14]=(t[14]<<6)/67;
+   t[13]=(t[13]<<6)/67;
+   t[12]=(t[12]<<6)/67;
+   t[11]=(t[11]<<6)/71;
+   t[10]=(t[10]<<6)/73;
+   t[9]=(t[9]<<5)/37;
+   t[8]=(t[8]<<5)/45;
    t[0]+=t[15]>>1;
    _x[0]=(od_coeff)t[0];
    t[1]+=t[14]>>1;
@@ -695,24 +685,37 @@ int main(void){
    int         i;
    int         j;
    for(j=0;j<16;j++)min[j]=max[j]=mini[j]=maxi[j]=0;
-   for(i=0;i<256;i++){
+   for(i=0;i<65536;i++){
       od_coeff x[16];
-      for(j=0;j<16;j++)x[j]=i>>j&1?1325:-1325;
-      od_pre_filter8(x,x);
-      /*for(j=0;j<16;j++){
-         if(x[j]<min[j]){
-            min[j]=x[j];
+      od_coeff y[16];
+      od_coeff x2[16];
+      for(j=0;j<16;j++)x[j]=i>>j&1?676:-676;
+      od_pre_filter16(y,x);
+      for(j=0;j<16;j++){
+         if(y[j]<min[j]){
+            min[j]=y[j];
             mini[j]=i;
          }
-         else if(x[j]>max[j]){
-            max[j]=x[j];
+         else if(y[j]>max[j]){
+            max[j]=y[j];
             maxi[j]=i;
          }
-      }*/
+      }
+      od_post_filter16(x2,y);
+      for(j=0;j<16;j++)if(x[j]!=x2[j]){
+         printf("Mismatch:\n");
+         printf("in:    ");
+         for(j=0;j<16;j++)printf(" %i",x[j]);
+         printf("\nxform: ");
+         for(j=0;j<16;j++)printf(" %i",y[j]);
+         printf("\nout:    ");
+         for(j=0;j<16;j++)printf(" %i",x2[j]);
+         printf("\n\n");
+      }
    }
-   for(j=0;j<18;j++){
-      printf("Min: %5i  ",minv[j]);
-      printf("Max: %5i\n",maxv[j]);
+   for(j=0;j<16;j++){
+      printf("Min: %5i  ",min[j]);
+      printf("Max: %5i\n",max[j]);
    }
    return 0;
 }