pink/brown noise generator
Okapi
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
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
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
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:
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
I took the liberty of using SI units, which allows Igor to be smarter about what units are displayed in axes.
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