<div dir="ltr">Hi Jean-Marc,<div><br></div><div>Thanks!</div><div><br></div><div>Please find the 2 updated patches which only optimize stride 2 case and keep the bit exactness. They have passed our internal tests as usual.</div><div><br></div><div>Thanks,</div><div>Linfeng</div></div><div class="gmail_extra"><br><div class="gmail_quote">On Mon, May 15, 2017 at 9:36 AM, Jean-Marc Valin <span dir="ltr"><<a href="mailto:jmvalin@jmvalin.ca" target="_blank">jmvalin@jmvalin.ca</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi Linfeng,<br>
<br>
Sorry for the delay -- I was actually trying to think of the best option<br>
here. For now, my preference would be to keep things bit-exact, but<br>
should there be more similar optimizations relying on 64-bit<br>
multiplication results, then we could consider having a special option<br>
to enable those (even in C).<br>
<br>
Cheers,<br>
<br>
  Â  Â  Â  Jean-Marc<br>
<span class="im HOEnZb"><br>
On 08/05/17 12:12 PM, Linfeng Zhang wrote:<br>
> Ping for comments.<br>
><br>
> Thanks,<br>
> Linfeng<br>
><br>
> On Wed, Apr 26, 2017 at 2:15 PM, Linfeng Zhang <<a href="mailto:linfengz@google.com">linfengz@google.com</a><br>
</span><span class="im HOEnZb">> <mailto:<a href="mailto:linfengz@google.com">linfengz@google.com</a>>> wrote:<br>
><br>
>  Â  Â On Tue, Apr 25, 2017 at 10:31 PM, Jean-Marc Valin<br>
</span><div class="HOEnZb"><div class="h5">>  Â  Â <<a href="mailto:jmvalin@jmvalin.ca">jmvalin@jmvalin.ca</a> <mailto:<a href="mailto:jmvalin@jmvalin.ca">jmvalin@jmvalin.ca</a>>> wrote:<br>
><br>
><br>
>  Â  Â  Â  Â > A_Q28 is split to 2 14-bit (or 16-bit, whatever) integers, to make the<br>
>  Â  Â  Â  Â > multiplication operation within 32-bits. NEON can do 32-bit x 32-bit =<br>
>  Â  Â  Â  Â > 64-bit using 'int64x2_t vmull_s32(int32x2_t a, int32x2_t b)', and it<br>
>  Â  Â  Â  Â > could possibly be faster and less rounding/shifting errors than above C<br>
>  Â  Â  Â  Â > code. But it may increase difficulties for other CPUs not supporting<br>
>  Â  Â  Â  Â > 32-bit multiplication.<br>
><br>
>  Â  Â  Â  Â OK, so I'm not totally opposed to that, but it increases the<br>
>  Â  Â  Â  Â testing/maintenance cost so it needs to be worth it. So the<br>
>  Â  Â  Â  Â question is<br>
>  Â  Â  Â  Â how much speedup can you get and how close you can make the<br>
>  Â  Â  Â  Â result to<br>
>  Â  Â  Â  Â the original function. If you can make the output be always<br>
>  Â  Â  Â  Â within one<br>
>  Â  Â  Â  Â of two LSBs of the C version, then the asm check can simply be a<br>
>  Â  Â  Â  Â little<br>
>  Â  Â  Â  Â bit more lax than usual. Otherwise it becomes more complicated. This<br>
>  Â  Â  Â  Â isn't a function that scares me too much about going<br>
>  Â  Â  Â  Â non-bitexact, but<br>
>  Â  Â  Â  Â it's also not one of the big complexity costs either. In any<br>
>  Â  Â  Â  Â case, let<br>
>  Â  Â  Â  Â me know what you find.<br>
><br>
><br>
>  Â  Â Tested the proposed NEON optimization change, it can increase the<br>
>  Â  Â whole encoder speed by about 1% at Complexity 8 and 10 for stride =<br>
>  Â  Â 2 (compared to the NEON optimization in the previous patch), and<br>
>  Â  Â 0.7% at Complexity 8 and 1.1% at Complexity 10 for stride = 1<br>
>  Â  Â (compared to the original C code).<br>
><br>
>  Â  Â Unfortunately, the truncating difference accumulates in S[] and its<br>
>  Â  Â difference cannot be easily bounded. The difference in out[] may<br>
>  Â  Â somehow be bounded to 5 in my quick testing, but is not guaranteed<br>
>  Â  Â to other inputs. So maybe comparing bit exactness with the<br>
>  Â  Â following silk_biquad_alt_c_<wbr>MulSingleAQ28() is better.<br>
><br>
>  Â  Â Please let me know the decision (whether keeping the original NEON<br>
>  Â  Â (stride 2 only) or choosing the new NEON (both stride 1 and 2) which<br>
>  Â  Â optimizes following silk_biquad_alt_c_<wbr>MulSingleAQ28()), and I'll<br>
>  Â  Â wrap up the patch.<br>
><br>
>  Â  Â Here attached the corresponding C<br>
>  Â  Â code silk_biquad_alt_c_<wbr>MulSingleAQ28() for your information. It uses<br>
>  Â  Â 64-bit multiplications and is about 0.6% - 0.9% slower (for the<br>
>  Â  Â whole encoder) than the original C code using 32-bit multiplications<br>
>  Â  Â (for both strides and the same stride 2 unrolling).<br>
><br>
>  Â  Â void silk_biquad_alt_c_<wbr>MulSingleAQ28(<br>
>  Â  Â  Â  Â const opus_int16  Â  Â  Â  Â  Â  *in,  Â  Â  Â  Â  Â  Â  Â  /* I  Â  Â input<br>
>  Â  Â signal  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â */<br>
>  Â  Â  Â  Â const opus_int32  Â  Â  Â  Â  Â  *B_Q28,  Â  Â  Â  Â  Â  Â /* I  Â  Â MA<br>
>  Â  Â coefficients [3]  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  */<br>
>  Â  Â  Â  Â const opus_int32  Â  Â  Â  Â  Â  *A_Q28,  Â  Â  Â  Â  Â  Â /* I  Â  Â AR<br>
>  Â  Â coefficients [2]  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  */<br>
>  Â  Â  Â  Â opus_int32  Â  Â  Â  Â  Â  Â  Â  Â  *S,  Â  Â  Â  Â  Â  Â  Â  Â /* I/O  Â State<br>
>  Â  Â vector [2*stride]  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  */<br>
>  Â  Â  Â  Â opus_int16  Â  Â  Â  Â  Â  Â  Â  Â  *out,  Â  Â  Â  Â  Â  Â  Â /* O  Â  Â output<br>
>  Â  Â signal  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  */<br>
>  Â  Â  Â  Â const opus_int32  Â  Â  Â  Â  Â  len,  Â  Â  Â  Â  Â  Â  Â  /* I  Â  Â signal<br>
>  Â  Â length (must be even)  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â */<br>
>  Â  Â  Â  Â opus_int  Â  Â  Â  Â  Â  Â  Â  Â  Â  stride  Â  Â  Â  Â  Â  Â  /* I  Â  Â Operate<br>
>  Â  Â on interleaved signal if > 1  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â */<br>
>  Â  Â )<br>
>  Â  Â {<br>
>  Â  Â  Â  Â /* DIRECT FORM II TRANSPOSED (uses 2 element state vector) */<br>
>  Â  Â  Â  Â opus_int  Â k;<br>
><br>
>  Â  Â  Â  Â silk_assert( ( stride == 1 ) || ( stride == 2 ) );<br>
><br>
>  Â  Â  Â  Â if( stride == 1) {<br>
>  Â  Â  Â  Â  Â  Â opus_int32 out32_Q14;<br>
>  Â  Â  Â  Â  Â  Â for( k = 0; k < len; k++ ) {<br>
>  Â  Â  Â  Â  Â  Â  Â  Â /* S[ 0 ], S[ 1 ]: Q12 */<br>
>  Â  Â  Â  Â  Â  Â  Â  Â out32_Q14 = silk_LSHIFT( silk_SMLAWB( S[ 0 ], B_Q28[ 0<br>
>  Â  Â ], in[ k ] ), 2 );<br>
><br>
>  Â  Â  Â  Â  Â  Â  Â  Â S[ 0 ] = S[ 1 ] + silk_RSHIFT_ROUND(<br>
>  Â  Â (opus_int64)out32_Q14 * (-A_Q28[ 0 ]), 30 );<br>
>  Â  Â  Â  Â  Â  Â  Â  Â S[ 0 ] = silk_SMLAWB( S[ 0 ], B_Q28[ 1 ], in[ k ] );<br>
><br>
>  Â  Â  Â  Â  Â  Â  Â  Â S[ 1 ] = silk_RSHIFT_ROUND( (opus_int64)out32_Q14 *<br>
>  Â  Â (-A_Q28[ 1 ]) , 30 );<br>
>  Â  Â  Â  Â  Â  Â  Â  Â S[ 1 ] = silk_SMLAWB( S[ 1 ], B_Q28[ 2 ], in[ k ] );<br>
><br>
>  Â  Â  Â  Â  Â  Â  Â  Â /* Scale back to Q0 and saturate */<br>
>  Â  Â  Â  Â  Â  Â  Â  Â out[ k ] = (opus_int16)silk_SAT16( silk_RSHIFT(<br>
>  Â  Â out32_Q14 + (1<<14) - 1, 14 ) );<br>
>  Â  Â  Â  Â  Â  Â }<br>
>  Â  Â  Â  Â } else {<br>
>  Â  Â  Â  Â  Â  Â opus_int32 out32_Q14[ 2 ];<br>
>  Â  Â  Â  Â  Â  Â for( k = 0; k < len; k++ ) {<br>
>  Â  Â  Â  Â  Â  Â  Â  Â /* S[ 0 ], S[ 1 ]: Q12 */<br>
>  Â  Â  Â  Â  Â  Â  Â  Â out32_Q14[ 0 ] = silk_LSHIFT( silk_SMLAWB( S[ 0 ],<br>
>  Â  Â B_Q28[ 0 ], in[ k * 2 + 0 ] ), 2 );<br>
>  Â  Â  Â  Â  Â  Â  Â  Â out32_Q14[ 1 ] = silk_LSHIFT( silk_SMLAWB( S[ 2 ],<br>
>  Â  Â B_Q28[ 0 ], in[ k * 2 + 1 ] ), 2 );<br>
><br>
>  Â  Â  Â  Â  Â  Â  Â  Â S[ 0 ] = S[ 1 ] + silk_RSHIFT_ROUND(<br>
>  Â  Â (opus_int64)out32_Q14[ 0 ] * (-A_Q28[ 0 ]), 30 );<br>
>  Â  Â  Â  Â  Â  Â  Â  Â S[ 2 ] = S[ 3 ] + silk_RSHIFT_ROUND(<br>
>  Â  Â (opus_int64)out32_Q14[ 1 ] * (-A_Q28[ 0 ]), 30 );<br>
>  Â  Â  Â  Â  Â  Â  Â  Â S[ 0 ] = silk_SMLAWB( S[ 0 ], B_Q28[ 1 ], in[ k * 2 + 0 ] );<br>
>  Â  Â  Â  Â  Â  Â  Â  Â S[ 2 ] = silk_SMLAWB( S[ 2 ], B_Q28[ 1 ], in[ k * 2 + 1 ] );<br>
><br>
>  Â  Â  Â  Â  Â  Â  Â  Â S[ 1 ] = silk_RSHIFT_ROUND( (opus_int64)out32_Q14[ 0 ] *<br>
>  Â  Â (-A_Q28[ 1 ]), 30 );<br>
>  Â  Â  Â  Â  Â  Â  Â  Â S[ 3 ] = silk_RSHIFT_ROUND( (opus_int64)out32_Q14[ 1 ] *<br>
>  Â  Â (-A_Q28[ 1 ]), 30 );<br>
>  Â  Â  Â  Â  Â  Â  Â  Â S[ 1 ] = silk_SMLAWB( S[ 1 ], B_Q28[ 2 ], in[ k * 2 + 0 ] );<br>
>  Â  Â  Â  Â  Â  Â  Â  Â S[ 3 ] = silk_SMLAWB( S[ 3 ], B_Q28[ 2 ], in[ k * 2 + 1 ] );<br>
><br>
>  Â  Â  Â  Â  Â  Â  Â  Â /* Scale back to Q0 and saturate */<br>
>  Â  Â  Â  Â  Â  Â  Â  Â out[ k * 2 + 0 ] = (opus_int16)silk_SAT16( silk_RSHIFT(<br>
>  Â  Â out32_Q14[ 0 ] + (1<<14) - 1, 14 ) );<br>
>  Â  Â  Â  Â  Â  Â  Â  Â out[ k * 2 + 1 ] = (opus_int16)silk_SAT16( silk_RSHIFT(<br>
>  Â  Â out32_Q14[ 1 ] + (1<<14) - 1, 14 ) );<br>
>  Â  Â  Â  Â  Â  Â }<br>
>  Â  Â  Â  Â }<br>
>  Â  Â }<br>
><br>
>  Â  Â Here is the NEON kernels which uses vqrdmulh_lane_s32() to do the<br>
>  Â  Â multiplication and rounding, where A_Q28_s32x{2,4} stores doubled<br>
>  Â  Â -A_Q28[]:<br>
><br>
>  Â  Â static inline void silk_biquad_alt_stride1_<wbr>kernel(const int32x2_t<br>
>  Â  Â A_Q28_s32x2, const int32x4_t t_s32x4, int32x2_t *S_s32x2, int32x2_t<br>
>  Â  Â *out32_Q14_s32x2)<br>
>  Â  Â {<br>
>  Â  Â  Â  Â int32x2_t t_s32x2;<br>
><br>
>  Â  Â  Â  Â *out32_Q14_s32x2 = vadd_s32(*S_s32x2, vget_low_s32(t_s32x4));<br>
>  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  /* silk_SMLAWB( S[ 0 ], B_Q28[ 0 ], in[ k ]<br>
>  Â  Â )  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  */<br>
>  Â  Â  Â  Â *S_s32x2  Â  Â  Â  Â =<br>
>  Â  Â vreinterpret_s32_u64(vshr_n_<wbr>u64(vreinterpret_u64_s32(*S_<wbr>s32x2),<br>
>  Â  Â 32)); /* S[ 0 ] = S[ 1 ]; S[ 1 ] = 0;<br>
>  Â  Â  Â  Â  Â  Â  Â  Â  Â  */<br>
>  Â  Â  Â  Â *out32_Q14_s32x2 = vshl_n_s32(*out32_Q14_s32x2, 2);<br>
>  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  /* out32_Q14 = silk_LSHIFT( silk_SMLAWB( S[<br>
>  Â  Â 0 ], B_Q28[ 0 ], in[ k ] ), 2 ); */<br>
>  Â  Â  Â  Â t_s32x2  Â  Â  Â  Â  = vqrdmulh_lane_s32(A_Q28_s32x2,<br>
>  Â  Â *out32_Q14_s32x2, 0);  Â  Â  Â  Â  Â  Â  Â  Â  /* silk_RSHIFT_ROUND(<br>
>  Â  Â (opus_int64)out32_Q14 * (-A_Q28[ {0,1} ]), 30 )  Â  Â  Â  */<br>
>  Â  Â  Â  Â *S_s32x2  Â  Â  Â  Â = vadd_s32(*S_s32x2, t_s32x2);<br>
>  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  /* S[ {0,1} ] = {S[ 1 ],0} +<br>
>  Â  Â silk_RSHIFT_ROUND( );  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â */<br>
>  Â  Â  Â  Â *S_s32x2  Â  Â  Â  Â = vadd_s32(*S_s32x2, vget_high_s32(t_s32x4));<br>
>  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â /* S[ {0,1} ] = silk_SMLAWB( S[ {0,1} ],<br>
>  Â  Â B_Q28[ {1,2} ], in[ k ] );  Â  Â  Â  Â  */<br>
>  Â  Â }<br>
><br>
>  Â  Â static inline void silk_biquad_alt_stride2_<wbr>kernel(const int32x4_t<br>
>  Â  Â A_Q28_s32x4, const int32x4_t B_Q28_s32x4, const int32x2_t t_s32x2,<br>
>  Â  Â const int32x4_t inval_s32x4, int32x4_t *S_s32x4, int32x2_t<br>
>  Â  Â *out32_Q14_s32x2)<br>
>  Â  Â {<br>
>  Â  Â  Â  Â int32x4_t t_s32x4, out32_Q14_s32x4;<br>
><br>
>  Â  Â  Â  Â *out32_Q14_s32x2 = vadd_s32(vget_low_s32(*S_<wbr>s32x4), t_s32x2);<br>
>  Â  Â  Â  Â  Â  Â  /* silk_SMLAWB( S{0,1}, B_Q28[ 0 ], in[ k * 2 + {0,1} ] )<br>
>  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â */<br>
>  Â  Â  Â  Â *S_s32x4  Â  Â  Â  Â = vcombine_s32(vget_high_s32(*S_<wbr>s32x4),<br>
>  Â  Â vdup_n_s32(0)); /* S{0,1} = S{2,3}; S{2,3} = 0;<br>
>  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â */<br>
>  Â  Â  Â  Â *out32_Q14_s32x2 = vshl_n_s32(*out32_Q14_s32x2, 2);<br>
>  Â  Â  Â  Â  Â  Â  /* out32_Q14_{0,1} = silk_LSHIFT( silk_SMLAWB( S{0,1},<br>
>  Â  Â B_Q28[ 0 ], in[ k * 2 + {0,1} ] ), 2 );  */<br>
>  Â  Â  Â  Â out32_Q14_s32x4  = vcombine_s32(*out32_Q14_s32x2,<br>
>  Â  Â *out32_Q14_s32x2);  Â  Â /* out32_Q14_{0,1,0,1}<br>
>  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  */<br>
>  Â  Â  Â  Â t_s32x4  Â  Â  Â  Â  = vqrdmulhq_s32(out32_Q14_s32x4, A_Q28_s32x4);<br>
>  Â  Â  Â  Â  Â  Â  /* silk_RSHIFT_ROUND( (opus_int64)out32_Q14[ {0,1,0,1} ] *<br>
>  Â  Â (-A_Q28[ {0,0,1,1} ]), 30 )  Â  Â  Â  Â  */<br>
>  Â  Â  Â  Â *S_s32x4  Â  Â  Â  Â = vaddq_s32(*S_s32x4, t_s32x4);<br>
>  Â  Â  Â  Â  Â  Â /* S[ {0,1,2,3} ] = {S[ {2,3} ],0,0} + silk_RSHIFT_ROUND( );<br>
>  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  */<br>
>  Â  Â  Â  Â t_s32x4  Â  Â  Â  Â  = vqdmulhq_s32(inval_s32x4, B_Q28_s32x4);<br>
>  Â  Â  Â  Â  Â  Â /* silk_SMULWB(B_Q28[ {1,1,2,2} ], in[ k * 2 + {0,1,0,1} ] )<br>
>  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  */<br>
>  Â  Â  Â  Â *S_s32x4  Â  Â  Â  Â = vaddq_s32(*S_s32x4, t_s32x4);<br>
>  Â  Â  Â  Â  Â  Â /* S[ {0,1,2,3} ] = silk_SMLAWB( S[ {0,1,2,3} ], B_Q28[<br>
>  Â  Â {1,1,2,2} ], in[ k * 2 + {0,1,0,1} ] ); */<br>
>  Â  Â }<br>
><br>
>  Â  Â Thanks,<br>
>  Â  Â Linfeng<br>
><br>
><br>
</div></div></blockquote></div><br></div>