pink/brown noise generator

Share icon

Hello -

I wrote a function to generate a pink (P = 1 / f) or brown (P = 1/ f ^ 0.5) noise wave (see below). The function is slow because I am using a loop to create many individual sin waves (columns in a single wave) that later get summed together into a single wave (using a matrixOP). I find myself stuck using the loop because I can't think of a different way to give each component sin wave it's own (random) phase and power.

I am seeking advice on how I could possibly use Matrix op or some other way to eliminate the loop from the function.

Note: If you are testing this function, you will find that the BuildWave wave gets big! and can lead to insufficient memory errors. One way to avoid this is to keep the Duration variable small (ex. 200 ms). I am planning to tackle this problem later by having the build wave sample rate to be twice the frequency of the LowPassCutOff. For now though, I am creating the wave a a sample interval that my hardware "likes".

Threadsafe Function WB_PinkAndBrownNoise(Amplitude, Duration, LowPassCutOff, HighPassCutOff, FrequencyIncrement, PinkOrBrown)
		variable Amplitude, Duration, LowPassCutOff, HighPassCutOff, frequencyIncrement, PinkOrBrown
		variable phase = (abs(enoise(2)) * Pi)
		variable NumberOfBuildWaves = floor((LowPassCutOff - HighPassCutOff) / FrequencyIncrement)
		make /free /n = (Duration / 0.005, NumberOfBuildWaves) BuildWave
		SetScale /P x 0,0.005,"ms", BuildWave
		variable Frequency = HighPassCutOff
		variable i = 0
		variable localAmplitude
		do
			phase = ((abs(enoise(2))) * Pi) // random phase generator
			if(PinkOrBrown == 0)
				localAmplitude = 1 / Frequency
			else
				localAmplitude = 1 / (Frequency ^ .5)
			endif
			
			MultiThread BuildWave[][i] = localAmplitude * sin(2 * Pi * (Frequency*1000) * (5 / 1000000000) * p + phase) // 
			Frequency += FrequencyIncrement
			i += 1
		while (i < NumberOfBuildWaves)
		
		MatrixOp /o /NTHR = 0   SegmentWave = sumRows(BuildWave)
		
		SetScale /P x 0, 0.005,"ms", SegmentWave
		SegmentWave /= NumberOfBuildWaves
		
		Wavestats /q SegmentWave
		variable scalefactor = Amplitude/(V_max - V_min)
		SegmentWave *= ScaleFactor
End
In all cases where speed is important you want to minimize the number of operations in the loop. In this case I'd start by factoring out

2 * Pi * (1000) * (5 / 1000000000)

as this does not change during the implied loop.

The second point I'd make is that there is no reason for you to generate the phase and test for pink/brown inside the loop. Do it once outside and store the phase and frequency in a single waves (together with the argument for the sine) so that you can replace the multithread statement with MatrixOP. Note that enoise() is thread-safe so you should be able to fill that wave with a multithread instruction.

HTH,

A.G.
WaveMetrics, Inc.

January 23, 2014 at 05:45 pm - Permalink
AG,

Thank you for your help! Adam also sent me a modified function code. Adam eliminated the loop entirely and the function got significantly slower!? I've found that moving the phase or frequency calculation out of the loop also slows down the function. Simplifying the "(2 * Pi * (Frequency*1000) * (5 / 1000000000) * p + phase)" line did seem to increase the speed significantly with each variable I removed, but a second round of testing didn't realize as large improvements. I was also surprised that the "IF" statement doesn't slow down the function.

