~noskcaj/ubuntu/wily/nootka/1.2

« back to all changes in this revision

Viewing changes to src/libs/sound/stouch/mmx_optimized.cpp

  • Committer: Jackson Doak
  • Date: 2015-07-03 05:03:13 UTC
  • mfrom: (1.1.3)
  • Revision ID: noskcaj@ubuntu.com-20150703050313-q07jgvw6ltsosx83
* New upstream release
* debian/control:
  - Switch build-deps to Qt5
  - Run wrap-and-sort -sa
  - Bump standards-version to 3.9.6 (no changes)
  - Bump debhelper version to 9

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
////////////////////////////////////////////////////////////////////////////////
 
2
///
 
3
/// MMX optimized routines. All MMX optimized functions have been gathered into 
 
4
/// this single source code file, regardless to their class or original source 
 
5
/// code file, in order to ease porting the library to other compiler and 
 
6
/// processor platforms.
 
7
///
 
8
/// The MMX-optimizations are programmed using MMX compiler intrinsics that
 
9
/// are supported both by Microsoft Visual C++ and GCC compilers, so this file
 
10
/// should compile with both toolsets.
 
11
///
 
12
/// NOTICE: If using Visual Studio 6.0, you'll need to install the "Visual C++ 
 
13
/// 6.0 processor pack" update to support compiler intrinsic syntax. The update
 
14
/// is available for download at Microsoft Developers Network, see here:
 
15
/// http://msdn.microsoft.com/en-us/vstudio/aa718349.aspx
 
16
///
 
17
/// Author        : Copyright (c) Olli Parviainen
 
18
/// Author e-mail : oparviai 'at' iki.fi
 
19
/// SoundTouch WWW: http://www.surina.net/soundtouch
 
20
///
 
21
////////////////////////////////////////////////////////////////////////////////
 
22
//
 
23
// Last changed  : $Date: 2014-01-07 20:25:40 +0200 (Tue, 07 Jan 2014) $
 
24
// File revision : $Revision: 4 $
 
25
//
 
26
// $Id: mmx_optimized.cpp 184 2014-01-07 18:25:40Z oparviai $
 
27
//
 
28
////////////////////////////////////////////////////////////////////////////////
 
29
//
 
30
// License :
 
31
//
 
32
//  SoundTouch audio processing library
 
33
//  Copyright (c) Olli Parviainen
 
34
//
 
35
//  This library is free software; you can redistribute it and/or
 
36
//  modify it under the terms of the GNU Lesser General Public
 
37
//  License as published by the Free Software Foundation; either
 
38
//  version 2.1 of the License, or (at your option) any later version.
 
39
//
 
40
//  This library is distributed in the hope that it will be useful,
 
41
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
42
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
43
//  Lesser General Public License for more details.
 
44
//
 
45
//  You should have received a copy of the GNU Lesser General Public
 
46
//  License along with this library; if not, write to the Free Software
 
47
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
48
//
 
49
////////////////////////////////////////////////////////////////////////////////
 
50
 
 
51
#include "STTypes.h"
 
52
 
 
53
#ifdef SOUNDTOUCH_ALLOW_MMX
 
54
// MMX routines available only with integer sample type
 
55
 
 
56
using namespace soundtouch;
 
57
 
 
58
//////////////////////////////////////////////////////////////////////////////
 
59
//
 
60
// implementation of MMX optimized functions of class 'TDStretchMMX'
 
61
//
 
62
//////////////////////////////////////////////////////////////////////////////
 
63
 
 
64
#include "TDStretch.h"
 
65
#include <mmintrin.h>
 
66
#include <limits.h>
 
67
#include <math.h>
 
68
 
 
69
 
 
70
// Calculates cross correlation of two buffers
 
71
double TDStretchMMX::calcCrossCorr(const short *pV1, const short *pV2, double &dnorm) const
 
