Direct pink noise synthesis with auto-correlated generator

notes
Canonical C++ class with minimum system dependencies, BUT you must provide your own
uniform random number generator. Accurate range is a little over 9 octaves, degrading
gracefully beyond this. Estimated deviations +-0.25 dB from ideal 1/f curve in range.
Scaled to fit signed 16-bit range.
code
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
// Pink noise class using the autocorrelated generator method.
// Method proposed and described by Larry Trammell "the RidgeRat" --
// see http://home.earthlink.net/~ltrammell/tech/newpink.htm
// There are no restrictions.
//
// ------------------------------------------------------------------
//
// This is a canonical, 16-bit fixed-point implementation of the
// generator in 32-bit arithmetic. There are only a few system
// dependencies.
//
//   -- access to an allocator 'malloc' for operator new
//   -- access to definition of 'size_t'
//   -- assumes 32-bit two's complement arithmetic
//   -- assumes long int is 32 bits, short int is 16 bits
//   -- assumes that signed right shift propagates the sign bit
//
// It needs a separate URand class to provide uniform 16-bit random
// numbers on interval [1,65535]. The assumed class must provide
// methods to query and set the current seed value, establish a
// scrambled initial seed value, and evaluate uniform random values.
//
//
// ----------- header -----------------------------------------------
// pinkgen.h

#ifndef  _pinkgen_h_
#define  _pinkgen_h_  1

#include  <stddef.h>
#include  <alloc.h>

// You must provide the uniform random generator class.
#ifndef  _URand_h_
#include  "URand.h"
#endif

class PinkNoise {
  private:
    // Coefficients (fixed)
    static long int const pA[5];
    static short int const pPSUM[5];

    // Internal pink generator state
    long int   contrib[5];   // stage contributions
    long int   accum;        // combined generators
    void       internal_clear( );

    // Include a UNoise component
    URand     ugen;

  public:
    PinkNoise( );
    PinkNoise( PinkNoise & );
    ~PinkNoise( );
    void *  operator new( size_t );
    void  pinkclear( );
    short int  pinkrand( );
} ;
#endif

// ----------- implementation ---------------------------------------
// pinkgen.cpp

#include  "pinkgen.h"

// Static class data
long int const PinkNoise::pA[5] =
    { 14055, 12759, 10733, 12273, 15716 };
short int const PinkNoise::pPSUM[5] =
    { 22347, 27917, 29523, 29942, 30007 };

// Clear generator to a zero state.
void   PinkNoise::pinkclear( )
{
    int  i;
    for  (i=0; i<5; ++i)  { contrib[i]=0L; }
    accum = 0L;
}

// PRIVATE, clear generator and also scramble the internal
// uniform generator seed.
void   PinkNoise::internal_clear( )
{
    pinkclear();
    ugen.seed(0);    // Randomizes the seed!
}

// Constructor. Guarantee that initial state is cleared
// and uniform generator scrambled.
PinkNoise::PinkNoise( )
{
    internal_clear();
}

// Copy constructor. Preserve generator state from the source
// object, including the uniform generator seed.
PinkNoise::PinkNoise( PinkNoise & Source )
{
    int  i;
    for (i=0; i<5; ++i)  contrib[i]=Source.contrib[i];
    accum = Source.accum;
    ugen.seed( Source.ugen.seed( ) );
}

// Operator new. Just fetch required object storage.
void *   PinkNoise::operator new( size_t size )
{
    return   malloc(size);
}

// Destructor. No special action required.
PinkNoise::~PinkNoise( )  { /* NIL */ }

// Coding artifact for convenience
#define   UPDATE_CONTRIB(n)  \
     {                                   \
       accum -= contrib[n];              \
       contrib[n] = (long)randv * pA[n]; \
       accum += contrib[n];              \
       break;                            \
     }

// Evaluate next randomized 'pink' number with uniform CPU loading.
short int   PinkNoise::pinkrand( )
{
    short int  randu = ugen.urand() & 0x7fff;     // U[0,32767]
    short int  randv = (short int) ugen.urand();  // U[-32768,32767]

    // Structured block, at most one update is performed
    while (1)
    {
      if (randu < pPSUM[0]) UPDATE_CONTRIB(0);
      if (randu < pPSUM[1]) UPDATE_CONTRIB(1);
      if (randu < pPSUM[2]) UPDATE_CONTRIB(2);
      if (randu < pPSUM[3]) UPDATE_CONTRIB(3);
      if (randu < pPSUM[4]) UPDATE_CONTRIB(4);
      break;
    }
    return (short int) (accum >> 16);
}

// ----------- application   -------------------------------

short int    pink_signal[1024];

void   example(void)
{
    PinkNoise    pinkgen;
    int  i;
    for  (i=0; i<1024; ++i)   pink_signal[i] = pinkgen.pinkrand();
}