Doing the spreading with a "pseudo-fractional-Hadamard" transform
authorJean-Marc Valin <jean-marc.valin@usherbrooke.ca>
Sat, 3 Apr 2010 13:23:29 +0000 (09:23 -0400)
committerJean-Marc Valin <jean-marc.valin@usherbrooke.ca>
Sat, 3 Apr 2010 13:23:29 +0000 (09:23 -0400)
libcelt/vq.c

index 7323561..56096a1 100644 (file)
 #define M_PI 3.141592653
 #endif
 
+static void frac_hadamard1(celt_norm *X, int len, int stride, celt_word16 c, celt_word16 s)
+{
+   int j;
+   celt_norm *x, *y;
+   celt_norm * end;
+
+   j = 0;
+   x = X;
+   y = X+stride;
+   end = X+len;
+   do
+   {
+      celt_norm x1, x2;
+      x1 = *x;
+      x2 = *y;
+      *x++ = EXTRACT16(SHR32(MULT16_16(c,x1) + MULT16_16(s,x2),15));
+      *y++ = EXTRACT16(SHR32(MULT16_16(s,x1) - MULT16_16(c,x2),15));
+      j++;
+      if (j>=stride)
+      {
+         j=0;
+         x+=stride;
+         y+=stride;
+      }
+   } while (y<end);
+
+   /* Reverse samples so that the next level starts from the other end */
+   for (j=0;j<len>>1;j++)
+   {
+      celt_norm tmp = X[j];
+      X[j] = X[len-j-1];
+      X[len-j-1] = tmp;
+   }
+}
+
+#define MAX_LEVELS 8
 static void exp_rotation(celt_norm *X, int len, int dir, int stride, int K)
 {
-   int i, k, iter;
-   celt_word16 c, s;
+   int i, N=0;
    celt_word16 gain, theta;
-   celt_norm *Xptr;
-   gain = celt_div((celt_word32)MULT16_16(Q15_ONE,len),(celt_word32)(3+len+6*K));
+   int istride[MAX_LEVELS];
+   celt_word16 c, s;
+
+   do {
+      istride[N] = stride;
+      stride *= 2;
+      N++;
+   } while (N<MAX_LEVELS && stride < len);
+
+   gain = celt_div((celt_word32)MULT16_16(Q15_ONE,len),(celt_word32)(3+len+4*K));
    /* FIXME: Make that HALF16 instead of HALF32 */
-   theta = SUB16(Q15ONE, HALF32(MULT16_16_Q15(gain,gain)));
-   /*if (len==30)
-   {
-   for (i=0;i<len;i++)
-   X[i] = 0;
-   X[14] = 1;
-}*/ 
+   theta = HALF32(MULT16_16_Q15(gain,gain));
    c = celt_cos_norm(EXTEND32(theta));
-   s = dir*celt_cos_norm(EXTEND32(SUB16(Q15ONE,theta))); /*  sin(theta) */
-   if (len > 8*stride)
-      stride *= len/(8*stride);
-   iter = 1;
-   for (k=0;k<iter;k++)
+   s = celt_cos_norm(EXTEND32(SUB16(Q15ONE,theta))); /*  sin(theta) */
+
+   if (dir > 0)
    {
-      /* We could use MULT16_16_P15 instead of MULT16_16_Q15 for more accuracy, 
-      but at this point, I really don't think it's necessary */
-      Xptr = X;
-      for (i=0;i<len-stride;i++)
-      {
-         celt_norm x1, x2;
-         x1 = Xptr[0];
-         x2 = Xptr[stride];
-         Xptr[stride] = EXTRACT16(SHR32(MULT16_16(c,x2) + MULT16_16(s,x1),15));
-         *Xptr++      = EXTRACT16(SHR32(MULT16_16(c,x1) - MULT16_16(s,x2),15));
-      }
-      Xptr = &X[len-2*stride-1];
-      for (i=len-2*stride-1;i>=0;i--)
-      {
-         celt_norm x1, x2;
-         x1 = Xptr[0];
-         x2 = Xptr[stride];
-         Xptr[stride] = EXTRACT16(SHR32(MULT16_16(c,x2) + MULT16_16(s,x1),15));
-         *Xptr--      = EXTRACT16(SHR32(MULT16_16(c,x1) - MULT16_16(s,x2),15));
-      }
+      for (i=0;i<N;i++)
+         frac_hadamard1(X, len, istride[i], c, s);
+   } else {
+      for (i=N-1;i>=0;i--)
+         frac_hadamard1(X, len, istride[i], c, s);
    }
-   /*if (len==30)
+
+   /* Undo last reversal */
+   for (i=0;i<len>>1;i++)
    {
-   for (i=0;i<len;i++)
-   printf ("%f ", X[i]);
-   printf ("\n");
-   exit(0);
-}*/
+      celt_norm tmp = X[i];
+      X[i] = X[len-i-1];
+      X[len-i-1] = tmp;
+   }
 }
 
-
 /** Takes the pitch vector and the decoded residual vector, computes the gain
     that will give ||p+g*y||=1 and mixes the residual with the pitch. */
 static void normalise_residual(int * restrict iy, celt_norm * restrict X, int N, int K, celt_word32 Ryy)