72
{
 
73
    const __m64 *pVec1, *pVec2;
 
74
    __m64 shifter;
 
75
    __m64 accu, normaccu;
 
76
    long corr, norm;
 
77
    int i;
 
78
   
 
79
    pVec1 = (__m64*)pV1;
 
80
    pVec2 = (__m64*)pV2;
 
81
 
 
82
    shifter = _m_from_int(overlapDividerBits);
 
83
    normaccu = accu = _mm_setzero_si64();
 
84
 
 
85
    // Process 4 parallel sets of 2 * stereo samples or 4 * mono samples 
 
86
    // during each round for improved CPU-level parallellization.
 
87
    for (i = 0; i < channels * overlapLength / 16; i ++)
 
88
    {
 
89
        __m64 temp, temp2;
 
90
 
 
91
        // dictionary of instructions:
 
92
        // _m_pmaddwd   : 4*16bit multiply-add, resulting two 32bits = [a0*b0+a1*b1 ; a2*b2+a3*b3]
 
93
        // _mm_add_pi32 : 2*32bit add
 
94
        // _m_psrad     : 32bit right-shift
 
95
 
 
96
        temp = _mm_add_pi32(_mm_sra_pi32(_mm_madd_pi16(pVec1[0], pVec2[0]), shifter),
 
97
                            _mm_sra_pi32(_mm_madd_pi16(pVec1[1], pVec2[1]), shifter));
 
98
        temp2 = _mm_add_pi32(_mm_sra_pi32(_mm_madd_pi16(pVec1[0], pVec1[0]), shifter),
 
99
                            _mm_sra_pi32(_mm_madd_pi16(pVec1[1], pVec1[1]), shifter));
 
100
        accu = _mm_add_pi32(accu, temp);
 
101
        normaccu = _mm_add_pi32(normaccu, temp2);
 
102
 
 
103
        temp = _mm_add_pi32(_mm_sra_pi32(_mm_madd_pi16(pVec1[2], pVec2[2]), shifter),
 
104
                            _mm_sra_pi32(_mm_madd_pi16(pVec1[3], pVec2[3]), shifter));
 
105
        temp2 = _mm_add_pi32(_mm_sra_pi32(_mm_madd_pi16(pVec1[2], pVec1[2]), shifter),
 
106
                            _mm_sra_pi32(_mm_madd_pi16(pVec1[3], pVec1[3]), shifter));
 
107
        accu = _mm_add_pi32(accu, temp);
 
108
        normaccu = _mm_add_pi32(normaccu, temp2);
 
109
 
 
110
        pVec1 += 4;
 
111
        pVec2 += 4;
 
112
    }
 
113
 
 
114
    // copy hi-dword of mm0 to lo-dword of mm1, then sum mmo+mm1
 
115
    // and finally store the result into the variable "corr"
 
116
 
 
117
    accu = _mm_add_pi32(accu, _mm_srli_si64(accu, 32));
 
118
    corr = _m_to_int(accu);
 
119
 
 
120
    normaccu = _mm_add_pi32(normaccu, _mm_srli_si64(normaccu, 32));
 
121
    norm = _m_to_int(normaccu);
 
122
 
 
123
    // Clear MMS state
 
124
    _m_empty();
 
125
 
 
126
    // Normalize result by dividing by sqrt(norm) - this step is easiest 
 
127
    // done using floating point operation
 
128
    dnorm = (double)norm;
 
129
 
 
130
    return (double)corr / sqrt(dnorm < 1e-9 ? 1.0 : dnorm);
 
131
    // Note: Warning about the missing EMMS instruction is harmless
 
132
    // as it'll be called elsewhere.
 
133
}
 
134
 
 
135
 
 
136
/// Update cross-correlation by accumulating "norm" coefficient by previously calculated value
 
137
double TDStretchMMX::calcCrossCorrAccumulate(const short *pV1, const short *pV2, double &dnorm) const
 
