cwMath.h/cpp : Populate cwMath from libcm and move some functions from cwUtility to cwMath.

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));

cwMath.cpp

@@ -2,8 +2,403 @@
 #include "cwLog.h"
 #include "cwCommonImpl.h"
 #include "cwMath.h"
+#include "cwMem.h"
+#include <algorithm>
 
-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<DBL_MIN )
+    return d < (eps * DBL_MIN);
+
+  return (d / std::min( fabs(x0) + fabs(x1), DBL_MAX)) < eps;
+}
+
+bool cw::math::isCloseF( float  x0, float  x1, double  eps_d )
+{
+  float eps = (float)eps_d;
+  float d = fabsf(x0-x1);
+  
+  if( x0 == x1 )
+    return true;
+
+  if( x0==0 || x1==0 || d<FLT_MIN )
+    return d < (eps * FLT_MIN);
+
+  return (d / std::min( fabsf(x0) + fabsf(x1), FLT_MAX)) < eps;
+}
+
+bool cw::math::isCloseI( int x0, int x1, double eps )
+{
+  if( x0 == x1 )
+    return true;
+  
+  return abs(x0-x1)/(abs(x0)+abs(x1)) < eps;
+}
+
+
+bool cw::math::isCloseU( unsigned x0, unsigned x1, double eps )
+{
+  if( x0 == x1 )
+    return true;
+  if( x0 > 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<yN; ++i)
+  {
+    if( (yV[i] = seed & 1)==1 )
+      seed = (seed >> 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<xN; ++i)
+        if( xV[i] == 1 )
+          ++a;
+      
+      return abs(a - (xN-a)) == 1;
+    }
+  }
+
+  unsigned _genGoldCopy( int* y, unsigned yi, unsigned yN, unsigned* x, unsigned xN)
+  {
+    unsigned i;
+    for(i=0; i<xN; ++i,++yi)
+      y[yi] = x[i]==1 ? -1 : 1;
+
+    assert(yi <= yN);
+    return yi;
+  }
+  
+}
+
+
+bool cw::math::genGoldCodes( unsigned lfsrN, unsigned poly_coeff0, unsigned poly_coeff1, unsigned goldN, int* yM, unsigned mlsN  )
+{
+  bool      retFl = true;
+  unsigned  yi    = 0;
+  unsigned  yN    = goldN * mlsN;
+  unsigned* mls0V = mem::allocZ<unsigned>(mlsN);
+  unsigned* mls1V = mem::allocZ<unsigned>(mlsN);
+  unsigned* xorV  = mem::allocZ<unsigned>(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<yN && mLS_IsBalanced(mls1V,mlsN) )
+    yi = _genGoldCopy(yM, yi, yN, mls1V, mlsN);
+
+  
+  for(i=0;  yi < yN && i<mlsN-1; ++i )
+  {
+    for(j=0; j<mlsN; ++j)
+      xorV[j] = (mls0V[j] + mls1V[ (i+j) % mlsN ]) % 2;
+    
+    if( mLS_IsBalanced(xorV,mlsN) )
+      yi = _genGoldCopy(yM,yi,yN,xorV,mlsN);
+  }
+
+  if(yi < yN )
+  {    
+    //rc = errMsg(err,kOpFailAtRC,"Gold code generation failed.  Insuffient balanced pairs.");
+    retFl = false;
+  }
+  
+  mem::release(mls0V);
+  mem::release(mls1V);
+  mem::release(xorV);
+
+  return retFl;
+
+}
+
+bool  cw::math::lFSR_Test()
+{
+  // lfsrN          = 5;   % 5    6    7;
+  // poly_coeff0    = 0x12;  % 0x12 0x21 0x41;
+  // poly_coeff1    = 0x1e;  % 0x1e 0x36 0x72;
+
+  unsigned lfsrN = 7;
+  unsigned pc0   = 0x41;
+  unsigned pc1   = 0x72;
+  unsigned mlsN    = (1 << lfsrN)-1;
+
+  unsigned yN = mlsN*2;
+  unsigned yV[ yN ];
+  unsigned i;
+
+  lFSR( lfsrN, pc0, 1 << (lfsrN-1), yV, yN );
+
+  for(i=0; i<mlsN; ++i)
+    if( yV[i] != yV[i+mlsN] )
+      return false;
+
+  //atVOU_PrintL(NULL,"0x12",yV,mlsN,2);
+
+  lFSR( lfsrN, pc1, 1 << (lfsrN-1), yV, yN );
+
+  //atVOU_PrintL(NULL,"0x17",yV,mlsN,2);
+
+  for(i=0; i<mlsN; ++i)
+    if( yV[i] != yV[i+mlsN] )
+      return false;
+
+  return true;
+}
+
+
+
+

cwMath.h

@@ -5,7 +5,277 @@ namespace cw
 {
   namespace math
   {
-    unsigned randUInt(unsigned minVal, unsigned maxVal );
+
+    double   x80ToDouble( unsigned char s[10] );
+    void     doubleToX80( double v, unsigned char s[10] );
+
+    bool     isPowerOfTwo(   unsigned i );
+    unsigned nextPowerOfTwo( unsigned i );
+    unsigned nearPowerOfTwo( unsigned i );
+
+    bool     isOddU(    unsigned v );
+    bool     isEvenU(   unsigned v );
+    unsigned nextOddU(  unsigned v );
+    unsigned prevOddU(  unsigned v );
+    unsigned nextEvenU( unsigned v );
+    unsigned prevEvenU( unsigned v );
+
+    /// Increment or decrement 'idx' by 'delta' always wrapping the result into the range
+    /// 0 to (maxN-1).
+    /// 'idx': initial value 
+    /// 'delta':  incremental amount
+    /// 'maxN' - 1 : maximum return value.
+    unsigned modIncr(int idx, int delta, int maxN );
+
+    // modified bessel function of first kind, order 0
+    // ref: orfandis appendix B io.m
+    template< typename T >
+    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<T>::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<T>::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<T>::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<T>::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  );
+
+    
     
   }    
 }

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