Pink noise generator
TonyWaite
Posts: 219
Hi,
I need to knock out a 'high-quality' pink noise generator for the acoustic engineers next door.
There are perfectly good circuits on the Web, based on provoking a diode and then amplifying it;
but I would prefer a Prop because it enables a glorious interface!
My first thought would be to use a WAV player via the SD card; playing a pink-noise sample downloaded
from the Web; later analysing it using one of the FFT/Fast Hartley functions and displaying the result on the QVGA
of Rayman's PSM. The Prop would also control the subsequent amplification stages of high-power audio output to
give pre-selected levels of sound.
Then I found pink-noise generation code for other processors and wondered whether anyone had thought about
this or done anything similar on the Propeller. Amongst other things, it would be useful to be able to alter the crest
factor.
Regards,
T o n y
I need to knock out a 'high-quality' pink noise generator for the acoustic engineers next door.
There are perfectly good circuits on the Web, based on provoking a diode and then amplifying it;
but I would prefer a Prop because it enables a glorious interface!
My first thought would be to use a WAV player via the SD card; playing a pink-noise sample downloaded
from the Web; later analysing it using one of the FFT/Fast Hartley functions and displaying the result on the QVGA
of Rayman's PSM. The Prop would also control the subsequent amplification stages of high-power audio output to
give pre-selected levels of sound.
Then I found pink-noise generation code for other processors and wondered whether anyone had thought about
this or done anything similar on the Propeller. Amongst other things, it would be useful to be able to alter the crest
factor.
Regards,
T o n y
Comments
Sample the output of the diode with the prop and then you could do
real time processing of the random stream of values before you send it
to the outputs for use.
This way you can make a pretty video display using the prop and lay claim
to having a true random noise source as well as a modifiable one.
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
24 bit LCD Breakout Board now in. $24.99 has backlight driver and touch sensitive decoder.
@mctrivia
Thank you very much for your advice; upon which I am pondering deeply.
I'm still thinking about how to use the 'genuinely random' facility to generate the pink noise; I might buy
Leon lunch next week and pick his brains!
To n y
Leon
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Amateur radio callsign: G1HSM
No sooner had I thought of lunch ...
Pseudorandom probably fails this test and has subharmonic peaks at 1/2 1/3 1/4 1/5 the sampling rate, that's a guess.
I think I'm right on the basis that this phenomenon seems to produce a sound with unexpected privacy uses.
It might even fool an FFT into seeing random noise. <<<lip zip>>>
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
VIRAND, If you spent as much time SPINNING as you do Trolling the Forums,
you'd have tons of awesome code to post! (Note to self)
http://forums.parallaxinc.com/ccrma.stanford.edu/~jos/sasp/Example_Synthesis_1_F_Noise.html
I converted it to a Scilab program (Scilab has a Matlab conversion utility):
clear; clf(); // fs=22050; T=2; N=round(fs*T); if modulo(N,2)==1, N=N+1; end rand('seed',sum(getdate())); s=rand(1,N,'normal'); // a=fs/10; F=zeros(1,N/2); for i=1:N/2 F(i)=1/(1+(a*i/N)^2)^0.25; end F=[noparse][[/noparse]F mtlb_fliplr(F)]; f=linspace(0,1/T,N/2+1); subplot(2,2,1); plot2d1('gll',f,F(1:N/2+1),color('blue')); // S=fft(s); S=S.*F; subplot(2,2,2); plot2d1('gll',f,abs(S(1:N/2+1)),color('blue')); // s=real(ifft(S)); s=s/max(abs(s)); subplot(2,2,3); plot(s); sound(s,fs); // subplot(2,2,4); histplot([noparse][[/noparse]-2:0.2:2],s);It ran OK, and played the noise through the PC speakers. It should be quite easy to implement the algorithm on a dsPIC or ARM chip with a DAC. I'm not sure if it is feasible on a Propeller, though.
Leon
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Amateur radio callsign: G1HSM
Post Edited (Leon) : 1/14/2010 2:36:54 PM GMT
Phil Pilgrim has just uploaded a PASM pink-noise generator, that I have been playing
with for the last few minutes.
I'm just looking up the Fast Hartley thread too ...
Regards,
T o n y
The day is looking up: it started without any pink-noise algorithm and suddenly there's two!
I downloaded the Scilab suite - and ran your code straight off.
I'd better go and do a bit more metalwork now and investigate tonight ...
T o n y
ultimately to drive an I2S DAC, currently experimental.
The three approaches I identified from researching pink noise are:
1) inverse FFT of a suitable 1/f randomized spectrum, compute intensive, tricky for continuous stream.
2) filtering Gaussian white noise with a suitable digital filter (expensive on the Prop due to the
number of multiplies involved. Basically one pole per 2 octaves combine to a 3dB/octave average slope,
with humps.
3) summing bandlimited gaussian noise sources, one per octave. Since noise sums by power this gives 3dB/octave.
The last one, if done naively, gives humps in the spectrum, per octave. The method was proposed by Voss
and McCartney. Put simply you maintain a vector of gaussian samples, one per octave, and update each
at a frequency fs/(2^i). So one is always new every sample, one changes every other, one changes every
fourth, etc. The sum of the current values is the output sample.
Several people(*) improved on the basic method by adding randomness to the timing of the updates to the samples
in the vector, maintaining the same average frequencies of update. This tends to even out the humps in
the noise spectrum.
(*) notably Tramell (http://www.ridgerat-tech.us/tech/newpink.htm)
and Allen Downey (https://www.dsprelated.com/showarticle/908.php)
I took the idea of Allen Downey, using the trailing zero count in a random integer to select a single
generator to update each step (means workload is constant per sample). The noise spectrum from this has a bit of a gap at middle frequencies
and peaks too high at the Nyquist frequency.
Playing around I found that weighting only the most frequently updated gaussian generator by a factor
around 1.4 or 1.5 improved the gap (but peaks the high frequencies more).
The spectrum is pretty good for the lower octaves, it just needs some treatment in the higher frequencies
to flatten it. Time for a carefully selected digital filter (using only add/sub/shifts) to fixup the
higher end (and provide cutoff before the Nyquist freqency).
1 + 15/16 z H(z) = ------------- 1 + 3/4 zWhich is a handful of PASM instructions to compute, but makes a much nicer spectrum:
Overall I use the Tyche RNG which is a fast modern design that passes BigCrush, but takes only 12 PASM instructions
to generate each 32 bit output value
Neves and Araujo: https://pdfs.semanticscholar.org/bfca/81c7fbc6b2c32430dc756b936a6b0c2d0585.pdf
Computing the gaussians from uniform random numbers can be done in several ways, I chose to pre-compute
a table of 1024 gaussian samples, and select 3 entries at random from the table to sum.
The distribution of 1024 gaussians:
The distribution of 500,000 values formed from summing 3 entries from a table of 1024:
This gives a universe of 2^30 possible gaussians to select at random with a fairly tolerable amount
of code - generating 3 table indexes, reading each from hub ram, summing.
Given 16 bits of randomness to determine the choice of stochastic gaussian, and two gaussian samples (one added every
sample, one added to the chosen stochastic gaussian, we get a single pink noise sample in about 4us.
The 16 stochastic guassians aren't directly summed for each sample, the sum is cached in a single variable
and this is updated every time one of them changes
Note that the plots above are all done using Python (numpy, scipy.signal, matplotlib.pyplot)
If there's interest I can post the PASM generator code (when its less in a state of flux maybe).
At some point when its all hooked up to a DAC real spectrum analyzer screen shots should follow.
frequency fs ]
A green 3dB/octave guide line is plotted.
Amazingly, when it comes to pink noise generation, by the last few weeks, since Chip has started the:
https://forums.parallax.com/discussion/169298/adc-sampling-breakthrough/p1
I've been pursuing almost the same objectives, but, in my case, l was looking for some way to improve the P2 ADC Spurious-free dynamic range (SFDR), stimulated by the many readings I was doing about that subject, untill I found this:
https://analog.com/en/analog-dialogue/articles/adc-input-noise.html
Then, following the track of pink noise I've hit:
firstpr.com.au/dsp/pink-noise/#Update
Thru it, I've found Allen Downey's, '2016 blog:
https://dsprelated.com/showarticle/908.php
, and the unique comment it received, at the bottom of the page, wroten by "Stenz".
Keeping the tracking process, I've deeply dived at "Stenz" Github:
https://github.com/Stenzel/newshadeofpink
Reading its terms of use (that sure didn't seem to be restritive, by any means) I've found his web site, with a page that surelly demostrates the efficacy of the whole pink noise generation process, that can be extended, usign the same maths, till fSample/2, if needed:
stenzel.waldorfmusic.de/post/pink/
Since them, a whole new concept catched my attention; and it's nice and clean, as for my eyes and for my ears too.
If you didn't yet had a chance to look at it, I suggest you to spent some minutes...
I hope it'll enlightenning to you as it was to me.
Henrique
I'm happy you've received it this well.
Henrique
{ ┌──────────────────────────────────────────────────────────┐ │ Pink Noise Generator. │ │(c) Copyright 2009 Philip C. Pilgrim (propeller@phipi.com)│ │ See end of file for terms of use. │ └──────────────────────────────────────────────────────────┘ This program generates pink noise using the Voss-McCartney algorithm. Version History ─────────────── 2010.01.13: Initial alpha release } CON _clkmode = xtal1 + pll16x _xinfreq = 5_000_000 NRAND = 30 'Number of random numbers to sum. == 1: white noise; > 1: getting "pinker". AOUT_PIN = 26 'Output pin number. SHIFT = >| (NRAND- 1) 'Amount to shift volume to keep from overflowing. PUB start | p p := clkfreq / 44000 cognew(@pink_noise, @p) DAT org 0 pink_noise rdlong period,par 'Get the period. mov ctra,ctra0 'Set up ctra for DUTY mode. mov dira,dira0 'Set DUTY mode pin as output. movd :seed,#rand 'Point :seed to rand array. mov count,#NRAND 'Need to seed NRAND rands. :seed mov 0-0,count 'Seed rand with count. add :seed,inc_dst 'Increment address to rand. djnz count,#:seed 'Next seed. mov time,cnt 'Initialize wait time. add time,period 'Now + period. mov sum,#NRAND*(NRAND+1)/2 'Sum is sum of all rands. mov epoch,#0 'Initialize epoch. mainlp mov count,#1 'Need to count trailing consecutive ones in epoch shl count,#NRAND-1 ' to get address of next rand. mov acc,#rand 'Initialize address to head of array. jmp #:next_bit 'Start counting... :count_ones test epoch,count wc 'Is bit a one? if_c add acc,#1 ' Yes: Increment rand array address. if_nc mov acc,#rand ' No: Start over. :next_bit shr count,#1 wz 'Done counting? if_nz jmp #:count_ones ' No: Back for another bit. movs :ldrand0,acc 'Set up loads and stores with rand pointer. movs :ldrand1,acc movd :strand,acc :ldrand0 mov acc,0-0 'Get old value. and acc,rand_mask 'Just the part that was added. sub sum,acc 'Subtract it. :ldrand1 mov acc,0-0 'Get old value again. mov count,#32-SHIFT 'LFSR is 32-SHIFT cycles. :randlp test acc,taps wc 'Do the LFSR tango... rcr acc,#1 djnz count,#:randlp '...32-SHIFT times :strand mov 0-0,acc 'Save the new value back to rand array. and acc,rand_mask 'Mask to get the low bits. add sum,acc 'Add the new value into the sum. waitcnt time,period 'Wait for the period to be over. mov frqa,sum 'Write the sum into the DUTY value. add epoch,#1 'Increment the epoch number. jmp #mainlp 'And back to do it all over again. ctra0 long %00110 << 26 | AOUT_PIN 'DUTY mode counter, outputting on AOUT_PIN. dira0 long 1 << AOUT_PIN 'Mask to set AOUT_PIN to output. taps long %11010000_00000001 'Taps for LFSR. rand_mask long $ffff_ffff >> SHIFT 'Use 32-SHIFT lsbs of random number. inc_dst long 1 << 9 'Added to instruction to increment dest address. rand res NRAND 'Array of random numbers. acc res 1 'General purpose accumulator. count res 1 'General purpose counter. sum res 1 'Amplitude that gets written to frqa. epoch res 1 'Number of times through loop. time res 1 'Time of next write to frqa. period res 1 'Period of main loop in clock ticks. {{ ┌──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐ │ TERMS OF USE: MIT License │ ├──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤ │Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation │ │files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, │ │modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software│ │is furnished to do so, subject to the following conditions: │ │ │ │The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.│ │ │ │THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE │ │WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR │ │COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, │ │ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. │ └──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘ }}it works now and its got some fiendishly clever techniques.
o 12 tap FIR filter implemented with 2 64-entry table lookups (since its on a bit stream)
o linearly interpolating many waveforms simultaneously, exploiting the power-of-two ratios between waveforms
o directly stealing a bit stream history from inside a suitably tapped LFSR
In short a definite candidate for PASM implementation
I do worry about its very parsimonious use of only one random bit per sample, and my next investigations will be looking at the autocorrelation behaviour to see if its great looking spectrum hides shortcomings in the randomness department.
Talking of its spectrum:
I'll be following all your findings with deep interest.
Sixteen years ago, I've decided to left all about C++, and everything else closely related, to my son and his generation of coders.
Assembly is ever fine, crystal clear as a timing diagram.
To my eyes, the rest seems like Spaghetti with Bolognese sauce; I can have it for a meal, but I don't want to follow its footprints inside my digestive tract.
Phil Pilgrim brought his own '2010 pasm version of a pink noise generator too.
At the end, we all should learn a little bit more in the direction of attaining better results.
It'll a win, win situation for everyone.
per sample. Core of the code below, also attached my test app which uses a board I made with
I2S DAC and ADC clocking at 98MHz (48kHz * 2048 to be precise)
stenzel ' get bit mask from counter's lowest 1 and bit reversing add counter, #1 and counter, mask mov bitmask, counter neg t, counter and bitmask, t ' find lowest set bit in counter rev bitmask, #32-POWER ' reverse so highest bit changes most rapidly shl lfsr, #1 wc ' progress the LFSR, leaving carry with latest bit if_c xor lfsr, H46000001 andn dec, bitmask ' clever stuff to update one of the bits in dec and inc mov t, inc ' the inc value of the bit is copied to the dec register, then and t, bitmask ' depending on the carry flag (a random bit), flip the inc's bit at that position or dec, t if_c xor inc, bitmask add accum, inc ' these two instructions do linear interpolation on 16 waveforms simultaneuously! sub accum, dec mov resul, accum mov t, lfsr ' lookup lower bits 5:0 - doing a 12 tap FIR filter on last 12 bits output shl t, #2 and t, #$FC add t, table_addr rdlong t, t add resul, t mov t, lfsr ' lookup next bits 11:6 shr t, #4 and t, #$FC add t, table_addr add t, #$100 rdlong t, t add resul, t stenzel_ret retHow do you intend to use the resultings you've got so far?
Sorry by my spaghetti-alike writing...
I was meanting to ask how do you intend to inject the pink noise at your setup:
- forwarding pink noise binary data to the DAC input, then summing the resulting analog pink noise to any other audio sample and connect it to the ADC input?