138
{
 
139
    const __m64 *pVec1, *pVec2;
 
140
    __m64 shifter;
 
141
    __m64 accu;
 
142
    long corr, lnorm;
 
143
    int i;
 
144
   
 
145
    // cancel first normalizer tap from previous round
 
146
    lnorm = 0;
 
147
    for (i = 1; i <= channels; i ++)
 
148
    {
 
149
        lnorm -= (pV1[-i] * pV1[-i]) >> overlapDividerBits;
 
150
    }
 
151
 
 
152
    pVec1 = (__m64*)pV1;
 
153
    pVec2 = (__m64*)pV2;
 
154
 
 
155
    shifter = _m_from_int(overlapDividerBits);
 
156
    accu = _mm_setzero_si64();
 
157
 
 
158
    // Process 4 parallel sets of 2 * stereo samples or 4 * mono samples 
 
159
    // during each round for improved CPU-level parallellization.
 
160
    for (i = 0; i < channels * overlapLength / 16; i ++)
 
161
    {
 
162
        __m64 temp;
 
163
 
 
164
        // dictionary of instructions:
 
165
        // _m_pmaddwd   : 4*16bit multiply-add, resulting two 32bits = [a0*b0+a1*b1 ; a2*b2+a3*b3]
 
166
        // _mm_add_pi32 : 2*32bit add
 
167
        // _m_psrad     : 32bit right-shift
 
168
 
 
169
        temp = _mm_add_pi32(_mm_sra_pi32(_mm_madd_pi16(pVec1[0], pVec2[0]), shifter),
 
170
                            _mm_sra_pi32(_mm_madd_pi16(pVec1[1], pVec2[1]), shifter));
 
171
        accu = _mm_add_pi32(accu, temp);
 
172
 
 
173
        temp = _mm_add_pi32(_mm_sra_pi32(_mm_madd_pi16(pVec1[2], pVec2[2]), shifter),
 
174
                            _mm_sra_pi32(_mm_madd_pi16(pVec1[3], pVec2[3]), shifter));
 
175
        accu = _mm_add_pi32(accu, temp);
 
176
 
 
177
        pVec1 += 4;
 
178
        pVec2 += 4;
 
179
    }
 
180
 
 
181
    // copy hi-dword of mm0 to lo-dword of mm1, then sum mmo+mm1
 
182
    // and finally store the result into the variable "corr"
 
183
 
 
184
    accu = _mm_add_pi32(accu, _mm_srli_si64(accu, 32));
 
185
    corr = _m_to_int(accu);
 
186
 
 
187
    // Clear MMS state
 
188
    _m_empty();
 
189
 
 
190
    // update normalizer with last samples of this round
 
191
    pV1 = (short *)pVec1;
 
192
    for (int j = 1; j <= channels; j ++)
 
193
    {
 
194
        lnorm += (pV1[-j] * pV1[-j]) >> overlapDividerBits;
 
195
    }
 
196
    dnorm += (double)lnorm;
 
197
 
 
198
    // Normalize result by dividing by sqrt(norm) - this step is easiest 
 
199
    // done using floating point operation
 
200
    return (double)corr / sqrt((dnorm < 1e-9) ? 1.0 : dnorm);
 
201
}
 
202
 
 
203
 
 
204
void TDStretchMMX::clearCrossCorrState()
 
205
{
 
206
    // Clear MMS state
 
207
    _m_empty();
 
208
    //_asm EMMS;
 
209
}
 
210
 
 
211
 
 
212
 
 
213
// MMX-optimized version of the function overlapStereo
 
214
void TDStretchMMX::overlapStereo(short *output, const short *input) const
 
