diff --git a/cwAudioFile.cpp b/cwAudioFile.cpp index 6ddb05a..267c994 100644 --- a/cwAudioFile.cpp +++ b/cwAudioFile.cpp @@ -5,7 +5,7 @@ #include "cwFile.h" #include "cwObject.h" #include "cwAudioFile.h" -#include "cwUtility.h" +#include "cwMath.h" #include "cwFileSys.h" // #define _24to32_aif( p ) ((int)( ((p[0]>127?255:0) << 24) + (((int)p[0]) << 16) + (((int)p[1]) <<8) + p[2])) // no-swap equivalent @@ -179,7 +179,7 @@ namespace cw if((rc = _read(p,s,10,1)) != kOkRC ) return rc; - *x80Ptr = x80ToDouble(s); + *x80Ptr = math::x80ToDouble(s); return kOkRC; } @@ -557,7 +557,7 @@ namespace cw rc_t rc = kOkRC; unsigned char srateX80[10]; - doubleToX80( p->info.srate, srateX80 ); + math::doubleToX80( p->info.srate, srateX80 ); unsigned hdrByteCnt = 54; unsigned ssndByteCnt = 8 + (p->info.chCnt * p->info.frameCnt * (p->info.bits/kBitsPerByte)); diff --git a/cwMath.cpp b/cwMath.cpp index 0ddfae8..9f855d3 100644 --- a/cwMath.cpp +++ b/cwMath.cpp @@ -2,8 +2,403 @@ #include "cwLog.h" #include "cwCommonImpl.h" #include "cwMath.h" +#include "cwMem.h" +#include -unsigned cw::math::randUInt( unsigned minVal, unsigned maxVal ) + +// TODO: rewrite to avoid copying +// this code comes via csound source ... +double cw::math::x80ToDouble( unsigned char rate[10] ) { - return std::max(minVal,std::min(maxVal,minVal + (unsigned)round(((maxVal - minVal) * rand())/RAND_MAX))); + char sign; + short exp = 0; + unsigned long mant1 = 0; + unsigned long mant0 = 0; + double val; + unsigned char* p = (unsigned char*)rate; + + exp = *p++; + exp <<= 8; + exp |= *p++; + sign = (exp & 0x8000) ? 1 : 0; + exp &= 0x7FFF; + + mant1 = *p++; + mant1 <<= 8; + mant1 |= *p++; + mant1 <<= 8; + mant1 |= *p++; + mant1 <<= 8; + mant1 |= *p++; + + mant0 = *p++; + mant0 <<= 8; + mant0 |= *p++; + mant0 <<= 8; + mant0 |= *p++; + mant0 <<= 8; + mant0 |= *p++; + + /* special test for all bits zero meaning zero + - else pow(2,-16383) bombs */ + if (mant1 == 0 && mant0 == 0 && exp == 0 && sign == 0) + return 0.0; + else { + val = ((double)mant0) * pow(2.0,-63.0); + val += ((double)mant1) * pow(2.0,-31.0); + val *= pow(2.0,((double) exp) - 16383.0); + return sign ? -val : val; + } } + +// TODO: rewrite to avoid copying +/* + * Convert double to IEEE 80 bit floating point + * Should be portable to all C compilers. + * 19aug91 aldel/dpwe covered for MSB bug in Ultrix 'cc' + */ + +void cw::math::doubleToX80(double val, unsigned char rate[10]) +{ + char sign = 0; + short exp = 0; + unsigned long mant1 = 0; + unsigned long mant0 = 0; + unsigned char* p = (unsigned char*)rate; + + if (val < 0.0) { sign = 1; val = -val; } + + if (val != 0.0) /* val identically zero -> all elements zero */ + { + exp = (short)(std::log(val)/std::log(2.0) + 16383.0); + val *= pow(2.0, 31.0+16383.0-(double)exp); + mant1 =((unsigned)val); + val -= ((double)mant1); + val *= pow(2.0, 32.0); + mant0 =((double)val); + } + + *p++ = ((sign<<7)|(exp>>8)); + *p++ = (u_char)(0xFF & exp); + *p++ = (u_char)(0xFF & (mant1>>24)); + *p++ = (u_char)(0xFF & (mant1>>16)); + *p++ = (u_char)(0xFF & (mant1>> 8)); + *p++ = (u_char)(0xFF & (mant1)); + *p++ = (u_char)(0xFF & (mant0>>24)); + *p++ = (u_char)(0xFF & (mant0>>16)); + *p++ = (u_char)(0xFF & (mant0>> 8)); + *p++ = (u_char)(0xFF & (mant0)); + +} + +bool cw::math::isPowerOfTwo( unsigned x ) +{ + return !( (x < 2) || (x & (x-1)) ); +} + +unsigned cw::math::nextPowerOfTwo( unsigned val ) +{ + unsigned i; + unsigned mask = 1; + unsigned msb = 0; + unsigned cnt = 0; + + // if val is a power of two return it + if( isPowerOfTwo(val) ) + return val; + + // next pow of zero is 2 + if( val == 0 ) + return 2; + + // if the next power of two can't be represented in 32 bits + if( val > 0x80000000) + { + assert(0); + return 0; + } + + // find most sig. bit that is set - the number with only the next msb set is next pow 2 + for(i=0; i<31; i++,mask<<=1) + if( mask & val ) + { + msb = i; + cnt++; + } + + + return 1 << (msb + 1); +} + +unsigned cw::math::nearPowerOfTwo( unsigned i ) +{ + unsigned vh = nextPowerOfTwo(i); + + if( vh == 2 ) + return vh; + + unsigned vl = vh / 2; + + if( vh - i < i - vl ) + return vh; + return vl; +} + +bool cw::math::isOddU( unsigned v ) { return v % 2 == 1; } +bool cw::math::isEvenU( unsigned v ) { return !isOddU(v); } +unsigned cw::math::nextOddU( unsigned v ) { return isOddU(v) ? v : v+1; } +unsigned cw::math::prevOddU( unsigned v ) { return isOddU(v) ? v : v-1; } +unsigned cw::math::nextEvenU( unsigned v ) { return isEvenU(v) ? v : v+1; } +unsigned cw::math::prevEvenU( unsigned v ) { return isEvenU(v) ? v : v-1; } + +unsigned cw::math::modIncr(int idx, int delta, int maxN ) +{ + int sum = idx + delta; + + if( sum >= maxN ) + return sum - maxN; + + if( sum < 0 ) + return maxN + sum; + + return sum; +} + + +unsigned cw::math::hzToMidi( double hz ) +{ + + float midi = 12.0 * std::log2(hz/13.75) + 9; + + if( midi < 0 ) + midi = 0; + if( midi > 127 ) + midi = 127; + + return (unsigned)lround(midi); +} + +float cw::math::midiToHz( unsigned midi ) +{ + double m = midi <= 127 ? midi : 127; + + return (float)( 13.75 * pow(2.0,(m - 9.0)/12.0)); +} + + + +//================================================================= +// Random numbers + +int cw::math::randInt( int min, int max ) +{ + assert( min <= max ); + int range = max - min; + return min + std::max(0,std::min(range,(int)round(range * (double)rand() / RAND_MAX))); +} + +unsigned cw::math::randUInt( unsigned min, unsigned max ) +{ + assert( min <= max ); + unsigned range = max - min; + unsigned val = (unsigned)round(range * (double)rand() / RAND_MAX); + return min + std::max((unsigned)0,std::min(range,val)); +} + +float cw::math::randFloat( float min, float max ) +{ + assert( min <= max ); + float range = max - min; + float val = (float)(range * (double)rand() / RAND_MAX); + return min + std::max(0.0f,std::min(range,val)); +} + +double cw::math::randDouble( double min, double max ) +{ + assert( min <= max ); + double range = max - min; + double val = range * (double)rand() / RAND_MAX; + return min + std::max(0.0,std::min(range,val)); +} + + +//================================================================= +// Base on: http://stackoverflow.com/questions/3874627/floating-point-comparison-functions-for-c-sharp + +bool cw::math::isCloseD( double x0, double x1, double eps ) +{ + double d = fabs(x0-x1); + + if( x0 == x1 ) + return true; + + if( x0==0 || x1==0 || d x1 ) + return (x0-x1)/(x0+x1) < eps; + else + return (x1-x0)/(x0+x1) < eps; +} + +//================================================================= + +// lFSR() implementation based on note at bottom of: +// http://www.ece.u.edu/~koopman/lfsr/index.html +void cw::math::lFSR( unsigned lfsrN, unsigned tapMask, unsigned seed, unsigned* yV, unsigned yN ) +{ + assert( 0 < lfsrN && lfsrN < 32 ); + + unsigned i; + for(i=0; i> 1) ^ tapMask; + else + seed = (seed >> 1); + + } +} + +namespace cw +{ + namespace math + { + bool mLS_IsBalanced( const unsigned* xV, int xN) + { + int a = 0; + int i; + + for(i=0; i(mlsN); + unsigned* mls1V = mem::allocZ(mlsN); + unsigned* xorV = mem::allocZ(mlsN); + + unsigned i,j; + + lFSR(lfsrN, poly_coeff0, 1 << (lfsrN-1), mls0V, mlsN); + + lFSR(lfsrN, poly_coeff1, 1 << (lfsrN-1), mls1V, mlsN); + + if( mLS_IsBalanced(mls0V,mlsN) ) + yi = _genGoldCopy(yM, yi, yN, mls0V, mlsN); + + if( yi + T bessel0( T x ) + { + T eps = pow(10.0,-9.0); + T n = 1.0; + T S = 1.0; + T D = 1.0; + + while(D > eps*S) + { + T t = x /(2.0*n); + n = n+1; + D = D * pow(t,2.0); + S = S + D; + } + + return S; + + } + + //================================================================= + // The following elliptic-related function approximations come from + // Parks & Burrus, Digital Filter Design, Appendix program 9, pp. 317-326 + // which in turn draws directly on other sources + + // Calculate complete elliptic integral (quarter period) K + // given *complimentary* modulus kc. + template< typename T > + T ellipK( T kc ) + { + T a = 1, b = kc, c = 1, tmp; + + while( c > std::numeric_limits::epsilon ) + { + c = 0.5*(a-b); + tmp = 0.5*(a+b); + b = sqrt(a*b); + a = tmp; + } + + return M_PI/(2*a); + } + + // Calculate elliptic modulus k + // given ratio of complete elliptic integrals r = K/K' + // (solves the "degree equation" for fixed N = K*K1'/K'K1) + template< typename T > + T ellipDeg( T r ) + { + T q,a,b,c,d; + a = b = c = 1; + d = q = exp(-M_PI*r); + + while( c > std::numeric_limits::epsilon ) + { + a = a + 2*c*d; + c = c*d*d; + b = b + c; + d = d*q; + } + + return 4*sqrt(q)*pow(b/a,2); + } + + // calculate arc elliptic tangent u (elliptic integral of the 1st kind) + // given argument x = sc(u,k) and *complimentary* modulus kc + template< typename T > + T ellipArcSc( T x, T kc ) + { + T a = 1, b = kc, y = 1/x, tmp; + unsigned L = 0; + + while( true ) + { + tmp = a*b; + a += b; + b = 2*sqrt(tmp); + y -= tmp/y; + if( y == 0 ) + y = sqrt(tmp) * 1E-10; + if( fabs(a-b)/a < std::numeric_limits::epsilon ) + break; + L *= 2; + if( y < 0 ) + L++; + } + + if( y < 0 ) + L++; + + return (atan(a/y) + M_PI*L)/a; + } + + // calculate Jacobi elliptic functions sn, cn, and dn + // given argument u and *complimentary* modulus kc + template< typename T > + rc_t ellipJ( T u, T kc, T* sn, T* cn, T* dn ) + { + assert( sn != NULL || cn != NULL || dn != NULL ); + + if( u == 0 ) + { + if( sn != NULL ) *sn = 0; + if( cn != NULL ) *cn = 1; + if( dn != NULL ) *dn = 1; + return kOkRC; + } + + int i; + T a,b,c,d,e,tmp,_sn,_cn,_dn; + T aa[16], bb[16]; + + a = 1; + b = kc; + + for( i = 0; i < 16; i++ ) + { + aa[i] = a; + bb[i] = b; + tmp = (a+b)/2; + b = sqrt(a*b); + a = tmp; + if( (a-b)/a < std::numeric_limits::epsilon ) + break; + } + + c = a/tan(u*a); + d = 1; + + for( ; i >= 0; i-- ) + { + e = c*c/a; + c = c*d; + a = aa[i]; + d = (e + bb[i]) / (e+a); + } + + _sn = 1/sqrt(1+c*c); + _cn = _sn*c; + _dn = d; + + if( sn != NULL ) *sn = _sn; + if( cn != NULL ) *cn = _cn; + if( dn != NULL ) *dn = _dn; + + return kOkRC; + } + + //================================================================= + // bilinear transform + // z = (2*sr + s)/(2*sr - s) + template< typename T > + rc_t blt( unsigned n, T sr, T* rp, T* ip ) + { + unsigned i; + T a = 2*sr, + tr, ti, td; + + for( i = 0; i < n; i++ ) + { + tr = rp[i]; + ti = ip[i]; + td = pow(a-tr, 2) + ti*ti; + rp[i] = (a*a - tr*tr - ti*ti)/td; + ip[i] = 2*a*ti/td; + if( tr < -1E15 ) + rp[i] = 0; + if( fabs(ti) > 1E15 ) + ip[i] = 0; + } + + return kOkRC; + } + + + + //================================================================= + // Pitch conversion + unsigned hzToMidi( double hz ); + float midiToHz( unsigned midi ); + + //================================================================= + // Floating point byte swapping + unsigned ffSwapFloatToUInt( float v ); + float ffSwapUIntToFloat( unsigned v ); + unsigned long long ffSwapDoubleToULLong( double v ); + double ffSwapULLongToDouble( unsigned long long v ); + + //================================================================= + template< typename T > + T rand_range(T min, T max ) + { + assert( min <= max ); + T range = max - min; + return min + std::max(0,std::min(range,(T)range * rand() / RAND_MAX)); + + } + + int randInt( int min, int max ); + unsigned randUInt( unsigned min, unsigned max ); + float randFloat( float min, float max ); + double randDouble( double min, double max ); + + //================================================================= + bool isCloseD( double x0, double x1, double eps ); + bool isCloseF( float x0, float x1, double eps ); + bool isCloseI( int x0, int x1, double eps ); + bool isCloseU( unsigned x0, unsigned x1, double eps ); + + //================================================================= + // Run a length 'lfsrN' linear feedback shift register (LFSR) for 'yN' iterations to + // produce a length 'yN' bit string in yV[yN]. + // 'lfsrN' count of bits in the shift register range: 2<= lfsrN <= 32. + // 'tapMask' is a bit mask which gives the tap indexes positions for the LFSR. + // The least significant bit corresponds to the maximum delay tap position. + // The min tap position is therefore denoted by the tap mask bit location 1 << (lfsrN-1). + // A minimum of two taps must exist. + // 'seed' sets the initial delay state. + // 'yV[yN]' is the the output vector + // 'yN' is count of elements in yV. + // The function resturn kOkAtRC on success or kInvalidArgsRCRC if any arguments are invalid. + // /sa lFSR_Test. + void lFSR( unsigned lfsrN, unsigned tapMask, unsigned seed, unsigned* yV, unsigned yN ); + + // Example and test code for lFSR() + bool lFSR_Test(); + + + // Generate a set of 'goldN' Gold codes using the Maximum Length Sequences (MLS) generated + // by a length 'lfsrN' linear feedback shift register. + // 'err' is an error object to be set if the the function fails. + // 'lfsrN' is the length of the Linear Feedback Shift Registers (LFSR) used to generate the MLS. + // 'poly_coeff0' tap mask for the first LFSR. + // 'coeff1' tap mask the the second LFSR. + // 'goldN' is the count of Gold codes to generate. + // 'yM[mlsN', goldN] is a column major output matrix where each column contains a Gold code. + // 'mlsN' is the length of the maximum length sequence for each Gold code which can be + // calculated as mlsN = (1 << a->lfsrN) - 1. + // Note that values of 'lfsrN' and the 'poly_coeffx' must be carefully selected such that + // they will produce a MLS. For example to generate a MLS with length 31 set 'lfsrN' to 5 and + // then select poly_coeff from two different elements of the set {0x12 0x14 0x17 0x1B 0x1D 0x1E}. + // See http://www.ece.u.edu/~koopman/lfsr/index.html for a complete set of MSL polynomial + // coefficients for given LFSR lengths. + // Returns false if insufficient balanced pairs exist. + bool genGoldCodes( unsigned lfsrN, unsigned poly_coeff0, unsigned poly_coeff1, unsigned goldN, int* yM, unsigned mlsN ); + + } } diff --git a/cwUtility.cpp b/cwUtility.cpp index 1597bae..505f0cb 100644 --- a/cwUtility.cpp +++ b/cwUtility.cpp @@ -36,6 +36,7 @@ void cw::printHex( const void* buf, unsigned bufByteN, bool asciiFl ) } } +#ifdef NOT_DEF // TODO: rewrite to avoid copying // this code comes via csound source ... @@ -174,3 +175,5 @@ unsigned cw::nearestPowerOfTwo( unsigned i ) return vh; return vl; } + +#endif