pink/brown noise generator

Share icon

Thu, 01/23/2014 - 08:58 am

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