Generating Pink Noise (Flicker, 1/f) in an FPGA

Intro

The Harmon Instruments signal generator provides simulated phase noise modulation. Typical signal sources include a flicker noise component at low offset frequencies. The 3 dB/octave slope is more complex to create than 6 dB/octave which can be produced with a simple integrator.

Stochastic Voss-McCartney

The Voss-Mcartney method of generating pink noise is a multirate algorithm. The output is the sum of many random variables updated at sample rates that are power of two multiples of each other. There are a few other useful references. In the stochastic version, rather than updating each random variable in the sum at fixed intervals, they are updated at randomized intervals averaging the update rates in the non stochastic version.

In this implementation, a sum of 32 values is used, numbered 0 to 31. Value 0 has a probability of update of 0.5 each clock cycle, value 1 0.25, value n 2^-(n+1), etc. It might be desirable to add a value that updates every cycle for improved high frequency conformance, but it's not required in this application.

Here's a simple Python model:

import numpy as np
import random
ram = np.zeros(32)

def get_pink():
    for i in range(len(ram)):
        if random.getrandbits(1):
            ram[i] = random.gauss(0, 1)
            break
    return np.sum(ram)

The low frequency deviation in the plot below is due to the number of samples used, not the generator.

Pink Noise plot

The plot below is the output of the phase noise source set to pure pink as measured with a spectrum analyzer. This shows good performance over at least 6 decades of frequency range.

Pink Noise plot

nMigen Implementation

I'm unaware of a closed form solution to the output spectral density. Numerical evaluation gives aproximately 48100 / sqrt(Hz) at 1 Hz assuming the gaussian input has a standard deviation of 10102.5.

The value in RAM at index 31 is updated twice as often as it should be in this code. That may be fixed at some point in the future. At 250 MSPS, that results in noise that should be below 0.058 Hz being at 0.058 Hz.

On Artix-7, usage is 58 LUTs, 96 registers. An external gaussian white noise source is required as well as 31 pseudo random bits per clock.

class PinkNoise(Elaboratable):
    def __init__(self, i, urnd):
        self.i = i # 20 bit signed white gaussian noise
        self.urnd = urnd # 31 pseudo random bits
        self.o = Signal(signed(25))
    def elaborate(self, platform):
        m = Module()
        # count trailing zeros
        bits_1 = self.urnd
        cond_2 = bits_1[:16] == 0
        bits_2 = Signal(15)
        result_2 = Signal()
        cond_3 = bits_2[:8] == 0
        bits_3 = Signal(7)
        result_3 = Signal(2)
        cond_4 = bits_3[:4] == 0
        bits_4 = Signal(3)
        result_4 = Signal(3)
        ptr = Signal(5)
        m.d.sync += [
            result_2.eq(cond_2),
            bits_2.eq(Mux(cond_2, bits_1[16:31], bits_1[:15])),
            result_3.eq(Cat(cond_3, result_2)),
            bits_3.eq(Mux(cond_3, bits_2[8:15], bits_2[:7])),
            result_4.eq(Cat(cond_4, result_3)),
            bits_4.eq(Mux(cond_4, bits_3[4:7], bits_3[:3])),
            ptr.eq(
                Mux(bits_4[0],
                    Cat(C(0,2),result_4),
                    Mux(bits_4[1],
                        Cat(C(1,2),result_4),
                        Mux(bits_4[2],
                            Cat(C(2,2),result_4),
                            Cat(C(3,2),result_4)
                        )
                    )
                )
            ),
        ]

        ram = Memory(width=len(self.o) - len(ptr), depth=2**len(ptr))
        wrport = m.submodules.wrport = ram.write_port(domain='sync')
        rdport = m.submodules.rdport = ram.read_port(domain='comb')
        m.d.comb += [
            wrport.en.eq(1),
            wrport.addr.eq(ptr),
            wrport.data.eq(self.i),
            rdport.addr.eq(ptr),
        ]
        i_pipe = Signal(signed(len(self.i)))
        ram_pipe = Signal(signed(len(rdport.data)))
        m.d.sync += [
            i_pipe.eq(self.i),
            ram_pipe.eq(rdport.data),
            self.o.eq(self.o + i_pipe - ram_pipe),
        ]

        return m