michael@0: ////////////////////////////////////////////////////////////////////////////////
michael@0: /// 
michael@0: /// Sample interpolation routine using 8-tap band-limited Shannon interpolation 
michael@0: /// with kaiser window.
michael@0: ///
michael@0: /// Notice. This algorithm is remarkably much heavier than linear or cubic
michael@0: /// interpolation, and not remarkably better than cubic algorithm. Thus mostly
michael@0: /// for experimental purposes
michael@0: ///
michael@0: /// Author        : Copyright (c) Olli Parviainen
michael@0: /// Author e-mail : oparviai 'at' iki.fi
michael@0: /// SoundTouch WWW: http://www.surina.net/soundtouch
michael@0: ///
michael@0: ////////////////////////////////////////////////////////////////////////////////
michael@0: //
michael@0: // $Id: InterpolateShannon.cpp 195 2014-04-06 15:57:21Z oparviai $
michael@0: //
michael@0: ////////////////////////////////////////////////////////////////////////////////
michael@0: //
michael@0: // License :
michael@0: //
michael@0: //  SoundTouch audio processing library
michael@0: //  Copyright (c) Olli Parviainen
michael@0: //
michael@0: //  This library is free software; you can redistribute it and/or
michael@0: //  modify it under the terms of the GNU Lesser General Public
michael@0: //  License as published by the Free Software Foundation; either
michael@0: //  version 2.1 of the License, or (at your option) any later version.
michael@0: //
michael@0: //  This library is distributed in the hope that it will be useful,
michael@0: //  but WITHOUT ANY WARRANTY; without even the implied warranty of
michael@0: //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
michael@0: //  Lesser General Public License for more details.
michael@0: //
michael@0: //  You should have received a copy of the GNU Lesser General Public
michael@0: //  License along with this library; if not, write to the Free Software
michael@0: //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
michael@0: //
michael@0: ////////////////////////////////////////////////////////////////////////////////
michael@0: 
michael@0: #include <math.h>
michael@0: #include "InterpolateShannon.h"
michael@0: #include "STTypes.h"
michael@0: 
michael@0: using namespace soundtouch;
michael@0: 
michael@0: 
michael@0: /// Kaiser window with beta = 2.0
michael@0: /// Values scaled down by 5% to avoid overflows
michael@0: static const double _kaiser8[8] = 
michael@0: {
michael@0:    0.41778693317814,
michael@0:    0.64888025049173,
michael@0:    0.83508562409944,
michael@0:    0.93887857733412,
michael@0:    0.93887857733412,
michael@0:    0.83508562409944,
michael@0:    0.64888025049173,
michael@0:    0.41778693317814
michael@0: };
michael@0: 
michael@0: 
michael@0: InterpolateShannon::InterpolateShannon()
michael@0: {
michael@0:     fract = 0;
michael@0: }
michael@0: 
michael@0: 
michael@0: void InterpolateShannon::resetRegisters()
michael@0: {
michael@0:     fract = 0;
michael@0: }
michael@0: 
michael@0: 
michael@0: #define PI 3.1415926536
michael@0: #define sinc(x) (sin(PI * (x)) / (PI * (x)))
michael@0: 
michael@0: /// Transpose mono audio. Returns number of produced output samples, and 
michael@0: /// updates "srcSamples" to amount of consumed source samples
michael@0: int InterpolateShannon::transposeMono(SAMPLETYPE *pdest, 
michael@0:                     const SAMPLETYPE *psrc, 
michael@0:                     int &srcSamples)
michael@0: {
michael@0:     int i;
michael@0:     int srcSampleEnd = srcSamples - 8;
michael@0:     int srcCount = 0;
michael@0: 
michael@0:     i = 0;
michael@0:     while (srcCount < srcSampleEnd)
michael@0:     {
michael@0:         double out;
michael@0:         assert(fract < 1.0);
michael@0: 
michael@0:         out  = psrc[0] * sinc(-3.0 - fract) * _kaiser8[0];
michael@0:         out += psrc[1] * sinc(-2.0 - fract) * _kaiser8[1];
michael@0:         out += psrc[2] * sinc(-1.0 - fract) * _kaiser8[2];
michael@0:         if (fract < 1e-6)
michael@0:         {
michael@0:             out += psrc[3] * _kaiser8[3];     // sinc(0) = 1
michael@0:         }
michael@0:         else
michael@0:         {
michael@0:             out += psrc[3] * sinc(- fract) * _kaiser8[3];
michael@0:         }
michael@0:         out += psrc[4] * sinc( 1.0 - fract) * _kaiser8[4];
michael@0:         out += psrc[5] * sinc( 2.0 - fract) * _kaiser8[5];
michael@0:         out += psrc[6] * sinc( 3.0 - fract) * _kaiser8[6];
michael@0:         out += psrc[7] * sinc( 4.0 - fract) * _kaiser8[7];
michael@0: 
michael@0:         pdest[i] = (SAMPLETYPE)out;
michael@0:         i ++;
michael@0: 
michael@0:         // update position fraction
michael@0:         fract += rate;
michael@0:         // update whole positions
michael@0:         int whole = (int)fract;
michael@0:         fract -= whole;
michael@0:         psrc += whole;
michael@0:         srcCount += whole;
michael@0:     }
michael@0:     srcSamples = srcCount;
michael@0:     return i;
michael@0: }
michael@0: 
michael@0: 
michael@0: /// Transpose stereo audio. Returns number of produced output samples, and 
michael@0: /// updates "srcSamples" to amount of consumed source samples
michael@0: int InterpolateShannon::transposeStereo(SAMPLETYPE *pdest, 
michael@0:                     const SAMPLETYPE *psrc, 
michael@0:                     int &srcSamples)
michael@0: {
michael@0:     int i;
michael@0:     int srcSampleEnd = srcSamples - 8;
michael@0:     int srcCount = 0;
michael@0: 
michael@0:     i = 0;
michael@0:     while (srcCount < srcSampleEnd)
michael@0:     {
michael@0:         double out0, out1, w;
michael@0:         assert(fract < 1.0);
michael@0: 
michael@0:         w = sinc(-3.0 - fract) * _kaiser8[0];
michael@0:         out0 = psrc[0] * w; out1 = psrc[1] * w;
michael@0:         w = sinc(-2.0 - fract) * _kaiser8[1];
michael@0:         out0 += psrc[2] * w; out1 += psrc[3] * w;
michael@0:         w = sinc(-1.0 - fract) * _kaiser8[2];
michael@0:         out0 += psrc[4] * w; out1 += psrc[5] * w;
michael@0:         w = _kaiser8[3] * ((fract < 1e-5) ? 1.0 : sinc(- fract));   // sinc(0) = 1
michael@0:         out0 += psrc[6] * w; out1 += psrc[7] * w;
michael@0:         w = sinc( 1.0 - fract) * _kaiser8[4];
michael@0:         out0 += psrc[8] * w; out1 += psrc[9] * w;
michael@0:         w = sinc( 2.0 - fract) * _kaiser8[5];
michael@0:         out0 += psrc[10] * w; out1 += psrc[11] * w;
michael@0:         w = sinc( 3.0 - fract) * _kaiser8[6];
michael@0:         out0 += psrc[12] * w; out1 += psrc[13] * w;
michael@0:         w = sinc( 4.0 - fract) * _kaiser8[7];
michael@0:         out0 += psrc[14] * w; out1 += psrc[15] * w;
michael@0: 
michael@0:         pdest[2*i]   = (SAMPLETYPE)out0;
michael@0:         pdest[2*i+1] = (SAMPLETYPE)out1;
michael@0:         i ++;
michael@0: 
michael@0:         // update position fraction
michael@0:         fract += rate;
michael@0:         // update whole positions
michael@0:         int whole = (int)fract;
michael@0:         fract -= whole;
michael@0:         psrc += 2*whole;
michael@0:         srcCount += whole;
michael@0:     }
michael@0:     srcSamples = srcCount;
michael@0:     return i;
michael@0: }
michael@0: 
michael@0: 
michael@0: /// Transpose stereo audio. Returns number of produced output samples, and 
michael@0: /// updates "srcSamples" to amount of consumed source samples
michael@0: int InterpolateShannon::transposeMulti(SAMPLETYPE *pdest, 
michael@0:                     const SAMPLETYPE *psrc, 
michael@0:                     int &srcSamples)
michael@0: {
michael@0:     // not implemented
michael@0:     assert(false);
michael@0:     return 0;
michael@0: }