Here is what I have settled on for now:
Threadsafe Function WB_PinkAndBrownNoise3(Amplitude, Duration, LowPassCutOff, HighPassCutOff, FrequencyIncrement, PinkOrBrown)// Pink = 0, Brown = 1
		variable Amplitude, Duration, LowPassCutOff, HighPassCutOff, frequencyIncrement, PinkOrBrown
		Variable startTime = StopMsTimer(-2)
		variable phase 
		variable NumberOfBuildWaves = floor((LowPassCutOff - HighPassCutOff) / FrequencyIncrement)
		make /free /n = (Duration / 0.005, NumberOfBuildWaves) BuildWave
		SetScale /P x 0,0.005,"ms", BuildWave
		variable Frequency = HighPassCutOff
		variable i = 0
		variable localAmplitude
		//make /free /n = (NumberOfBuildWaves) PhaseWave
		//Multithread PhaseWave[] = ((abs(enoise(2))) * Pi)
		do
			phase = ((abs(enoise(2))) * Pi) // random phase generator
			if(PinkOrBrown == 0)
				localAmplitude = 1 / Frequency
			else
				localAmplitude = 1 / (Frequency ^ .5)
			endif
			
			MultiThread BuildWave[][i] = localAmplitude * sin( Pi * Frequency * 1e-05 * p + phase) //  factoring out Pi * 1e-05 actually makes it a tiny bit slower
			//MultiThread BuildWave[][i] = localAmplitude * sin(2 * Pi * (Frequency*1000) * (5 / 1000000000) * p + phase) // Multithread of sin funciton is the BOMB!!
			Frequency += FrequencyIncrement
			i += 1
		while (i < NumberOfBuildWaves)
		
		MatrixOp /o /NTHR = 0   SegmentWave = sumRows(BuildWave)
		
		SetScale /P x 0, 0.005,"ms", SegmentWave
		SegmentWave /= NumberOfBuildWaves
		
		Wavestats /q SegmentWave
		variable scalefactor = Amplitude / (V_max - V_min)
		SegmentWave *= ScaleFactor
		printf "WB_PinkAndBrownNoise3() took %g ms.\r", (StopMsTimer(-2) - startTime)/ 1000

End

January 23, 2014 at 09:43 pm - Permalink
Is there any reason you are not doing this sum as an (inverse) FFT? I'm not sure the commented-out version of the switch statement does exactly what you want for the frequency cutoffs and increment, but presuming those were only there to limit NumberofBuildWaves and noting that having zeros in the FFT input doesn't speed it up, you can just get rid of them. The code as is is equivalent to using zero LowPassCutOff and FrequencyIncrement, and infinite HighPassCutoff. 10 minutes worth of data sampled at 200 kHz is about as much as fits in memory on my laptop, and takes about 40 seconds to generate. 200 ms worth is generated in less than 10 ms. How fast does it need to be?

I took the liberty of using SI units, which allows Igor to be smarter about what units are displayed in axes.


 Function WB_PinkAndBrownNoise4(Amplitude, Duration, LowPassCutOff, HighPassCutOff, FrequencyIncrement, PinkOrBrown)// Pink = 0, Brown = 1
	variable Amplitude, Duration, LowPassCutOff, HighPassCutOff, frequencyIncrement, PinkOrBrown
	
	Variable startTime = StopMsTimer(-2)

	Variable samplerate = 1e3/0.005 // in Hz
	Variable samples = Duration * samplerate * 1e-3 // if duration is in ms
	samples = 2 * ceil(samples/2) // even number of points for FFT
	Make/O/C/N=(samples/2+1) magphase
	SetScale/P x 0,samplerate/samples,"Hz" magphase

	switch(PinkOrBrown)
		case 0:		// pink noise
			magphase[1,] = cmplx(1/x, enoise(Pi))
			break
		case 1:		// brown noise
			magphase[1,] = cmplx(1/sqrt(x), enoise(Pi))
			break
		default:	// brown noise
			magphase[1,] = cmplx(1/sqrt(x), enoise(Pi))
	endswitch
//	Variable lowsample = min(max(1,pnt2x(magphase, LowPassCutoff*1e3)),samples/2) // if LowPassCutOff is in kHz
//	Variable highsample = min(max(1,pnt2x(magphase, HighPassCutoff*1e3)),samples/2) // if HighPassCutOff is in kHz
//	Variable skipsample = min(1,pnt2x(magphase, frequencyIncrement*1e3)) // if frequencyIncrement is in kHz
//	
//	switch(PinkOrBrown)
//		case 0:		// pink noise
//			magphase[lowsample,highsample;skipsample] = cmplx(1/x, enoise(Pi))
//			break
//		case 1:		// brown noise
//			magphase[lowsample,highsample;skipsample] = cmplx(1/sqrt(x), enoise(Pi))
//			break
//		default:	// brown noise
//			magphase[lowsample,highsample;skipsample] = cmplx(1/sqrt(x), enoise(Pi))
//	endswitch
	MatrixOp/O SegmentWave = IFFT(p2Rect(magphase),3)
	MatrixOp/FREE ScaleFactor = Amplitude/(maxVal(SegmentWave)-minVal(SegmentWave)))
	MatrixOp/O SegmentWave = SegmentWave * ScaleFactor[0] // ScaleFactor is a 1x1 matrix

	SetScale /P x 0, 1/samplerate,"s", SegmentWave

	printf "WB_PinkAndBrownNoise4() took %g ms.\r", (StopMsTimer(-2) - startTime)/ 1000
End

January 24, 2014 at 08:16 pm - Permalink