/*
  patest_pink.c

  Generate Pink Noise using Gardner method.
  Optimization suggested by James McCartney uses a tree
  to select which random value to replace.

  x x x x x x x x x x x x x x x x 
  x   x   x   x   x   x   x   x   
  x       x       x       x       
  x               x               
  x                               

  Tree is generated by counting trailing zeros in an increasing index.
  When the index is zero, no random number is selected.

  This program uses the Portable Audio library which is under development.
  For more information see:   http://www.audiomulch.com/portaudio/

  Author: Phil Burk, http://www.softsynth.com
	
  Revision History:

  Copyleft 1999 Phil Burk - No rights reserved.
*/

#include <stdio.h>
#include <math.h>

/************************************************************/
/* Calculate pseudo-random 32 bit number based on linear congruential method. */
static unsigned long GenerateRandomNumber( void )
{
    static unsigned long randSeed = 22222;  /* Change this for different random sequences. */
    randSeed = (randSeed * 196314165) + 907633515;
    return randSeed;
}

#define PINK_MAX_RANDOM_ROWS   (30)
#define PINK_RANDOM_BITS       (24)
#define PINK_RANDOM_SHIFT      ((sizeof(long)*8)-PINK_RANDOM_BITS)

typedef struct
{
    long      pink_Rows[PINK_MAX_RANDOM_ROWS];
    long      pink_RunningSum;   /* Used to optimize summing of generators. */
    int       pink_Index;        /* Incremented each sample. */
    int       pink_IndexMask;    /* Index wrapped by ANDing with this mask. */
    float     pink_Scalar;       /* Used to scale within range of -1.0 to +1.0 */
} PinkNoise;

/* Setup PinkNoise structure for N rows of generators. */
void InitializePinkNoise( PinkNoise *pink, int numRows )
{
    int i;
    long pmax;
    pink->pink_Index = 0;
    pink->pink_IndexMask = (1<<numRows) - 1;
/* Calculate maximum possible signed random value. Extra 1 for white noise always added. */
    pmax = (numRows + 1) * (1<<(PINK_RANDOM_BITS-1));
    pink->pink_Scalar = 1.0f / pmax;
/* Initialize rows. */
    for( i=0; i<numRows; i++ ) pink->pink_Rows[i] = 0;
    pink->pink_RunningSum = 0;
}

/* Generate Pink noise values between -1.0 and +1.0 */
float GeneratePinkNoise( PinkNoise *pink )
{
    long newRandom;
    long sum;
    float output;

/* Increment and mask index. */
    pink->pink_Index = (pink->pink_Index + 1) & pink->pink_IndexMask;

/* If index is zero, don't update any random values. */
    if( pink->pink_Index != 0 )
    {
	/* Determine how many trailing zeros in PinkIndex. */
	/* This algorithm will hang if n==0 so test first. */
	int numZeros = 0;
	int n = pink->pink_Index;
	while( (n & 1) == 0 )
	{
	    n = n >> 1;
	    numZeros++;
	}

	/* Replace the indexed ROWS random value.
	 * Subtract and add back to RunningSum instead of adding all the random
	 * values together. Only one changes each time.
	 */
	pink->pink_RunningSum -= pink->pink_Rows[numZeros];
	newRandom = ((long)GenerateRandomNumber()) >> PINK_RANDOM_SHIFT;
	pink->pink_RunningSum += newRandom;
	pink->pink_Rows[numZeros] = newRandom;
    }
	
/* Add extra white noise value. */
    newRandom = ((long)GenerateRandomNumber()) >> PINK_RANDOM_SHIFT;
    sum = pink->pink_RunningSum + newRandom;

/* Scale to range of -1.0 to 0.9999. */
    output = pink->pink_Scalar * sum;

    return output;
}