https://ccrma.stanford.edu/~jos/resample/Theory_Ideal_Bandlimited_Interpolation.html

static inline float interpolate_product_single_align(const float *a,
const float *b, unsigned int len, const spx_uint32_t oversample, float *frac)
{
    /* calculate the offset of array starting point 
       (1) when overlap is 12 bytes, we need start at 16-byte boundary, so offset is always 0      
       address=>|  0  |  1  |  2  |  3  |
       offset =>|  0  |  0  |  0  |  0  |
       (1) when overlap is 8 bytes, we cannot peel to 16-byte boundar, so set offset to 1 or 0 for later special handing
       address=>|  0  |  1  |  2  |  3  |
       offset =>|  0  |  1  |  0  |  1  |
     */
    int i=0,offset=(oversample==2?((16-(((guint32)b)&15))>>2)%2:0);
    float ret;
    __m128 sum = _mm_setzero_ps();
    __m128 f = _mm_load_ps(frac);    
    /* remove SSE misalignment in the worst cases (half of) 
       namely each iteration will have overlap part as 12 bytes or 8 bytes 
       (1) when overlap is 1
       |  i  |  i  |  i  |  i  | ... |
       |     | i+1 | i+1 | i+1 | i+1 |
       (2) when overlap is 2
       |  i  |  i  |  i  |  i  | ... | ... |
       | ... | ... | i+1 | i+1 | i+1 | i+1 |
     */
    if (oversample==1||oversample==2) {
        __m128 p, q, r, s;
        /* peeling the leading items, make sure array "b" begins at the 16-byte boundary only when overlap is 12 bytes*/
        for(i=0; i<len&&((guint32)(b+i*oversample))&15!=offset; i++) {
            sum = _mm_add_ps(sum, _mm_mul_ps(_mm_load1_ps(a+i), _mm_loadu_ps(b+i*oversample)));
        } 
        /* special handing for first scenario when overlap is 12 bytes */
        if (oversample==1) {
            /* each iteration can calculate 4 values */
            for (r=_mm_load_ps(b+i*oversample); i+4<len; i+=4) {
                /* calculate first sum, load constant s for this iteration and pre-load constant p for the next value */
                s=_mm_load_ps(b+(i+4)*oversample);
                p = _mm_move_ss(r, s);
                p = _mm_shuffle_ps(p, p, _MM_SHUFFLE(0,3,2,1));
                sum = _mm_add_ps(sum, _mm_mul_ps(_mm_load1_ps(a+i), r));
        
                /* calculate second sum and pre-load constant q for the next value */
                q = _mm_shuffle_ps(p, s, _MM_SHUFFLE(1,0,2,1));
                sum = _mm_add_ps(sum, _mm_mul_ps(_mm_load1_ps(a+i+1), p));
        
                /* calculate third sum and pre-load constant p for the next value */
                p = _mm_move_ss(p, _mm_shuffle_ps(s, s, _MM_SHUFFLE(1,0,3,2)));
                p = _mm_shuffle_ps(p, p, _MM_SHUFFLE(0,3,2,1));
                sum = _mm_add_ps(sum, _mm_mul_ps(_mm_load1_ps(a+i+2), q));
        
                /* calculate fourth sum and set constant r as s for the next iteration */
                r = s;
                sum = _mm_add_ps(sum, _mm_mul_ps(_mm_load1_ps(a+i+3), p));
            }
        }
        /* special handing for part of second scenario when overlap is 8 bytes and start from even address, namely 0 and 2 (offset is 0) */
        else if (offset==0) {
            for(r=_mm_load_ps(b+i*oversample);i+2<len;i+=2) {
                s = _mm_load_ps(b+(i+2)*oversample);                
                p = _mm_shuffle_ps(r, s, _MM_SHUFFLE(1,0,3,2));
                sum = _mm_add_ps(sum, _mm_mul_ps(_mm_load1_ps(a+i), r));

                sum = _mm_add_ps(sum, _mm_mul_ps(_mm_load1_ps(a+i+1), p));
                r = s;
            }   
        }
        /* special handing for part of second scenario when overlap is 8 bytes and start from odd address, namely 1 and 3 (offset is 1) */
        else {
            for(r=_mm_load_ps(b+i*oversample-1);i+2<len;i+=2) {
                s = _mm_load_ps(b+(i+2)*oversample-1);
                p = _mm_move_ss(r, s);
                p = _mm_shuffle_ps(p, p, _MM_SHUFFLE(0,3,2,1));
                sum = _mm_add_ps(sum, _mm_mul_ps(_mm_load1_ps(a+i), p));
        
                q = _mm_move_ss(_mm_shuffle_ps(s, s, _MM_SHUFFLE(2,1,0,3)), _mm_shuffle_ps(r, r, _MM_SHUFFLE(2,1,0,3)));
                sum = _mm_add_ps(sum, _mm_mul_ps(_mm_load1_ps(a+i+1), q));
                r = s;
            }
        }
        /* handle remaining items */
		for(;i<len;i+=2) {
		    sum = _mm_add_ps(sum, _mm_mul_ps(_mm_load1_ps(a+i), _mm_loadu_ps(b+i*oversample)));
		    sum = _mm_add_ps(sum, _mm_mul_ps(_mm_load1_ps(a+i+1), _mm_loadu_ps(b+(i+1)*oversample)));
		}            
    }    
    /* not optimize for remaining scenarios and calculate the tailing values if there is any (not devided by 4) */
    for(i=0;i<len;i+=2) {
        sum = _mm_add_ps(sum, _mm_mul_ps(_mm_load1_ps(a+i), _mm_loadu_ps(b+i*oversample)));
        sum = _mm_add_ps(sum, _mm_mul_ps(_mm_load1_ps(a+i+1), _mm_loadu_ps(b+(i+1)*oversample)));
    }    
    /* calcuate the return value, collapse result vector to single scalar result */
    sum = _mm_mul_ps(f, sum);
    sum = _mm_add_ps(sum, _mm_movehl_ps(sum, sum));
    sum = _mm_add_ss(sum, _mm_shuffle_ps(sum, sum, 0x55));
    _mm_store_ss(&ret, sum);
    return ret;
}