randomized white noise = white noise?

** Surely you know how to post properly by now ?

The Google Groups way is not acceptable on usenet.

** Google Gropers are always good for their unintentional entertainment value.

......... Phil

As I posted in your other thread (I think it was), you can't. At least not unless you are _very_ clever. The best ideas I've seen here are:

1. To have a big memory with the noise you want pre-stored. At 44100 samples/second and one byte per sample, and considering you can get some pretty huge flash memories pretty cheap these days, you could store 23 seconds per megabyte, or several minutes in a "tiny" 16 megabyte memory. I don't think anyone would notice it repeating in that time. (They tend to want to leave the room after a couple minutes, or less, of loud pink or white noise! ;-) You could even put a few different "types" of noise: white, pink, arbitrarily filtered, some test sinewaves, ... in different segments of the memory and select them with high order address bits.

or

1. Use the Atmel to generate a pseudorandom bitstream, which will be "white" noise over the audio range if the bits come out at several times the highest frequency of the range--say 100k to 200k bits per second. Then filter that with a relatively simple RC filter with poles and zeros spaced (roughly) geometrically along the negative real axis of the s-plane...a pole closest to the origin, then a zero, then a pole, etc. With 1% tolerance caps and 1% or better resistors, I suspect you will get adequately "pink" noise with low enough manufacturing variation for it to be quite useful, though you've never said just what the spec is as far as I've seen. Note: you're likely to need an amplifier at the output end of that. If the white noise has a power density of 1 volt^2/R per 20kHz, then in the 10Hz-20Hz octave, the power is 5e-4 volts^2/R. If the pink is rolling off at 3dB/octave above that, each octave has the same power, and there are about 10 octaves from 20Hz-20kHz. The filter output, then, would be about 5e-3 volts^2/R or 70 millivolts RMS. If you want line-level output, you'll need some gain, and that initial amplitude I assumed is undoubtedly high.

It's possible if you are very clever to build a filter whose coefficients are integer powers of 2, so that multiplication becomes shifting, but even then, you can't do very much with a processor that doesn't do multi-bit shifts in a single cycle...I don't know the details of the Atmel parts but suspect it just does single bit shifts and only on 8 bits at a time.

In even a pretty minimal DSP, I'd say you wouldn't have much trouble doing both the noise generation and the pink filtering, but in an 8-bit part with no hardware multiplication, I think your time will be much better spent looking at other ways to do it. As I recall, I did a pink filter from 10Hz to 22kHz that was accurate to better than 0.05dB over the whole range, in an old Analog Devices ADSP2181...but there wasn't a whole lot of processing time left over.

Cheers, Tom

"Kiviranta, Mikko"

** The term in common use is "red noise".

** That is just no so.

What you "think" is irrelevant.

....... Phil

Thanks for your very interesting posts Tom. I think I have found a "target" so to speak, to focus my efforts on. This algorithm:

Takes white noise samples and outputs pink noise samples. If I can find a way to perform it using the 400 or so instructions I have to work with on the ATTINY13, I have a way to generate pink noise real time using it! Maybe I can use a look up table or somehow do it so I dont need floating point....hmm

I.e. if white is a signed byte, then I can make a lookup table for all possible white * the coefficients. Then I can do the same for b0, b1, and b2. But how will this affect the output? Not too well I imagine.

If the samples themselves are random, it doesn't much matter if you reorder them. But if the samples represent a digitization of noise that's bandlimited below the Nyquist frequency, scrambling the samples will reduce the amplitude and increase the bandwidth, essentially squashing the spectrum. Random samples are uncorrelated, but the samples of bandlimited noise will be correlated, so trashing the correlation changes the nature of the signal.

John

OK, if you implement that as shown it's a filter with three poles, three zeros. The frequency response is roughly "pink" -- if I set the response at 1kHz (assuming 44100 samples/second) to match a perfect filter, then the response error (that is, this filter minus perfect) has a bit less than 1dB ripple at low frequencies, and some loss at high: there's a local min of +0.012dB at 68Hz, a max of +0.915dB at

264Hz, a local min of -0.002dB at 1036Hz, a local max of +0.675dB at 3316Hz, and a min of -3.10dB at 18876Hz.

