media/libsoundtouch/src/mmx_optimized.cpp

Thu, 22 Jan 2015 13:21:57 +0100

author
Michael Schloh von Bennewitz <michael@schloh.com>
date
Thu, 22 Jan 2015 13:21:57 +0100
branch
TOR_BUG_9701
changeset 15
b8a032363ba2
permissions
-rw-r--r--

Incorporate requested changes from Mozilla in review:
https://bugzilla.mozilla.org/show_bug.cgi?id=1123480#c6

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

mercurial