iSAC: Move global trig tables into the codec instance

webrtc/modules/audio_coding/codecs/isac/main/source/transform.c

 /*
  *  Copyright (c) 2011 The WebRTC project authors. All Rights Reserved.
  *
  *  Use of this source code is governed by a BSD-style license
  *  that can be found in the LICENSE file in the root of the source
  *  tree. An additional intellectual property rights grant can be found
  *  in the file PATENTS.  All contributing project authors may
  *  be found in the AUTHORS file in the root of the source tree.
  */
 
 #include "settings.h"
 #include "fft.h"
 #include "codec.h"
 #include "os_specific_inline.h"
 #include <math.h>
 
-static double costab1[FRAMESAMPLES_HALF];
-static double sintab1[FRAMESAMPLES_HALF];
-static double costab2[FRAMESAMPLES_QUARTER];
-static double sintab2[FRAMESAMPLES_QUARTER];
-
-void WebRtcIsac_InitTransform()
-{
+void WebRtcIsac_InitTransform(TransformTables* tables) {
   int k;
   double fact, phase;
 
   fact = PI / (FRAMESAMPLES_HALF);
   phase = 0.0;
   for (k = 0; k < FRAMESAMPLES_HALF; k++) {
-    costab1[k] = cos(phase);
-    sintab1[k] = sin(phase);
+    tables->costab1[k] = cos(phase);
+    tables->sintab1[k] = sin(phase);
     phase += fact;
   }
 
   fact = PI * ((double) (FRAMESAMPLES_HALF - 1)) / ((double) FRAMESAMPLES_HALF);
   phase = 0.5 * fact;
   for (k = 0; k < FRAMESAMPLES_QUARTER; k++) {
-    costab2[k] = cos(phase);
-    sintab2[k] = sin(phase);
+    tables->costab2[k] = cos(phase);
+    tables->sintab2[k] = sin(phase);
     phase += fact;
   }
 }
 
-
-void WebRtcIsac_Time2Spec(double *inre1,
-                         double *inre2,
-                         int16_t *outreQ7,
-                         int16_t *outimQ7,
-                         FFTstr *fftstr_obj)
-{
-
+void WebRtcIsac_Time2Spec(const TransformTables* tables,
+                          double* inre1,
+                          double* inre2,
+                          int16_t* outreQ7,
+                          int16_t* outimQ7,
+                          FFTstr* fftstr_obj) {
   int k;
   int dims[1];
   double tmp1r, tmp1i, xr, xi, yr, yi, fact;
   double tmpre[FRAMESAMPLES_HALF], tmpim[FRAMESAMPLES_HALF];
 
 
   dims[0] = FRAMESAMPLES_HALF;
 
 
   /* Multiply with complex exponentials and combine into one complex vector */
   fact = 0.5 / sqrt(FRAMESAMPLES_HALF);
   for (k = 0; k < FRAMESAMPLES_HALF; k++) {
-    tmp1r = costab1[k];
-    tmp1i = sintab1[k];
+    tmp1r = tables->costab1[k];

[the sun, 2015/06/15 15:08:51]

Premature pessimization?

I hope the dereference is outside the loop. What code is generated?

[kwiberg-webrtc, 2015/06/15 21:57:28]

I'm glad you asked. Given this small test program ("double-1") which
reproduces the essential feature, indexing into an array of ~700
doubles that isn't the first member in the struct we have a pointer
to:

  struct X {
    double wooz[700];
    double frob[700];
  };

  double ss(struct X* x) {
    double p = 0;
    int i;
    for (i = 0; i < 700; ++i)
      p += x->frob[i];
    return p;
  }

gcc-4.8 -O2 generates the following ("gcc-double-1"):

  0000000000000000 <ss>:
     0:	66 0f 57 c0          	xorpd  %xmm0,%xmm0
     4:	31 c0                	xor    %eax,%eax
     6:	66 2e 0f 1f 84 00 00 	nopw   %cs:0x0(%rax,%rax,1)
     d:	00 00 00 
    10:	f2 0f 58 84 07 e0 15 	addsd  0x15e0(%rdi,%rax,1),%xmm0
    17:	00 00 
    19:	48 83 c0 08          	add    $0x8,%rax
    1d:	48 3d e0 15 00 00    	cmp    $0x15e0,%rax
    23:	75 eb                	jne    10 <ss+0x10>
    25:	f3 c3                	repz retq 

The actual loop is the four instructions at 0x10-0x23. %rdi is x, %rax
is 8*i, %xmm0 is p, and 0x15e0 is 700*8; the entire loop body (p +=
x->frob[i]) is done by the single instruction at 0x10. With the
following variation, where the C code manually tries to lift the
addition out of the loop ("double-2"):

  double* d = x->frob;
  for (i = 0; i < 700; ++i)
    p += d[i];

we get pretty much exactly the same inner-loop code ("gcc-double-2"),
except that the instruction that dereferences and accumulates uses a
slightly simpler addresing mode:

    10:	f2 0f 58 04 02       	addsd  (%rdx,%rax,1),%xmm0

(Here %rdx is d, and %rax and %eax are 8*i and p as before.)

The exact same programs with arrays of ints instead of doubles gave
pretty much identically structured code in both cases (gcc-int-1
structured as gcc-double-1, gcc-int-2 structured as gcc-double-2).

So what about clang? I'm glad you asked. clang-3.4 -O2 doesn't care
that we tried to micro-optimize the code, and generates identical code
for the two programs both in the int and double case (i.e.,
clang-double-1 == clang-double-2 and clang-int-1 ==
clang-int-2). clang-double-1/2 looks pretty similar to gcc-double-1;
in particular, the addition of the constant offset has not been lifted
out of the loop:

  0000000000000000 <ss>:
     0:	31 c0                	xor    %eax,%eax
     2:	0f 57 c0             	xorps  %xmm0,%xmm0
     5:	66 66 2e 0f 1f 84 00 	data32 nopw %cs:0x0(%rax,%rax,1)
     c:	00 00 00 00 
    10:	f2 0f 58 84 c7 e0 15 	addsd  0x15e0(%rdi,%rax,8),%xmm0
    17:	00 00 
    19:	48 ff c0             	inc    %rax
    1c:	48 3d bc 02 00 00    	cmp    $0x2bc,%rax
    22:	75 ec                	jne    10 <ss+0x10>
    24:	c3                   	retq   

And clang-int-1/2 is a thing of terrifying beauty:

  0000000000000000 <ss>:
     0:	31 c0                	xor    %eax,%eax
     2:	66 0f ef c0          	pxor   %xmm0,%xmm0
     6:	66 0f ef c9          	pxor   %xmm1,%xmm1
     a:	66 0f 1f 44 00 00    	nopw   0x0(%rax,%rax,1)
    10:	66 0f 6f d1          	movdqa %xmm1,%xmm2
    14:	66 0f 6f d8          	movdqa %xmm0,%xmm3
    18:	f3 0f 6f 84 87 f0 0a 	movdqu 0xaf0(%rdi,%rax,4),%xmm0
    1f:	00 00 
    21:	f3 0f 6f 8c 87 00 0b 	movdqu 0xb00(%rdi,%rax,4),%xmm1
    28:	00 00 
    2a:	66 0f fe c3          	paddd  %xmm3,%xmm0
    2e:	66 0f fe ca          	paddd  %xmm2,%xmm1
    32:	48 83 c0 08          	add    $0x8,%rax
    36:	48 3d b8 02 00 00    	cmp    $0x2b8,%rax
    3c:	75 d2                	jne    10 <ss+0x10>
    3e:	66 0f fe c8          	paddd  %xmm0,%xmm1
    42:	66 0f 6f c1          	movdqa %xmm1,%xmm0
    46:	0f 12 c0             	movhlps %xmm0,%xmm0
    49:	66 0f fe c1          	paddd  %xmm1,%xmm0
    4d:	66 0f 70 c8 01       	pshufd $0x1,%xmm0,%xmm1
    52:	66 0f fe c8          	paddd  %xmm0,%xmm1
    56:	66 0f 7e c8          	movd   %xmm1,%eax
    5a:	31 c9                	xor    %ecx,%ecx
    5c:	0f 1f 40 00          	nopl   0x0(%rax)
    60:	03 84 8f d0 15 00 00 	add    0x15d0(%rdi,%rcx,4),%eax
    67:	48 ff c1             	inc    %rcx
    6a:	48 83 f9 04          	cmp    $0x4,%rcx
    6e:	75 f0                	jne    60 <ss+0x60>
    70:	c3                   	retq   

The inner loop is the instructions at 0x10-0x3c; it reads two vectors
of four ints at a time and sums them in two separate four-int
accumulators that it adds up after the loop is done. Because 700 = 87
* 8 + 4, after 87 iteration of the main loop there are four more
values to add up in a separate loop at the end.

So how fast are all these variants? Calling each function 1e7 times
and taking the fastest of 3 runs, we see that the double programs are
all pretty much identical in performance:

  6.43s clang-double-1/2
  6.44s gcc-double-1
  6.45s gcc-double-2

And with ints, clang's bad-ass code runs circles around gcc:

  0.85s clang-int-1/2
  3.45s gcc-int-1
  3.49s gcc-int-2

We can also see that the micro-optimization (the difference between
programs 1 and 2) in both the double and int cases hurts performance
marginally on gcc; clang, as we saw above, ignored it entirely and
produced identical code for the two programs.

[the sun, 2015/06/16 08:19:54]

Hey, now I AM glad I asked. That was super informative! Thanks so much for the
lesson!

So I guess with f32 we should also see auto-vectorization kick in with clang.
Why didn't it do it in the double case? Not enough gain? Target only SSE1?

[kwiberg-webrtc, 2015/06/16 08:31:44]

I know just enough assembly to decrypt what code like the above is doing. In
order to be able to write it, I'd have to spend a few weeks with instruction set
reference manuals and running benchmarks. So I don't really have an idea if
vectorization is inherently harder with doubles than with ints, or if clang's
optimizer just isn't general enough. Or whether float code is more like int or
double code in this regard.

[the sun, 2015/06/16 08:46:36]

I was just thinking that it did auto-vectorization in the int case because it
could fit 4 ints in a register, but maybe the overhead wasn't deemed a good
tradeoff when only 2 doubles per register was possible.

Since it would be able to fit 4 f32:s too, hopefully it would do auto-vec there
as well. I don't see that the compiler would have a harder time doing it with
floats, especially since it is already using SSE for fp ops.

+    tmp1i = tables->sintab1[k];
     tmpre[k] = (inre1[k] * tmp1r + inre2[k] * tmp1i) * fact;
     tmpim[k] = (inre2[k] * tmp1r - inre1[k] * tmp1i) * fact;
   }
 
 
   /* Get DFT */
   WebRtcIsac_Fftns(1, dims, tmpre, tmpim, -1, 1.0, fftstr_obj);
 
   /* Use symmetry to separate into two complex vectors and center frames in time around zero */
   for (k = 0; k < FRAMESAMPLES_QUARTER; k++) {
     xr = tmpre[k] + tmpre[FRAMESAMPLES_HALF - 1 - k];
     yi = -tmpre[k] + tmpre[FRAMESAMPLES_HALF - 1 - k];
     xi = tmpim[k] - tmpim[FRAMESAMPLES_HALF - 1 - k];
     yr = tmpim[k] + tmpim[FRAMESAMPLES_HALF - 1 - k];
 
-    tmp1r = costab2[k];
-    tmp1i = sintab2[k];
+    tmp1r = tables->costab2[k];
+    tmp1i = tables->sintab2[k];
     outreQ7[k] = (int16_t)WebRtcIsac_lrint((xr * tmp1r - xi * tmp1i) * 128.0);
     outimQ7[k] = (int16_t)WebRtcIsac_lrint((xr * tmp1i + xi * tmp1r) * 128.0);
     outreQ7[FRAMESAMPLES_HALF - 1 - k] = (int16_t)WebRtcIsac_lrint((-yr * tmp1i - yi * tmp1r) * 128.0);
     outimQ7[FRAMESAMPLES_HALF - 1 - k] = (int16_t)WebRtcIsac_lrint((-yr * tmp1r + yi * tmp1i) * 128.0);
   }
 }
 
-
-void WebRtcIsac_Spec2time(double *inre, double *inim, double *outre1, double *outre2, FFTstr *fftstr_obj)
-{
-
+void WebRtcIsac_Spec2time(const TransformTables* tables,
+                          double* inre,
+                          double* inim,
+                          double* outre1,
+                          double* outre2,
+                          FFTstr* fftstr_obj) {
   int k;
   double tmp1r, tmp1i, xr, xi, yr, yi, fact;
 
   int dims;
 
   dims = FRAMESAMPLES_HALF;
 
   for (k = 0; k < FRAMESAMPLES_QUARTER; k++) {
     /* Move zero in time to beginning of frames */
-    tmp1r = costab2[k];
-    tmp1i = sintab2[k];
+    tmp1r = tables->costab2[k];
+    tmp1i = tables->sintab2[k];
     xr = inre[k] * tmp1r + inim[k] * tmp1i;
     xi = inim[k] * tmp1r - inre[k] * tmp1i;
     yr = -inim[FRAMESAMPLES_HALF - 1 - k] * tmp1r - inre[FRAMESAMPLES_HALF - 1 - k] * tmp1i;
     yi = -inre[FRAMESAMPLES_HALF - 1 - k] * tmp1r + inim[FRAMESAMPLES_HALF - 1 - k] * tmp1i;
 
     /* Combine into one vector,  z = x + j * y */
     outre1[k] = xr - yi;
     outre1[FRAMESAMPLES_HALF - 1 - k] = xr + yi;
     outre2[k] = xi + yr;
     outre2[FRAMESAMPLES_HALF - 1 - k] = -xi + yr;
   }
 
 
   /* Get IDFT */
   WebRtcIsac_Fftns(1, &dims, outre1, outre2, 1, FRAMESAMPLES_HALF, fftstr_obj);
 
 
   /* Demodulate and separate */
   fact = sqrt(FRAMESAMPLES_HALF);
   for (k = 0; k < FRAMESAMPLES_HALF; k++) {
-    tmp1r = costab1[k];
-    tmp1i = sintab1[k];
+    tmp1r = tables->costab1[k];
+    tmp1i = tables->sintab1[k];
     xr = (outre1[k] * tmp1r - outre2[k] * tmp1i) * fact;
     outre2[k] = (outre2[k] * tmp1r + outre1[k] * tmp1i) * fact;
     outre1[k] = xr;
   }
 }