AudioResamplerDyn.cpp revision 771386e6e6e79697e2d839ef0f25a242946ba1e5
1/* 2 * Copyright (C) 2013 The Android Open Source Project 3 * 4 * Licensed under the Apache License, Version 2.0 (the "License"); 5 * you may not use this file except in compliance with the License. 6 * You may obtain a copy of the License at 7 * 8 * http://www.apache.org/licenses/LICENSE-2.0 9 * 10 * Unless required by applicable law or agreed to in writing, software 11 * distributed under the License is distributed on an "AS IS" BASIS, 12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13 * See the License for the specific language governing permissions and 14 * limitations under the License. 15 */ 16 17#define LOG_TAG "AudioResamplerDyn" 18//#define LOG_NDEBUG 0 19 20#include <malloc.h> 21#include <string.h> 22#include <stdlib.h> 23#include <dlfcn.h> 24#include <math.h> 25 26#include <cutils/compiler.h> 27#include <cutils/properties.h> 28#include <utils/Debug.h> 29#include <utils/Log.h> 30 31#include "AudioResamplerFirOps.h" // USE_NEON and USE_INLINE_ASSEMBLY defined here 32#include "AudioResamplerFirProcess.h" 33#include "AudioResamplerFirProcessNeon.h" 34#include "AudioResamplerFirGen.h" // requires math.h 35#include "AudioResamplerDyn.h" 36 37//#define DEBUG_RESAMPLER 38 39namespace android { 40 41// generate a unique resample type compile-time constant (constexpr) 42#define RESAMPLETYPE(CHANNELS, LOCKED, STRIDE) \ 43 ((((CHANNELS)-1)&1) | !!(LOCKED)<<1 \ 44 | ((STRIDE)==8 ? 1 : (STRIDE)==16 ? 2 : 0)<<2) 45 46/* 47 * InBuffer is a type agnostic input buffer. 48 * 49 * Layout of the state buffer for halfNumCoefs=8. 50 * 51 * [rrrrrrppppppppnnnnnnnnrrrrrrrrrrrrrrrrrrr.... rrrrrrr] 52 * S I R 53 * 54 * S = mState 55 * I = mImpulse 56 * R = mRingFull 57 * p = past samples, convoluted with the (p)ositive side of sinc() 58 * n = future samples, convoluted with the (n)egative side of sinc() 59 * r = extra space for implementing the ring buffer 60 */ 61 62template<typename TC, typename TI, typename TO> 63AudioResamplerDyn<TC, TI, TO>::InBuffer::InBuffer() 64 : mState(NULL), mImpulse(NULL), mRingFull(NULL), mStateCount(0) 65{ 66} 67 68template<typename TC, typename TI, typename TO> 69AudioResamplerDyn<TC, TI, TO>::InBuffer::~InBuffer() 70{ 71 init(); 72} 73 74template<typename TC, typename TI, typename TO> 75void AudioResamplerDyn<TC, TI, TO>::InBuffer::init() 76{ 77 free(mState); 78 mState = NULL; 79 mImpulse = NULL; 80 mRingFull = NULL; 81 mStateCount = 0; 82} 83 84// resizes the state buffer to accommodate the appropriate filter length 85template<typename TC, typename TI, typename TO> 86void AudioResamplerDyn<TC, TI, TO>::InBuffer::resize(int CHANNELS, int halfNumCoefs) 87{ 88 // calculate desired state size 89 int stateCount = halfNumCoefs * CHANNELS * 2 * kStateSizeMultipleOfFilterLength; 90 91 // check if buffer needs resizing 92 if (mState 93 && stateCount == mStateCount 94 && mRingFull-mState == mStateCount-halfNumCoefs*CHANNELS) { 95 return; 96 } 97 98 // create new buffer 99 TI* state; 100 (void)posix_memalign(reinterpret_cast<void**>(&state), 32, stateCount*sizeof(*state)); 101 memset(state, 0, stateCount*sizeof(*state)); 102 103 // attempt to preserve state 104 if (mState) { 105 TI* srcLo = mImpulse - halfNumCoefs*CHANNELS; 106 TI* srcHi = mImpulse + halfNumCoefs*CHANNELS; 107 TI* dst = state; 108 109 if (srcLo < mState) { 110 dst += mState-srcLo; 111 srcLo = mState; 112 } 113 if (srcHi > mState + mStateCount) { 114 srcHi = mState + mStateCount; 115 } 116 memcpy(dst, srcLo, (srcHi - srcLo) * sizeof(*srcLo)); 117 free(mState); 118 } 119 120 // set class member vars 121 mState = state; 122 mStateCount = stateCount; 123 mImpulse = state + halfNumCoefs*CHANNELS; // actually one sample greater than needed 124 mRingFull = state + mStateCount - halfNumCoefs*CHANNELS; 125} 126 127// copy in the input data into the head (impulse+halfNumCoefs) of the buffer. 128template<typename TC, typename TI, typename TO> 129template<int CHANNELS> 130void AudioResamplerDyn<TC, TI, TO>::InBuffer::readAgain(TI*& impulse, const int halfNumCoefs, 131 const TI* const in, const size_t inputIndex) 132{ 133 TI* head = impulse + halfNumCoefs*CHANNELS; 134 for (size_t i=0 ; i<CHANNELS ; i++) { 135 head[i] = in[inputIndex*CHANNELS + i]; 136 } 137} 138 139// advance the impulse pointer, and load in data into the head (impulse+halfNumCoefs) 140template<typename TC, typename TI, typename TO> 141template<int CHANNELS> 142void AudioResamplerDyn<TC, TI, TO>::InBuffer::readAdvance(TI*& impulse, const int halfNumCoefs, 143 const TI* const in, const size_t inputIndex) 144{ 145 impulse += CHANNELS; 146 147 if (CC_UNLIKELY(impulse >= mRingFull)) { 148 const size_t shiftDown = mRingFull - mState - halfNumCoefs*CHANNELS; 149 memcpy(mState, mState+shiftDown, halfNumCoefs*CHANNELS*2*sizeof(TI)); 150 impulse -= shiftDown; 151 } 152 readAgain<CHANNELS>(impulse, halfNumCoefs, in, inputIndex); 153} 154 155template<typename TC, typename TI, typename TO> 156void AudioResamplerDyn<TC, TI, TO>::Constants::set( 157 int L, int halfNumCoefs, int inSampleRate, int outSampleRate) 158{ 159 int bits = 0; 160 int lscale = inSampleRate/outSampleRate < 2 ? L - 1 : 161 static_cast<int>(static_cast<uint64_t>(L)*inSampleRate/outSampleRate); 162 for (int i=lscale; i; ++bits, i>>=1) 163 ; 164 mL = L; 165 mShift = kNumPhaseBits - bits; 166 mHalfNumCoefs = halfNumCoefs; 167} 168 169template<typename TC, typename TI, typename TO> 170AudioResamplerDyn<TC, TI, TO>::AudioResamplerDyn(int bitDepth, 171 int inChannelCount, int32_t sampleRate, src_quality quality) 172 : AudioResampler(bitDepth, inChannelCount, sampleRate, quality), 173 mResampleFunc(0), mFilterSampleRate(0), mFilterQuality(DEFAULT_QUALITY), 174 mCoefBuffer(NULL) 175{ 176 mVolumeSimd[0] = mVolumeSimd[1] = 0; 177 // The AudioResampler base class assumes we are always ready for 1:1 resampling. 178 // We reset mInSampleRate to 0, so setSampleRate() will calculate filters for 179 // setSampleRate() for 1:1. (May be removed if precalculated filters are used.) 180 mInSampleRate = 0; 181 mConstants.set(128, 8, mSampleRate, mSampleRate); // TODO: set better 182} 183 184template<typename TC, typename TI, typename TO> 185AudioResamplerDyn<TC, TI, TO>::~AudioResamplerDyn() 186{ 187 free(mCoefBuffer); 188} 189 190template<typename TC, typename TI, typename TO> 191void AudioResamplerDyn<TC, TI, TO>::init() 192{ 193 mFilterSampleRate = 0; // always trigger new filter generation 194 mInBuffer.init(); 195} 196 197template<typename TC, typename TI, typename TO> 198void AudioResamplerDyn<TC, TI, TO>::setVolume(int16_t left, int16_t right) 199{ 200 AudioResampler::setVolume(left, right); 201 // volume is applied on the output type. 202 if (is_same<TO, float>::value || is_same<TO, double>::value) { 203 const TO scale = 1. / (1UL << 12); 204 mVolumeSimd[0] = static_cast<TO>(left) * scale; 205 mVolumeSimd[1] = static_cast<TO>(right) * scale; 206 } else { 207 mVolumeSimd[0] = static_cast<int32_t>(left) << 16; 208 mVolumeSimd[1] = static_cast<int32_t>(right) << 16; 209 } 210} 211 212template<typename T> T max(T a, T b) {return a > b ? a : b;} 213 214template<typename T> T absdiff(T a, T b) {return a > b ? a - b : b - a;} 215 216template<typename TC, typename TI, typename TO> 217void AudioResamplerDyn<TC, TI, TO>::createKaiserFir(Constants &c, 218 double stopBandAtten, int inSampleRate, int outSampleRate, double tbwCheat) 219{ 220 TC* buf; 221 static const double atten = 0.9998; // to avoid ripple overflow 222 double fcr; 223 double tbw = firKaiserTbw(c.mHalfNumCoefs, stopBandAtten); 224 225 (void)posix_memalign(reinterpret_cast<void**>(&buf), 32, (c.mL+1)*c.mHalfNumCoefs*sizeof(TC)); 226 if (inSampleRate < outSampleRate) { // upsample 227 fcr = max(0.5*tbwCheat - tbw/2, tbw/2); 228 } else { // downsample 229 fcr = max(0.5*tbwCheat*outSampleRate/inSampleRate - tbw/2, tbw/2); 230 } 231 // create and set filter 232 firKaiserGen(buf, c.mL, c.mHalfNumCoefs, stopBandAtten, fcr, atten); 233 c.mFirCoefs = buf; 234 if (mCoefBuffer) { 235 free(mCoefBuffer); 236 } 237 mCoefBuffer = buf; 238#ifdef DEBUG_RESAMPLER 239 // print basic filter stats 240 printf("L:%d hnc:%d stopBandAtten:%lf fcr:%lf atten:%lf tbw:%lf\n", 241 c.mL, c.mHalfNumCoefs, stopBandAtten, fcr, atten, tbw); 242 // test the filter and report results 243 double fp = (fcr - tbw/2)/c.mL; 244 double fs = (fcr + tbw/2)/c.mL; 245 double passMin, passMax, passRipple; 246 double stopMax, stopRipple; 247 testFir(buf, c.mL, c.mHalfNumCoefs, fp, fs, /*passSteps*/ 1000, /*stopSteps*/ 100000, 248 passMin, passMax, passRipple, stopMax, stopRipple); 249 printf("passband(%lf, %lf): %.8lf %.8lf %.8lf\n", 0., fp, passMin, passMax, passRipple); 250 printf("stopband(%lf, %lf): %.8lf %.3lf\n", fs, 0.5, stopMax, stopRipple); 251#endif 252} 253 254// recursive gcd. Using objdump, it appears the tail recursion is converted to a while loop. 255static int gcd(int n, int m) 256{ 257 if (m == 0) { 258 return n; 259 } 260 return gcd(m, n % m); 261} 262 263static bool isClose(int32_t newSampleRate, int32_t prevSampleRate, 264 int32_t filterSampleRate, int32_t outSampleRate) 265{ 266 267 // different upsampling ratios do not need a filter change. 268 if (filterSampleRate != 0 269 && filterSampleRate < outSampleRate 270 && newSampleRate < outSampleRate) 271 return true; 272 273 // check design criteria again if downsampling is detected. 274 int pdiff = absdiff(newSampleRate, prevSampleRate); 275 int adiff = absdiff(newSampleRate, filterSampleRate); 276 277 // allow up to 6% relative change increments. 278 // allow up to 12% absolute change increments (from filter design) 279 return pdiff < prevSampleRate>>4 && adiff < filterSampleRate>>3; 280} 281 282template<typename TC, typename TI, typename TO> 283void AudioResamplerDyn<TC, TI, TO>::setSampleRate(int32_t inSampleRate) 284{ 285 if (mInSampleRate == inSampleRate) { 286 return; 287 } 288 int32_t oldSampleRate = mInSampleRate; 289 int32_t oldHalfNumCoefs = mConstants.mHalfNumCoefs; 290 uint32_t oldPhaseWrapLimit = mConstants.mL << mConstants.mShift; 291 bool useS32 = false; 292 293 mInSampleRate = inSampleRate; 294 295 // TODO: Add precalculated Equiripple filters 296 297 if (mFilterQuality != getQuality() || 298 !isClose(inSampleRate, oldSampleRate, mFilterSampleRate, mSampleRate)) { 299 mFilterSampleRate = inSampleRate; 300 mFilterQuality = getQuality(); 301 302 // Begin Kaiser Filter computation 303 // 304 // The quantization floor for S16 is about 96db - 10*log_10(#length) + 3dB. 305 // Keep the stop band attenuation no greater than 84-85dB for 32 length S16 filters 306 // 307 // For s32 we keep the stop band attenuation at the same as 16b resolution, about 308 // 96-98dB 309 // 310 311 double stopBandAtten; 312 double tbwCheat = 1.; // how much we "cheat" into aliasing 313 int halfLength; 314 if (mFilterQuality == DYN_HIGH_QUALITY) { 315 // 32b coefficients, 64 length 316 useS32 = true; 317 stopBandAtten = 98.; 318 if (inSampleRate >= mSampleRate * 4) { 319 halfLength = 48; 320 } else if (inSampleRate >= mSampleRate * 2) { 321 halfLength = 40; 322 } else { 323 halfLength = 32; 324 } 325 } else if (mFilterQuality == DYN_LOW_QUALITY) { 326 // 16b coefficients, 16-32 length 327 useS32 = false; 328 stopBandAtten = 80.; 329 if (inSampleRate >= mSampleRate * 4) { 330 halfLength = 24; 331 } else if (inSampleRate >= mSampleRate * 2) { 332 halfLength = 16; 333 } else { 334 halfLength = 8; 335 } 336 if (inSampleRate <= mSampleRate) { 337 tbwCheat = 1.05; 338 } else { 339 tbwCheat = 1.03; 340 } 341 } else { // DYN_MED_QUALITY 342 // 16b coefficients, 32-64 length 343 // note: > 64 length filters with 16b coefs can have quantization noise problems 344 useS32 = false; 345 stopBandAtten = 84.; 346 if (inSampleRate >= mSampleRate * 4) { 347 halfLength = 32; 348 } else if (inSampleRate >= mSampleRate * 2) { 349 halfLength = 24; 350 } else { 351 halfLength = 16; 352 } 353 if (inSampleRate <= mSampleRate) { 354 tbwCheat = 1.03; 355 } else { 356 tbwCheat = 1.01; 357 } 358 } 359 360 // determine the number of polyphases in the filterbank. 361 // for 16b, it is desirable to have 2^(16/2) = 256 phases. 362 // https://ccrma.stanford.edu/~jos/resample/Relation_Interpolation_Error_Quantization.html 363 // 364 // We are a bit more lax on this. 365 366 int phases = mSampleRate / gcd(mSampleRate, inSampleRate); 367 368 // TODO: Once dynamic sample rate change is an option, the code below 369 // should be modified to execute only when dynamic sample rate change is enabled. 370 // 371 // as above, #phases less than 63 is too few phases for accurate linear interpolation. 372 // we increase the phases to compensate, but more phases means more memory per 373 // filter and more time to compute the filter. 374 // 375 // if we know that the filter will be used for dynamic sample rate changes, 376 // that would allow us skip this part for fixed sample rate resamplers. 377 // 378 while (phases<63) { 379 phases *= 2; // this code only needed to support dynamic rate changes 380 } 381 382 if (phases>=256) { // too many phases, always interpolate 383 phases = 127; 384 } 385 386 // create the filter 387 mConstants.set(phases, halfLength, inSampleRate, mSampleRate); 388 createKaiserFir(mConstants, stopBandAtten, 389 inSampleRate, mSampleRate, tbwCheat); 390 } // End Kaiser filter 391 392 // update phase and state based on the new filter. 393 const Constants& c(mConstants); 394 mInBuffer.resize(mChannelCount, c.mHalfNumCoefs); 395 const uint32_t phaseWrapLimit = c.mL << c.mShift; 396 // try to preserve as much of the phase fraction as possible for on-the-fly changes 397 mPhaseFraction = static_cast<unsigned long long>(mPhaseFraction) 398 * phaseWrapLimit / oldPhaseWrapLimit; 399 mPhaseFraction %= phaseWrapLimit; // should not do anything, but just in case. 400 mPhaseIncrement = static_cast<uint32_t>(static_cast<double>(phaseWrapLimit) 401 * inSampleRate / mSampleRate); 402 403 // determine which resampler to use 404 // check if locked phase (works only if mPhaseIncrement has no "fractional phase bits") 405 int locked = (mPhaseIncrement << (sizeof(mPhaseIncrement)*8 - c.mShift)) == 0; 406 int stride = (c.mHalfNumCoefs&7)==0 ? 16 : (c.mHalfNumCoefs&3)==0 ? 8 : 2; 407 if (locked) { 408 mPhaseFraction = mPhaseFraction >> c.mShift << c.mShift; // remove fractional phase 409 } 410 411 setResampler(RESAMPLETYPE(mChannelCount, locked, stride)); 412#ifdef DEBUG_RESAMPLER 413 printf("channels:%d %s stride:%d %s coef:%d shift:%d\n", 414 mChannelCount, locked ? "locked" : "interpolated", 415 stride, useS32 ? "S32" : "S16", 2*c.mHalfNumCoefs, c.mShift); 416#endif 417} 418 419template<typename TC, typename TI, typename TO> 420void AudioResamplerDyn<TC, TI, TO>::resample(int32_t* out, size_t outFrameCount, 421 AudioBufferProvider* provider) 422{ 423 (this->*mResampleFunc)(reinterpret_cast<TO*>(out), outFrameCount, provider); 424} 425 426template<typename TC, typename TI, typename TO> 427void AudioResamplerDyn<TC, TI, TO>::setResampler(unsigned resampleType) 428{ 429 // stride 16 (falls back to stride 2 for machines that do not support NEON) 430 switch (resampleType) { 431 case RESAMPLETYPE(1, true, 16): 432 mResampleFunc = &AudioResamplerDyn<TC, TI, TO>::resample<1, true, 16>; 433 return; 434 case RESAMPLETYPE(2, true, 16): 435 mResampleFunc = &AudioResamplerDyn<TC, TI, TO>::resample<2, true, 16>; 436 return; 437 case RESAMPLETYPE(1, false, 16): 438 mResampleFunc = &AudioResamplerDyn<TC, TI, TO>::resample<1, false, 16>; 439 return; 440 case RESAMPLETYPE(2, false, 16): 441 mResampleFunc = &AudioResamplerDyn<TC, TI, TO>::resample<2, false, 16>; 442 return; 443 default: 444 LOG_ALWAYS_FATAL("Invalid resampler type: %u", resampleType); 445 mResampleFunc = NULL; 446 return; 447 } 448} 449 450template<typename TC, typename TI, typename TO> 451template<int CHANNELS, bool LOCKED, int STRIDE> 452void AudioResamplerDyn<TC, TI, TO>::resample(TO* out, size_t outFrameCount, 453 AudioBufferProvider* provider) 454{ 455 const Constants& c(mConstants); 456 const TC* const coefs = mConstants.mFirCoefs; 457 TI* impulse = mInBuffer.getImpulse(); 458 size_t inputIndex = mInputIndex; 459 uint32_t phaseFraction = mPhaseFraction; 460 const uint32_t phaseIncrement = mPhaseIncrement; 461 size_t outputIndex = 0; 462 size_t outputSampleCount = outFrameCount * 2; // stereo output 463 size_t inFrameCount = getInFrameCountRequired(outFrameCount); 464 const uint32_t phaseWrapLimit = c.mL << c.mShift; 465 466 // NOTE: be very careful when modifying the code here. register 467 // pressure is very high and a small change might cause the compiler 468 // to generate far less efficient code. 469 // Always sanity check the result with objdump or test-resample. 470 471 // the following logic is a bit convoluted to keep the main processing loop 472 // as tight as possible with register allocation. 473 while (outputIndex < outputSampleCount) { 474 // buffer is empty, fetch a new one 475 while (mBuffer.frameCount == 0) { 476 mBuffer.frameCount = inFrameCount; 477 provider->getNextBuffer(&mBuffer, 478 calculateOutputPTS(outputIndex / 2)); 479 if (mBuffer.raw == NULL) { 480 goto resample_exit; 481 } 482 if (phaseFraction >= phaseWrapLimit) { // read in data 483 mInBuffer.template readAdvance<CHANNELS>( 484 impulse, c.mHalfNumCoefs, 485 reinterpret_cast<TI*>(mBuffer.raw), inputIndex); 486 phaseFraction -= phaseWrapLimit; 487 while (phaseFraction >= phaseWrapLimit) { 488 inputIndex++; 489 if (inputIndex >= mBuffer.frameCount) { 490 inputIndex -= mBuffer.frameCount; 491 provider->releaseBuffer(&mBuffer); 492 break; 493 } 494 mInBuffer.template readAdvance<CHANNELS>( 495 impulse, c.mHalfNumCoefs, 496 reinterpret_cast<TI*>(mBuffer.raw), inputIndex); 497 phaseFraction -= phaseWrapLimit; 498 } 499 } 500 } 501 const TI* const in = reinterpret_cast<const TI*>(mBuffer.raw); 502 const size_t frameCount = mBuffer.frameCount; 503 const int coefShift = c.mShift; 504 const int halfNumCoefs = c.mHalfNumCoefs; 505 const TO* const volumeSimd = mVolumeSimd; 506 507 // reread the last input in. 508 mInBuffer.template readAgain<CHANNELS>(impulse, halfNumCoefs, in, inputIndex); 509 510 // main processing loop 511 while (CC_LIKELY(outputIndex < outputSampleCount)) { 512 // caution: fir() is inlined and may be large. 513 // output will be loaded with the appropriate values 514 // 515 // from the input samples in impulse[-halfNumCoefs+1]... impulse[halfNumCoefs] 516 // from the polyphase filter of (phaseFraction / phaseWrapLimit) in coefs. 517 // 518 fir<CHANNELS, LOCKED, STRIDE>( 519 &out[outputIndex], 520 phaseFraction, phaseWrapLimit, 521 coefShift, halfNumCoefs, coefs, 522 impulse, volumeSimd); 523 outputIndex += 2; 524 525 phaseFraction += phaseIncrement; 526 while (phaseFraction >= phaseWrapLimit) { 527 inputIndex++; 528 if (inputIndex >= frameCount) { 529 goto done; // need a new buffer 530 } 531 mInBuffer.template readAdvance<CHANNELS>(impulse, halfNumCoefs, in, inputIndex); 532 phaseFraction -= phaseWrapLimit; 533 } 534 } 535done: 536 // often arrives here when input buffer runs out 537 if (inputIndex >= frameCount) { 538 inputIndex -= frameCount; 539 provider->releaseBuffer(&mBuffer); 540 // mBuffer.frameCount MUST be zero here. 541 } 542 } 543 544resample_exit: 545 mInBuffer.setImpulse(impulse); 546 mInputIndex = inputIndex; 547 mPhaseFraction = phaseFraction; 548} 549 550/* instantiate templates used by AudioResampler::create */ 551template class AudioResamplerDyn<float, float, float>; 552template class AudioResamplerDyn<int16_t, int16_t, int32_t>; 553template class AudioResamplerDyn<int32_t, int16_t, int32_t>; 554 555// ---------------------------------------------------------------------------- 556}; // namespace android 557