Effect of sample rate? Also Java port

[johnmgithub]

Hi Stefan,

Just came across this, nice work! I do my programming in Java, so I made a Java version to investigate the algorithm. Code is below if of interest and contributed freely for any use. I made a slight change from your standard version, replacing the first zero entry in the pnmasks table by 1 << 7 to slightly reduce the low frequency extension and flatten off the spectrum at 5 Hz or so rather than 2.5 Hz, the aim being to keep generated levels a little more consistent across different interface low end rolloffs. I also adapted it to return single values rather than writing blocks of 16 samples.

Should there be an adaptation for different sample rates, different filter coefficient tables, perhaps?

/**
 * Pink noise generation using the algorithm by Stefan Stenzel https://github.com/Stenzel/newshadeofpink
 */
class PinkNoise {
    private static final float PINK_BIAS = 3.0f;

    // lookup for bitreversed masks
    // Replaced 1 << 7 by 0 to reduce low frequency extension slightly, 
    // spectrum levels off at approx 5 Hz rather than 2.5 Hz
    private static final int[] pnmask = {
           0,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,0x0800,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,
        0x08 << 7,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,0x0800,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,
        0x04 << 7,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,0x0800,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,
        0x08 << 7,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,0x0800,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,
           2 << 7,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,0x0800,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,
        0x08 << 7,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,0x0800,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,
        0x04 << 7,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,0x0800,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,
        0x08 << 7,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,0x0800,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000, 
           0,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,0x0800,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,
        0x08 << 7,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,0x0800,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,
        0x04 << 7,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,0x0800,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,
        0x08 << 7,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,0x0800,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,
           2 << 7,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,0x0800,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,
        0x08 << 7,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,0x0800,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,
        0x04 << 7,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,0x0800,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,
        0x08 << 7,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000,0x0800,0x4000,0x2000,0x4000,0x1000,0x4000,0x2000,0x4000};
    private static final int[] otherMasks = {
        0, 0x40000,0x20000,0x40000,0x10000,0x40000,0x20000,0x40000,0x08000,0x40000,0x20000,0x40000,0x10000,0x40000,0x20000,0x40000};
    // Lookup tables for 12-tap FIR filter with these coefficients:
    private static final float[] FIR = {
        1.190566f, 0.162580f, 0.002208f, 0.025475f, -0.001522f, 0.007322f, 
        0.001774f, 0.004529f, -0.001561f, 0.000776f, -0.000486f, 0.002017f};
    // 1st precalculated FIR lookup table, also for bias correction
    private static final float[] pfira = new float[64];
    // 2nd precalculated FIR lookup table
    private static final float[] pfirb = new float[64]; 
    private static int instance_cnt; // used for decorrelation in case of multiple instances

    private int lfsr; // linear feedback shift register
    private int inc;  // increment for all noise sources (bits)
    private int dec;  // decrement for all noise sources
    private int accui; // int interpretation of accu
    private int pnCount; // index to pnmask[]
    private int otherCount; // index to otherMasks

    private static float F(float cf, int m, int shift) {
        return 0.0625f * cf * (2 * ((m) >>> shift & 1) - 1);
    }

    private static float FA(int n) {
        float result = F(FIR[0], n, 0);
        for (int i = 1; i < 6; i++) {
            result += F(FIR[i], n, i);
        }
        return result - PINK_BIAS;
    }

    private static float FB(int n) {
        float result = F(FIR[6], n, 0);
        for (int i = 1; i < 6; i++) {
            result += F(FIR[6+i], n, i);
        }
        return result;
    }

    static {
        for (int n = 0; n < 64; n++) {
            pfira[n] = FA(n);
            pfirb[n] = FB(n);
        }
    }

    private float pink(int bitmask){
        int bit = lfsr >> 31;            // spill random to all bits
        dec &= ~bitmask;                 // blank old decrement bit
        lfsr <<= 1;                      // shift lfsr
        dec |= inc & bitmask;            // copy increment to decrement bit
        inc ^= bit & bitmask;            // new random bit
        float sample = Float.intBitsToFloat(accui); // get value as float  
        accui  += inc - dec;             // integrate   
        lfsr ^= bit & 0x46000001;        // update lfsr
        return sample                    // save output and ...
             + pfira[lfsr & 0x3F]        // add 1st half FIR & subtract bias
             + pfirb[lfsr>>6 & 0x3F];    // add 2nd half precalculated FIR
    }

    float nextValue(){
        int mask;
        if (otherCount == 0){
            mask = pnmask[pnCount];
            pnCount = (pnCount + 1) & 0x0ff;
        }else{
            mask = otherMasks[otherCount];
        }
        otherCount = (otherCount + 1) & 0x0f;
        return pink(mask);
    }

    PinkNoise(){
        lfsr  = 0x5EED41F5 + instance_cnt++;        // seed for lfsr, decorrelate multiple instances    
        accui = Float.floatToRawIntBits(PINK_BIAS);
        pnCount = 0;                                // counter from zero    
        inc   = 0x04CCCC;                           // balance initial states to avoid DC 
        dec   = 0x04CCCC;        
    }

}

[Stenzel]

Hello John,

Thanks for using the code! I am not sure if the change to the pnmask is a good idea, if I remember correctly this could lead to long term instability because of DC buildup from the CIC integrator. There is also no need to change coefficients, spectrum remains steady at 1/f independent of sample rate.

Best, Stefan

[johnmgithub]

Hmm, didn't notice any issues but I haven't run the generator for more than 10 minutes or so at a time. Would it be better to have two zero entries in pnmask instead of the two 1 << 7 entries?

[Stenzel]

I'd recommend zeros, that should be sufficient to disable noise at a specific rate.

[johnmgithub]

OK, thanks very much. I've updated the code section above as well.