// A fast implementation of Larry Trammell's pink-noise generator. // Written by Andreas Robinson in 2007. // This code is placed in the public domain. /*** References ***/ // White noise generation: // http://en.wikipedia.org/wiki/LFSR // Pink noise generation: // http://www.ridgerat-tech.us/pink/pinkalg.htm // http://www.ridgerat-tech.us/pink/newpink.htm // Declare htons() and htonl(). #ifdef WIN32 #define _CRT_SECURE_NO_WARNINGS #include // Remember to link with Ws2_32.lib #else // GNU #include #endif #include #include #include /*** Utility functions ***/ // Create a vector of noise samples using given noise_fn() short int* GetNoise(int length, double (*noise_fn)()) { int i; short int* buf; buf = malloc(sizeof(short int) * length); if (buf == NULL) { fprintf(stderr, "Out of memory! Exiting...\n"); exit(1); } for (i = 0; i < length; i++) { buf[i] = (short int)(noise_fn() * 32000.0); } return buf; } // Write a 32-bit value to file, in big-endian order. size_t Write32(FILE* file, int v) { v = htonl(v); return fwrite(&v, 1, 4, file); } // Write a 16-bit value to file, in big-endian order. size_t Write16(FILE* file, short int v) { v = htons(v); return fwrite(&v, 1, 2, file); } // Write 16-bit samples in data[] to a .au sound file. void WriteAU(const char* filename, int fs, short int* data, int length) { int i; FILE* file; file = fopen(filename, "wb"); if (!file) { printf("Could not open %s: %s. Exiting.\n", filename, strerror(errno)); exit(1); } Write32(file, 0x2e736e64); // Magic number Write32(file, 24); // Header size, in bytes. Write32(file, 0xffffffff); // Sound length (samples) Write32(file, 3); // Encoding. 3 = 16-bit PCM Write32(file, fs); // Sample rate Write32(file, 1); // Number of channels for (i = 0; i < length; i++) { Write16(file, data[i]); } fclose(file); } /*** The pink-noise generator. ***/ unsigned int lfsr; // Generate a pseudo-random 0 or 1. int Rand1() { int x = lfsr & 1; lfsr = (lfsr >> 1) ^ (-(signed int)(lfsr & 1) & 0xd0000001u); return x; } // Generate a random value between 0 and 32767 (inclusive) int Rand2() { int i; int x = 0; // You may shorten this loop as needed by the pSum comparison. for (i = 0; i < 15; i++) { x = (x << 1) | Rand1(); } return x; } // Generate a random double value between 0.0 and 1.0 (inclusive) double Rand() { return ((double) Rand2()) / 32767.0; } // PRNG test double GetWhiteValue() { int i; int x = 0; for (i = 0; i < 16; i++) { x = (x << 1) | Rand1(); } return 2*(((double)x / 65535) - 0.5); } // Original pink-noise algorithm by Larry Trammell. double contrib[5]; double genout; void Setup0() { int i; lfsr = 1; for (i = 0; i < 5; i++) contrib[i] = 0; } double GetPinkValue0() { const double pA[5] = { 3.80240, 2.96940, 2.59700, 3.08700, 3.40060 }; const double pSum[5] = { 0.00198, 0.01478, 0.06378, 0.23378, 0.91578 }; const double pASum = 15.8564; int i; double ur1 = Rand(); double ur2 = Rand(); for (i = 0; i < 5; i++) { if (ur1 < pSum[i]) { genout -= contrib[i]; contrib[i] = 2 * (ur2 - 0.5) * pA[i]; genout += contrib[i]; break; } } return genout / pASum; } // The modified algorithm double gensum[32]; int gsidx; void Setup1() { const double pA[] = { 3.80240, 2.96940, 2.59700, 3.08700, 3.40060 }; double pASum = 15.8564; int i; lfsr = 1; // Build lookup-table for (i = 0; i < 32; i++) { gensum[i] = ((((i & 1) != 0) ? pA[0] : 0) + (((i & 2) != 0) ? pA[1] : 0) + (((i & 4) != 0) ? pA[2] : 0) + (((i & 8) != 0) ? pA[3] : 0) + (((i & 16) != 0) ? pA[4] : 0)); gensum[i] = (gensum[i] / pASum) - 0.5; } gsidx = 0; } double GetPinkValue1() { // Floating point pSum: // const double pSum[] = { 0.00198, 0.01478, 0.06378, 0.23378, 0.91578 }; // double ur1 = Rand(); // 13-bit pSum: const int pSum[] = { 16, 121, 522, 1915, 7501 }; int ur1 = Rand2() >> 2; //ur1 < 8192 // 8-bit pSum: // const int pSum[] = { 1, 4, 16, 60, 234 }; // int ur1 = Rand2() >> 7; //ur1 < 256 int i; int ur2 = Rand1(); for (i = 0; i < 5; i++) { if (ur1 < pSum[i]) { gsidx &= ~(1 << i); gsidx |= ur2 << i; break; } } return gensum[gsidx]; } /* Hardware implementation ======================= The pA[] lookup-table/summing can be replaced with a custom 5-bit DAC, like this: R0 U0 o--XXXX--*--*---o Uout R1 | | U1 o--XXXX--* X | R2 | X Rs | Is U2 o--XXXX--* X v R3 | _|_ U3 o--XXXX--* GND R4 | U4 o--XXXX-- Ik ------> where Uk, k = 0..4 are the bits in "gsidx". Is, Ik = currents, defined for node analysis below. Node analysis ------------- Uout = Uk - Rk*Ik => Ik = (Uk - Uout) / Rk Is = Sum[Ik] = Uout / Rs => Sum[(Uk - Uout) / Rk] = Uout / Rs <=> Sum[Rs * (Uk - Uout) / Rk] = Uout <=> Sum[(Rs/Rk) * (Uk - Uout)] = Uout <=> Sum[(Rs/Rk) * Uk] - Uout*Sum[Rs/Rk] = Uout <=> Sum[(Rs/Rk) * Uk] = Uout + Uout*Sum[Rs/Rk] <=> Sum[(Rs/Rk) * Uk] = Uout * (1 + Sum[Rs/Rk]) Choose Rs/Rk = pA[k] / pASum => Sum[pA[k]/pASum * Uk] = Uout * (1 + pA[k]/pASum) Because Sum[pA[k]/pASum] = 1 => Sum[pA[k]/pASum * Uk] = 2 * Uout which, besides the factor 2, gives us the same function used in the pink-noise algorithm. Resistor calculation -------------------- It is given that Rs/Rk = pA[k] / pASum, k = 0..4 => Rk = Rs / (pA[k] / pASum) Choose Rs freely and then calculate Rk e.g with CalcResistors(). Results ------- Results for Rs = 10 kohm, first column matching pA[] exactly. Exact Nearest E24 Nearest E96 resistance resistance error resistance error R0 41701.030928 39.0k + 2.7k 0.00% 42.2k 1.20% R1 53399.339934 51.0k + 2.4k 0.00% 53.6k 0.38% R2 61056.603774 51.0k + 10.0k 0.09% 60.4k -1.07% R3 51365.079365 51.0k + 360 0.01% 51.1k -0.52% R4 46628.242075 43.0k + 3.6k -0.06% 46.4k -0.49% Rout = 5k */ void CalcResistors() { const double pA[] = { 3.80240, 2.96940, 2.59700, 3.08700, 3.40060 }; const double pASum = 15.8564; double R[5], Rs, Rout; int k; const double R_E96[] = { 42200, 53600, 60400, 51100, 46400 }; double pA_E96[5]; Rs = 10000; Rout = 1/Rs; for (k = 0; k < 5; k++) { R[k] = Rs / (pA[k] / pASum); Rout += 1/R[k]; printf("R[%d] = %f\n", k, R[k]); } Rout = 1 / Rout; // Output impedance printf("Rout = %f\n", Rout); // Reverse calculation: pA_E96[] from R_E96[]. // Rs = 10000, pASum = 15.8564 printf("pA_E96[] = { "); for (k = 0; k < 5; k++) { pA_E96[k] = Rs * pASum / R[k]; printf("%f", pA_E96[k]); if (k != 4) printf(", "); else printf(" };\n"); } } /*** The main program. ***/ int main() { int len; // Number of samples to generate int fs = 44100; // Sample rate short int* data; // Calculate resistor values for MCU implementation // CalcResistors(); // Test the white noise algorithm. #if 0 lfsr = 1; len = 100 * fs; data = GetNoise(len, GetWhiteValue); WriteAU("test_white.au", fs, data, len); free(data); #endif // Test the original algorithm. Setup0(); len = 100 * fs; data = GetNoise(len, GetPinkValue0); WriteAU("test0.au", fs, data, len); free(data); // Test the modified algorithm. Setup1(); len = 100 * fs; data = GetNoise(len, GetPinkValue1); WriteAU("test1.au", fs, data, len); free(data); return 0; }