michael@0: /* Copyright (C) 2007-2008 Jean-Marc Valin
michael@0:  * Copyright (C) 2008 Thorvald Natvig
michael@0:  */
michael@0: /**
michael@0:    @file resample_sse.h
michael@0:    @brief Resampler functions (SSE version)
michael@0: */
michael@0: /*
michael@0:    Redistribution and use in source and binary forms, with or without
michael@0:    modification, are permitted provided that the following conditions
michael@0:    are met:
michael@0:    
michael@0:    - Redistributions of source code must retain the above copyright
michael@0:    notice, this list of conditions and the following disclaimer.
michael@0:    
michael@0:    - Redistributions in binary form must reproduce the above copyright
michael@0:    notice, this list of conditions and the following disclaimer in the
michael@0:    documentation and/or other materials provided with the distribution.
michael@0:    
michael@0:    - Neither the name of the Xiph.org Foundation nor the names of its
michael@0:    contributors may be used to endorse or promote products derived from
michael@0:    this software without specific prior written permission.
michael@0:    
michael@0:    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
michael@0:    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
michael@0:    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
michael@0:    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
michael@0:    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
michael@0:    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
michael@0:    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
michael@0:    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
michael@0:    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
michael@0:    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
michael@0:    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
michael@0: */
michael@0: 
michael@0: #include <xmmintrin.h>
michael@0: 
michael@0: #define OVERRIDE_INNER_PRODUCT_SINGLE
michael@0: static inline float inner_product_single(const float *a, const float *b, unsigned int len)
michael@0: {
michael@0:    int i;
michael@0:    float ret;
michael@0:    if (1)
michael@0:    {
michael@0:       __m128 sum = _mm_setzero_ps();
michael@0:       for (i=0;i<len;i+=8)
michael@0:       {
michael@0:          sum = _mm_add_ps(sum, _mm_mul_ps(_mm_loadu_ps(a+i), _mm_loadu_ps(b+i)));
michael@0:          sum = _mm_add_ps(sum, _mm_mul_ps(_mm_loadu_ps(a+i+4), _mm_loadu_ps(b+i+4)));
michael@0:       }
michael@0:       sum = _mm_add_ps(sum, _mm_movehl_ps(sum, sum));
michael@0:       sum = _mm_add_ss(sum, _mm_shuffle_ps(sum, sum, 0x55));
michael@0:       _mm_store_ss(&ret, sum);
michael@0:    }
michael@0:    else
michael@0:    {
michael@0:       ret = 0;
michael@0:       for (i=0;i<len;i++) ret += a[i] * b[i];
michael@0:    }
michael@0:    return ret;
michael@0: }
michael@0: 
michael@0: #define OVERRIDE_INTERPOLATE_PRODUCT_SINGLE
michael@0: static inline float interpolate_product_single(const float *a, const float *b, unsigned int len, const spx_uint32_t oversample, float *frac) {
michael@0:   int i;
michael@0:   float ret;
michael@0:   if (1)
michael@0:   {
michael@0:     __m128 sum = _mm_setzero_ps();
michael@0:     __m128 f = _mm_loadu_ps(frac);
michael@0:     for(i=0;i<len;i+=2)
michael@0:     {
michael@0:       sum = _mm_add_ps(sum, _mm_mul_ps(_mm_load1_ps(a+i), _mm_loadu_ps(b+i*oversample)));
michael@0:       sum = _mm_add_ps(sum, _mm_mul_ps(_mm_load1_ps(a+i+1), _mm_loadu_ps(b+(i+1)*oversample)));
michael@0:     }
michael@0:     sum = _mm_mul_ps(f, sum);
michael@0:     sum = _mm_add_ps(sum, _mm_movehl_ps(sum, sum));
michael@0:     sum = _mm_add_ss(sum, _mm_shuffle_ps(sum, sum, 0x55));
michael@0:     _mm_store_ss(&ret, sum);
michael@0:   }
michael@0:   else
michael@0:   {
michael@0:     float accum[4] = {0,0,0,0};
michael@0:     for(i=0;i<len;i++)
michael@0:     {
michael@0:       const float curr_in=a[i];
michael@0:       accum[0] += curr_in * b[i * oversample + 0];
michael@0:       accum[1] += curr_in * b[i * oversample + 1];
michael@0:       accum[2] += curr_in * b[i * oversample + 2];
michael@0:       accum[3] += curr_in * b[i * oversample + 3];
michael@0:     }
michael@0:     ret = accum[0] * frac[0] + accum[1] * frac[1] + accum[2] * frac[2] + accum[3] * frac[3];
michael@0:   }
michael@0:   return ret;
michael@0: }
michael@0: 
michael@0: #ifdef __SSE2__
michael@0: #include <emmintrin.h>
michael@0: #define OVERRIDE_INNER_PRODUCT_DOUBLE
michael@0: 
michael@0: static inline double inner_product_double(const float *a, const float *b, unsigned int len)
michael@0: {
michael@0:    int i;
michael@0:    double ret;
michael@0:    __m128d sum = _mm_setzero_pd();
michael@0:    __m128 t;
michael@0:    for (i=0;i<len;i+=8)
michael@0:    {
michael@0:       t = _mm_mul_ps(_mm_loadu_ps(a+i), _mm_loadu_ps(b+i));
michael@0:       sum = _mm_add_pd(sum, _mm_cvtps_pd(t));
michael@0:       sum = _mm_add_pd(sum, _mm_cvtps_pd(_mm_movehl_ps(t, t)));
michael@0: 
michael@0:       t = _mm_mul_ps(_mm_loadu_ps(a+i+4), _mm_loadu_ps(b+i+4));
michael@0:       sum = _mm_add_pd(sum, _mm_cvtps_pd(t));
michael@0:       sum = _mm_add_pd(sum, _mm_cvtps_pd(_mm_movehl_ps(t, t)));
michael@0:    }
michael@0:    sum = _mm_add_sd(sum, (__m128d) _mm_movehl_ps((__m128) sum, (__m128) sum));
michael@0:    _mm_store_sd(&ret, sum);
michael@0:    return ret;
michael@0: }
michael@0: 
michael@0: #define OVERRIDE_INTERPOLATE_PRODUCT_DOUBLE
michael@0: static inline double interpolate_product_double(const float *a, const float *b, unsigned int len, const spx_uint32_t oversample, float *frac) {
michael@0:   int i;
michael@0:   double ret;
michael@0:   __m128d sum;
michael@0:   __m128d sum1 = _mm_setzero_pd();
michael@0:   __m128d sum2 = _mm_setzero_pd();
michael@0:   __m128 f = _mm_loadu_ps(frac);
michael@0:   __m128d f1 = _mm_cvtps_pd(f);
michael@0:   __m128d f2 = _mm_cvtps_pd(_mm_movehl_ps(f,f));
michael@0:   __m128 t;
michael@0:   for(i=0;i<len;i+=2)
michael@0:   {
michael@0:     t = _mm_mul_ps(_mm_load1_ps(a+i), _mm_loadu_ps(b+i*oversample));
michael@0:     sum1 = _mm_add_pd(sum1, _mm_cvtps_pd(t));
michael@0:     sum2 = _mm_add_pd(sum2, _mm_cvtps_pd(_mm_movehl_ps(t, t)));
michael@0: 
michael@0:     t = _mm_mul_ps(_mm_load1_ps(a+i+1), _mm_loadu_ps(b+(i+1)*oversample));
michael@0:     sum1 = _mm_add_pd(sum1, _mm_cvtps_pd(t));
michael@0:     sum2 = _mm_add_pd(sum2, _mm_cvtps_pd(_mm_movehl_ps(t, t)));
michael@0:   }
michael@0:   sum1 = _mm_mul_pd(f1, sum1);
michael@0:   sum2 = _mm_mul_pd(f2, sum2);
michael@0:   sum = _mm_add_pd(sum1, sum2);
michael@0:   sum = _mm_add_sd(sum, (__m128d) _mm_movehl_ps((__m128) sum, (__m128) sum));
michael@0:   _mm_store_sd(&ret, sum);
michael@0:   return ret;
michael@0: }
michael@0: 
michael@0: #endif