215
{
 
216
    const __m64 *pVinput, *pVMidBuf;
 
217
    __m64 *pVdest;
 
218
    __m64 mix1, mix2, adder, shifter;
 
219
    int i;
 
220
 
 
221
    pVinput  = (const __m64*)input;
 
222
    pVMidBuf = (const __m64*)pMidBuffer;
 
223
    pVdest   = (__m64*)output;
 
224
 
 
225
    // mix1  = mixer values for 1st stereo sample
 
226
    // mix1  = mixer values for 2nd stereo sample
 
227
    // adder = adder for updating mixer values after each round
 
228
    
 
229
    mix1  = _mm_set_pi16(0, overlapLength,   0, overlapLength);
 
230
    adder = _mm_set_pi16(1, -1, 1, -1);
 
231
    mix2  = _mm_add_pi16(mix1, adder);
 
232
    adder = _mm_add_pi16(adder, adder);
 
233
 
 
234
    // Overlaplength-division by shifter. "+1" is to account for "-1" deduced in
 
235
    // overlapDividerBits calculation earlier.
 
236
    shifter = _m_from_int(overlapDividerBits + 1);
 
237
 
 
238
    for (i = 0; i < overlapLength / 4; i ++)
 
239
    {
 
240
        __m64 temp1, temp2;
 
241
                
 
242
        // load & shuffle data so that input & mixbuffer data samples are paired
 
243
        temp1 = _mm_unpacklo_pi16(pVMidBuf[0], pVinput[0]);     // = i0l m0l i0r m0r
 
244
        temp2 = _mm_unpackhi_pi16(pVMidBuf[0], pVinput[0]);     // = i1l m1l i1r m1r
 
245
 
 
246
        // temp = (temp .* mix) >> shifter
 
247
        temp1 = _mm_sra_pi32(_mm_madd_pi16(temp1, mix1), shifter);
 
248
        temp2 = _mm_sra_pi32(_mm_madd_pi16(temp2, mix2), shifter);
 
249
        pVdest[0] = _mm_packs_pi32(temp1, temp2); // pack 2*2*32bit => 4*16bit
 
250
 
 
251
        // update mix += adder
 
252
        mix1 = _mm_add_pi16(mix1, adder);
 
253
        mix2 = _mm_add_pi16(mix2, adder);
 
254
 
 
255
        // --- second round begins here ---
 
256
 
 
257
        // load & shuffle data so that input & mixbuffer data samples are paired
 
258
        temp1 = _mm_unpacklo_pi16(pVMidBuf[1], pVinput[1]);       // = i2l m2l i2r m2r
 
259
        temp2 = _mm_unpackhi_pi16(pVMidBuf[1], pVinput[1]);       // = i3l m3l i3r m3r
 
260
 
 
261
        // temp = (temp .* mix) >> shifter
 
262
        temp1 = _mm_sra_pi32(_mm_madd_pi16(temp1, mix1), shifter);
 
263
        temp2 = _mm_sra_pi32(_mm_madd_pi16(temp2, mix2), shifter);
 
264
        pVdest[1] = _mm_packs_pi32(temp1, temp2); // pack 2*2*32bit => 4*16bit
 
265
 
 
266
        // update mix += adder
 
267
        mix1 = _mm_add_pi16(mix1, adder);
 
268
        mix2 = _mm_add_pi16(mix2, adder);
 
269
 
 
270
        pVinput  += 2;
 
271
        pVMidBuf += 2;
 
272
        pVdest   += 2;
 
273
    }
 
274
 
 
275
    _m_empty(); // clear MMS state
 
276
}
 
277
 
 
278
 
 
279
//////////////////////////////////////////////////////////////////////////////
 
280
//
 
281
// implementation of MMX optimized functions of class 'FIRFilter'
 
282
//
 
283
//////////////////////////////////////////////////////////////////////////////
 
284
 
 
285
#include "FIRFilter.h"
 
286
 
 
287
 
 
288
FIRFilterMMX::FIRFilterMMX() : FIRFilter()
 
289
{
 
290
    filterCoeffsUnalign = NULL;
 
291
}
 
292
 
 
293
 
 
294
FIRFilterMMX::~FIRFilterMMX()
 
295
{
 
296
    delete[] filterCoeffsUnalign;
 
297
}
 
298
 
 
299
 
 
300
// (overloaded) Calculates filter coefficients for MMX routine
 
301
void FIRFilterMMX::setCoefficients(const short *coeffs, uint newLength, uint uResultDivFactor)
 
302
{
 
303
    uint i;
 
304
    FIRFilter::setCoefficients(coeffs, newLength, uResultDivFactor);
 
305
 
 
306
    // Ensure that filter coeffs array is aligned to 16-byte boundary
 
307
    delete[] filterCoeffsUnalign;
 
308
    filterCoeffsUnalign = new short[2 * newLength + 8];
 
309
    filterCoeffsAlign = (short *)SOUNDTOUCH_ALIGN_POINTER_16(filterCoeffsUnalign);
 
310
 
 
311
    // rearrange the filter coefficients for mmx routines 
 
312
    for (i = 0;i < length; i += 4) 
 
313
    {
 
314
        filterCoeffsAlign[2 * i + 0] = coeffs[i + 0];
 
315
        filterCoeffsAlign[2 * i + 1] = coeffs[i + 2];
 
316
        filterCoeffsAlign[2 * i + 2] = coeffs[i + 0];
 
317
        filterCoeffsAlign[2 * i + 3] = coeffs[i + 2];
 
318
 
 
319
        filterCoeffsAlign[2 * i + 4] = coeffs[i + 1];
 
320
        filterCoeffsAlign[2 * i + 5] = coeffs[i + 3];
 
321
        filterCoeffsAlign[2 * i + 6] = coeffs[i + 1];
 
322
        filterCoeffsAlign[2 * i + 7] = coeffs[i + 3];
 
323
    }
 
324
}
 
