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

Sum of two shifted LUT output probability mass function

LUT output

class GaussLUT(Elaboratable):
    def __init__(self, i, o):
        self.i = i # 6 bits uniform random
        self.o = o # signed, 3 bits range is -3 to 3
    def elaborate(self, platform):
        m = Module()
        with m.Switch(self.i):
            for i in range(64):
                with m.Case(i):
                    sign = -1 if (i&1) else 1
                    if i < 20:
                        m.d.sync += self.o.eq(0)
                    elif i < 50:
                        m.d.sync += self.o.eq(1*sign)
                    elif i < 62:
                        m.d.sync += self.o.eq(2*sign)
                    else:
                        m.d.sync += self.o.eq(3*sign)
        return m

class GaussGen(Elaboratable):
    def __init__(self, prbs, enable, nbits=18):
    self.i = prbs.get(96 + 4*(nbits-7))
        self.o = Signal(signed(nbits))
    def elaborate(self, platform):
        m = Module()

        gluto = Array(Signal(signed(3), name='gluto_{}'.format(a)) for a in range(16))
        for i in range(16):
            m.submodules["glut{}".format(i)] = GaussLUT(i=self.i[6*i:6*i+6], o=gluto[i], clk=self.clk)

        add1 = Array(Signal(signed(4), name='add1_{}'.format(a)) for a in range(8))
        for i in range(8):
            m.d.sync += add1[i].eq(gluto[2*i] + gluto[2*i+1])

        add2 = Array(Signal(signed(5), name='add2_{}'.format(a)) for a in range(4))
        for i in range(4):
            m.d.sync += add2[i].eq(add1[2*i] + add1[2*i+1])

        add3 = Array(Signal(signed(6), name='add3_{}'.format(a)) for a in range(2))
        for i in range(2):
            m.d.sync += add3[i].eq(add2[2*i] + add2[2*i+1])

        add4 = Signal(signed(7)) # +- 48 range
        m.d.sync += add4.eq(add3[0] + add3[1])

        # two triangular PDF signals:
        tpdf = Array(Signal(signed(len(self.o)-6), name='tpdf_{}'.format(a)) for a in range(2))
        i_tri = self.i[96:]
        bits_tripart = len(self.o) - 7
        for i in range(2):
            i_tri2 = i_tri[2*bits_tripart*i:2*bits_tripart*i+2*bits_tripart]
            m.d.sync += tpdf[i].eq(i_tri2[:bits_tripart] - i_tri2[bits_tripart:])

        m.d.sync += self.o.eq((add4 << bits_tripart) + tpdf[0] + tpdf[1])

        return m
Histogram from simulation

The orange in the middle is the triangular PDF created from two uniform distributions.

Distribution

Ideal and actual PDF

PDF

PDF