BUT...consider that first relationship, for b0. If you round 0.99765 to an 8-bit signed number (or even an unsigned one), what do you get? The problem is that you need to represent that number accurately. Since you are dealing with, presumably, a uniformly distributed amplitude in your "white" input, the result may not be too bad, but I'd run a simulation on whatever you propose to do to see just what sort of spectral shape results. You might, for example, set your lookup table so that _on_average_ it represents multiplication by as close as you can to the correct value, but clearly any one individual entry won't be very close to the right value. A good way to think of the error as a percentage is not the .99765 coefficient itself, but rather 1-0.99765

-- the pole it represents on the z plane is at .99765, and its position very strongly affects the frequency response of signals on the unit circle near z=1.

It may be that you can use the fact that you are dealing with large amplitude white noise to save your hind end; if you were implementing a filter that needed to deal with wide dynamic range music, you'd be in deep doo-doo because of roundoff errors.

Cheers, Tom

...

You might try a crude approximation like the following, which uses

16-bit fixed point arithmetic, with binary point at middle. unsigned int b0, b1, b2, tmp, white, tw; do { tw = white/16; tmp = tw*3; b0 = b0 - b0/256 + tmp/2; b1 = b1 - b1/16 + tw/4; b2 = b2/2 + white + tw; tmp += b0 + b1 + b2; outp (tmp/256, PORTB); while (1); In code produced by the avr gcc port, this loop is around 75 instructions and 80 cycles if is quite crude, eg white = ADDP + 16 * white; and around 75 instructions and 80+m+d cycles if is less crude, eg white = (ADDP + MULP * white) % MODP; where m=cycles for rcall __mulhi3, d=cycles for rcall __udivmodhi4, and ADDP, MULP, MODP are appropriate numbers, eg numbers that might look like 2539, 27943, and 65521. [I don't know if the avr gcc port has a random() function; if it does, it's probably better and faster.]

You can calculate the values the above code corresponds to. Eg, b0-b0/256 is 0.99610*b0, vs. your desired 0.99765*b0, while 3*white/16 is white*0.1875 vs. desired white*0.1848, and so forth. Also, the loop is fast enough that it should be ok with more filter stages (b3=..., b4=...). Perhaps that would give you more leeway with the constants.

-jiw

Assuming you meant for b1 to have white/4, not tw/4, that doesn't do too bad, and is a good illustration of being clever with coefficients. The ripples in the error are not quite as well behaved, but the peak-to-peak error is actually less than with the original coefficients. If I match up to zero error at 1kHz, then 20Hz is

-0.907dB, peak at 31Hz of -0.552dB, valley at 111Hz of -1.420dB, peak at 4196Hz of +1.156dB, valley at 18872Hz of -1.837dB.

But we still don't know what's good enough for the application...

Cheers, Tom

Tom Bruhns wrote: [re code, tw = white/16; tmp = tw*3; b0 = b0 - b0/256 + tmp/2; b1 = b1 - b1/16 + tw/4; b2 = b2/2 + white + tw; tmp += b0 + b1 + b2; outp (tmp/256, PORTB); approximating Asa's b0 = 0.99765 * b0 + white * 0.0990460; b1 = 0.96300 * b1 + white * 0.2965164; b2 = 0.57000 * b2 + white * 1.0526913; tmp = b0 + b1 + b2 + white * 0.1848; ]

Right - should have been b1 = b1 - b1/16 + white/4;

Thanks for the analysis! What program, or spreadsheet, do you use to compute response?

Hi James -- and Asa too of course...

Scilab is "freeware" and has a rich set of functions to deal with linear system analysis. I suppose you can do the same in Matlab. It's about the same number of lines of code just typed into the interpreter as in your program to generate the response, generate a matching "true pink" response, and plot the difference. I use a state-space notation to put in the design:

a=[1-1/256,0,0;0,1-1/16,0;0,0,.5]; b=[3/32;1/4;1+1/16] c=[1,1,1]; d=3/16; pinkfilt=syslin('d',a,b,c,d) hz=ss2tf(pinkfilt) // z-domain transfer function [pinkmag,fr]=frmag(hz,22051) // magnitude response // 22051 points linearly spaced means a point every Hz // Now drop 0Hz point which just causes a divide by // zero below if you don't. Also makes the index // of each point be the freq in Hz pinkmag=pinkmag(2:22051); fr=fr(2:22051); const=pinkmag(1000)^2*fr(1000) // Gen. true pink response truepinkmag=sqrt(const ./ fr); pinkdb=20*log10(pinkmag); truepinkdb=20*log10(truepinkmag); errordb=pinkdb-truepinkdb; plot(errordb(20:20000))

You can search "by hand" for min and max values, or get an idea where to look from the graph and use min and max functions on ranges from errordb.

[peak,idx]=max(errordb(2001:8000)) //remember to add idx to 2001 for Hz where the peak occurs

Cheers, Tom

James Waldby wrote:

Wow thats amazing. I am going to have to spend a few days to get my head wrapped around this. Thank you for the very interesting response and juicy code snippet! I will let you know how it goes when I get it implemented.

To answer the question re: what the app needs: this pink noise will be listened to by a person, and is for relaxation purposes (to be general) so I would imagine its not too critical if there is 1db ripple here or there. I could be wrong. But its not for an experiment or something like that which requires a very accurate slope.

Asa

If you want to do relaxation, then find an abandoned stretch of beach and record a couple of hours of waves breaking at the shore. :-)

Good Luck! Rich

Hm. To get started and to get a better understanding of whats happening I tried to implement the original floating point un-efficient filter in open watcom running on my laptop, to generate 10 seconds of noise hence the 441000:

float b0, b1, b2;

signed int white;

for(sampnum = 0 ; sampnum < 441000 ; sampnum++) { white = (rand() * 2) - (RAND_MAX); b0 = 0.99765 * b0 + white * 0.0990460; b1 = 0.96300 * b1 + white * 0.2965164; b2 = 0.57000 * b2 + white * 1.0526913; noise[sampnum] = b0 + b1 + b2 + white * 0.1848; }

But the output array noise[], according the debugger, gets filled with numbers greater than and less than 32767/-32767. Is this because the noise is literally trying to be gaussian and has peaks and valleys which can go to infinity? The large numbers are quite frequent, maybe 3 or 4 out of ten is beyond a 16bit integer. Should I normalize it somehow? Or am I doing something wrong here.

The input is random between -32767 and 32767, (white)

The filter does its thing.

The output array and the b0, b1, b2 coeff's are float. So it should be able to handle everything.

I tried scaling it by dividing the final number by 10. And now the spectra of the output appears to be a straight 3db slope! now I am going to try it with more coefficients and see how good it gets. noise[sampnum] = (b0 + b1 + b2 + white * 0.1848) / 10;

Hi Tom,

Lets say I need to modify the filter slightly so that from 20 to 9000 hz its got a -3db pink rolloff, but above that its has a -6db roll off? Is there a program I can use to generate the coefficients? Will it look like the last filter, i.e. the b0, b1, b2, but with maybe different numbers? Or does the whole structure change?

I'll try to answer two posts at once here...

Note that when a rolloff slope changes from 3dB/octave to 6dB/octave, it's a gradual thing. You can draw lines on a Bode plot that are straight and intersect at a sharp corner, but in practice, the transition is gradual. So the difference in response between just the

3dB/octave and the added 6dB/octave in the top octave will be relatively small. In fact, if you get Scilab and run the simple code I suggested, you can see a plot of the response--easy to plot the response of your filter, or the difference between it and "true pink." And what you will see with your original coefficients is that they cause the high end to drop off faster than "pink" anyway, by just about what you are suggesting with that "6dB slope in the top octave". So you may have that already. In fact, in Scilab (or Matlab) it's easy to make most of what I posted be a routine that will plot the response or whatever, and you just change the coefficients and run the routine to see the effect. Even with a slow processor, it should do the calcs in a second or less, especially if you used fewer points logarithmically spaced, but even if you do that one point per Hz hack I posted, it's still fast.

OK, have a look at the relationships in your filter to see why you get large values pretty often. Notice that b2, for example, remembers over half its previous value, and adds MORE than one times the input to itself. So let's say you have three white samples in a row that are in the range 0.5 to 1.0 (assuming that white runs from -1 to +1). Note that:

b2(k+1) = 0.57*b2(k) + 1.05*white(k) = 0.57* ( 0.57*b2(k-1) + 1.05*white(k-1)) + 1.05*white(k) =

0.57*(0.57*(0.57*b2(k-2)+1.05*white(k-2))+1.05*white(k-1))+1.05*white(k) = 0.185*b(k-2)+0.34*white(k-2)+0.60*white(k-1)+1.05*white(k)

If the average value of white = 0 (no DC bias), then on average b(k-2) will be zero. A string of three "whites" in a row at 0.5 or above will then, on average, result in an output from just THAT equation that's greater than 1.

The way around this is to scale "white" down, and don't let it cover the whole range from -1 to +1. How much do you scale it? The absolute worst case would be if white stayed at +1 (or -1) for a long time. For your filter (which technically is three one-pole resonators), you can easily get this value: first you find the value for each of the "b" values, then add them plus 0.1848 for the current "white" value that's in the output. How do you get the value for each "b"? Simple: it will reach steady state when the contribution from the "old" b decays by exactly the amount of the new input. So for b0, white*0.099 = (1-0.99765)*b0(steady-state). So b0(steady-state)=42.12*white. That is, if the input is locked at 1, b0 will eventually get to 42.12. So an absolutely safe scaling would be to not let "white" get above

1/42.12. Now in practice, the probability of "white" staying that close to full scale in one direction for long enough to do that is so tiny (and in fact would be zero for any practical pseudorandom number generator) that you can safely make the scaling larger than that. For b2, the safe scaling (by a similar calculation) would be 0.408.

Note that since you are dealing with NOISE here, if you can gracefully handle the overflow, letting it happen once every few thousand samples is likely not going to even be noticable. I suspect that if you scale the white input by 0.25, you will be fine. Try it and see how that works. Your filter probably does a very decent job of turning the uniform amplitude distribution into a Gaussian amplitude distribution, so the mere fact that you see about 3 or 4 out of 10 that's larger than the range of the WHITE input tells me that the standard deviation is quite close to full scale. Then the scaling I suggested, 0.25, should get you to about one in 10,000 that overflows. Note that you do not have to SHIFT your white; you can just duplicate the top bit into the next two bits (for a signed value), and it will still be white. The bits in a uniform-amplitude-distribution random number should be all uncorrelated, and it doesn't matter which bits you drop.

Cheers, Tom

Also, note that if you scale the whole set of values (0.099046,

0.296516, 1.052691, 0.1848) by the same amount, it's equivalent to scaling WHITE by that amount. So you may be able to make some of those coefficients easier to deal with at the same time you do the scaling. If you scale the whole set the same, the filter frequency response will stay the same and only the amplitude will change. Don't change the coefficients applied to b0, b1, b2 or it will change the frequency response. Beware that you will NEED to keep b0, especially, as a 16 bit number at least. 8 bits will almost certainly NOT be enough for it, especially as you scale the input smaller. (Welcome to the world of digital filters with poles at a tiny fraction of the sample frequency--poles close to z=1.)

Cheers, Tom

Thanks for the fascinating reply as usual Tom! The code snippet Waldby posted is intended to receive random numbers in the range -32767 to

+32767 right? And his "unsigned int" on the atmel is a 16bit number correct? So when he says 16bits fixed point with the decimal point in the middle, hes talking about an 8bit whole number with 8 fractional bits? Seems kind of small to me. So that makes me think he meant 32bits split whole/fractional. Does this sound right to you? Or James? :)

I meant 16 bits, with randoms in the range 0 to 65535. (More precisely, up to MODP, which I gave as 65521.) I doubt that you will be able to hear any difference between 16 and 32 bit results, and suggest that you go ahead and try it both ways. For the 32 bit case, substitute "unsigned long int" for each occurrence of "unsigned int", and change "outp (tmp/256, PORTB);" to instead output either the most significant byte or the most significant 2 bytes, depending on your DAC arrangement.

The number of instructions almost doubles (from 76 to 141) and the total cycles more than doubles (from about 77 to about 160, due to 4 brne and 16 ldd 2-cycle instructions) but still is fast enough, at 62500 loops per second (given a 10MHz clock), to produce your desired 44100 outputs per second. As before, the instruction counts are based on a gcc-avr-port compilation, with a crude random number generator. If you use a statistically better RNG you probably can get better measurements on a spectrum analyzer, but perhaps no one can hear the difference.

-jiw