Makes surround_analysis() work in fixed-point

src/opus_multistream_encoder.c

@@ -179,6 +179,45 @@ static void channel_pos(int channels, int pos[8])
    }
 }
 
+#if 1
+/* Computes a rough approximation of log2(2^a + 2^b) */
+static opus_val16 logSum(opus_val16 a, opus_val16 b)
+{
+   opus_val16 max;
+   opus_val32 diff;
+   opus_val16 frac;
+   static const opus_val16 diff_table[17] = {
+         QCONST16(0.5000000f, DB_SHIFT), QCONST16(0.2924813f, DB_SHIFT), QCONST16(0.1609640f, DB_SHIFT), QCONST16(0.0849625f, DB_SHIFT),
+         QCONST16(0.0437314f, DB_SHIFT), QCONST16(0.0221971f, DB_SHIFT), QCONST16(0.0111839f, DB_SHIFT), QCONST16(0.0056136f, DB_SHIFT),
+         QCONST16(0.0028123f, DB_SHIFT)
+   };
+   int low;
+   if (a>b)
+   {
+      max = a;
+      diff = SUB32(EXTEND32(a),EXTEND32(b));
+   } else {
+      max = b;
+      diff = SUB32(EXTEND32(b),EXTEND32(a));
+   }
+   if (diff >= QCONST16(8.f, DB_SHIFT))
+      return max;
+#ifdef FIXED_POINT
+   low = SHR32(diff, DB_SHIFT-1);
+   frac = SHL16(diff - SHL16(low, DB_SHIFT-1), 16-DB_SHIFT);
+#else
+   low = floor(2*diff);
+   frac = 2*diff - low;
+#endif
+   return max + diff_table[low] + MULT16_16_Q15(frac, SUB16(diff_table[low+1], diff_table[low]));
+}
+#else
+opus_val16 logSum(opus_val16 a, opus_val16 b)
+{
+   return log2(pow(4, a)+ pow(4, b))/2;
+}
+#endif
+
 void surround_analysis(const CELTMode *celt_mode, const void *pcm, opus_val16 *bandLogE, opus_val32 *mem, opus_val32 *preemph_mem,
       int len, int overlap, int channels, int rate, opus_copy_channel_in_func copy_channel_in
 )
@@ -190,7 +229,6 @@ void surround_analysis(const CELTMode *celt_mode, const void *pcm, opus_val16 *b
    int pos[8] = {0};
    int upsample;
    opus_val32 bandE[21];
-   opus_val32 maskE[3][21];
    opus_val16 maskLogE[3][21];
    VARDECL(opus_val32, in);
    VARDECL(opus_val16, x);
@@ -202,9 +240,9 @@ void surround_analysis(const CELTMode *celt_mode, const void *pcm, opus_val16 *b
 
    channel_pos(channels, pos);
 
-   for (c=0;c<2;c++)
+   for (c=0;c<3;c++)
       for (i=0;i<21;i++)
-         maskE[c][i] = 0;
+         maskLogE[c][i] = -QCONST16(28.f, DB_SHIFT);
 
    upsample = resampling_factor(rate);
    for (c=0;c<channels;c++)
@@ -224,24 +262,23 @@ void surround_analysis(const CELTMode *celt_mode, const void *pcm, opus_val16 *b
       }
 
       compute_band_energies(celt_mode, freq, bandE, 21, 1, 1<<LM);
-      /* FIXME: Figure out how to square bandE[] in fixed-point */
+      amp2Log2(celt_mode, 21, 21, bandE, bandLogE+21*c, 1);
       if (pos[c]==1)
       {
          for (i=0;i<21;i++)
-            maskE[0][i] += bandE[i]*bandE[i];
+            maskLogE[0][i] = logSum(maskLogE[0][i], bandLogE[21*c+i]);
       } else if (pos[c]==3)
       {
          for (i=0;i<21;i++)
-            maskE[1][i] += bandE[i]*bandE[i];
+            maskLogE[2][i] = logSum(maskLogE[2][i], bandLogE[21*c+i]);
       } else if (pos[c]==2)
       {
          for (i=0;i<21;i++)
          {
-            maskE[0][i] += HALF32(bandE[i]*bandE[i]);
-            maskE[1][i] += HALF32(bandE[i]*bandE[i]);
+            maskLogE[0][i] = logSum(maskLogE[0][i], bandLogE[21*c+i]-QCONST16(.5f, DB_SHIFT));
+            maskLogE[2][i] = logSum(maskLogE[2][i], bandLogE[21*c+i]-QCONST16(.5f, DB_SHIFT));
          }
       }
-      amp2Log2(celt_mode, 21, 21, bandE, bandLogE+21*c, 1);
 #if 0
       for (i=0;i<21;i++)
          printf("%f ", bandLogE[21*c+i]);
@@ -254,16 +291,10 @@ void surround_analysis(const CELTMode *celt_mode, const void *pcm, opus_val16 *b
       OPUS_COPY(mem+c*overlap, in+len, overlap);
    }
    for (i=0;i<21;i++)
-      maskE[2][i] = MIN32(maskE[0][i],maskE[1][i]);
+      maskLogE[1][i] = MIN32(maskLogE[0][i],maskLogE[2][i]);
    for (c=0;c<3;c++)
       for (i=0;i<21;i++)
-         maskE[c][i] = sqrt(maskE[c][i]*2/(channels-1));
-   /* Left mask */
-   amp2Log2(celt_mode, 21, 21, &maskE[0][0], &maskLogE[0][0], 1);
-   /* Right mask */
-   amp2Log2(celt_mode, 21, 21, &maskE[1][0], &maskLogE[2][0], 1);
-   /* Centre mask */
-   amp2Log2(celt_mode, 21, 21, &maskE[2][0], &maskLogE[1][0], 1);
+         maskLogE[c][i] += QCONST16(.5f, DB_SHIFT)*log2(2.f/(channels-1));
 #if 0
    for (c=0;c<3;c++)
    {