325
 
 
326
 
 
327
 
 
328
// mmx-optimized version of the filter routine for stereo sound
 
329
uint FIRFilterMMX::evaluateFilterStereo(short *dest, const short *src, uint numSamples) const
 
330
{
 
331
    // Create stack copies of the needed member variables for asm routines :
 
332
    uint i, j;
 
333
    __m64 *pVdest = (__m64*)dest;
 
334
 
 
335
    if (length < 2) return 0;
 
336
 
 
337
    for (i = 0; i < (numSamples - length) / 2; i ++)
 
338
    {
 
339
        __m64 accu1;
 
340
        __m64 accu2;
 
341
        const __m64 *pVsrc = (const __m64*)src;
 
342
        const __m64 *pVfilter = (const __m64*)filterCoeffsAlign;
 
343
 
 
344
        accu1 = accu2 = _mm_setzero_si64();
 
345
        for (j = 0; j < lengthDiv8 * 2; j ++)
 
346
        {
 
347
            __m64 temp1, temp2;
 
348
 
 
349
            temp1 = _mm_unpacklo_pi16(pVsrc[0], pVsrc[1]);  // = l2 l0 r2 r0
 
350
            temp2 = _mm_unpackhi_pi16(pVsrc[0], pVsrc[1]);  // = l3 l1 r3 r1
 
351
 
 
352
            accu1 = _mm_add_pi32(accu1, _mm_madd_pi16(temp1, pVfilter[0]));  // += l2*f2+l0*f0 r2*f2+r0*f0
 
353
            accu1 = _mm_add_pi32(accu1, _mm_madd_pi16(temp2, pVfilter[1]));  // += l3*f3+l1*f1 r3*f3+r1*f1
 
354
 
 
355
            temp1 = _mm_unpacklo_pi16(pVsrc[1], pVsrc[2]);  // = l4 l2 r4 r2
 
356
 
 
357
            accu2 = _mm_add_pi32(accu2, _mm_madd_pi16(temp2, pVfilter[0]));  // += l3*f2+l1*f0 r3*f2+r1*f0
 
358
            accu2 = _mm_add_pi32(accu2, _mm_madd_pi16(temp1, pVfilter[1]));  // += l4*f3+l2*f1 r4*f3+r2*f1
 
359
 
 
360
            // accu1 += l2*f2+l0*f0 r2*f2+r0*f0
 
361
            //       += l3*f3+l1*f1 r3*f3+r1*f1
 
362
 
 
363
            // accu2 += l3*f2+l1*f0 r3*f2+r1*f0
 
364
            //          l4*f3+l2*f1 r4*f3+r2*f1
 
365
 
 
366
            pVfilter += 2;
 
367
            pVsrc += 2;
 
368
        }
 
369
        // accu >>= resultDivFactor
 
370
        accu1 = _mm_srai_pi32(accu1, resultDivFactor);
 
371
        accu2 = _mm_srai_pi32(accu2, resultDivFactor);
 
372
 
 
373
        // pack 2*2*32bits => 4*16 bits
 
374
        pVdest[0] = _mm_packs_pi32(accu1, accu2);
 
375
        src += 4;
 
376
        pVdest ++;
 
377
    }
 
378
 
 
379
   _m_empty();  // clear emms state
 
380
 
 
381
    return (numSamples & 0xfffffffe) - length;
 
382
}
 
383
 
 
384
#endif  // SOUNDTOUCH_ALLOW_MMX