noise.py

from scipy.fftpack import fft
from matplotlib import pyplot

N = 65536
f_s = 48000
f_c = 1000

prng_lfsr = 0x32B71700
prng_tap = 0xb874

def prng_lfsr_next():
  global prng_lfsr
  lsb = prng_lfsr & 1
  prng_lfsr >>= 1

  if (lsb):
      prng_lfsr ^= prng_tap

  return prng_lfsr


prng_xorshift_state = 0x32B71700

def prng_xorshift_next():
  global prng_xorshift_state
  x = prng_xorshift_state
  x ^= (x << 13) & 0xffffffff
  x ^= (x >> 17) & 0xffffffff
  x ^= (x <<  5) & 0xffffffff
  prng_xorshift_state = x
  return x

def sample(fn):
  return ((fn() & 0xffff) - 0x7fff) / 0x8000

def filter(input, f_s, f_c, integer = False):
  RC = 1.0 / (2 * 3.1415926 * f_c)
  dt = 1.0 / f_s
  a = dt / (RC + dt)
  if integer:
    # There's a missing 65536/f_s term, but it's so small it's not worthwhile
    # ... and also should probably round-to-nearest on the // operation
    a = int(2 * 3.1415926 * 65536) * f_c // f_s
  x = 0
  y = []
  for x_n in input:
    x += (x_n - x) * a
    y.append(x)
  return y

def freq(input):
  f = fft(input)
  N = len(f)
  scale = N ** -0.5
  return [abs(x) * scale for x in f[:N // 2]]

xordata = [sample(prng_xorshift_next) for _ in range(N)]
lfsrdata = [sample(prng_lfsr_next) for _ in range(N)]
bins = [i * f_s / N for i in range(N //2)]

pyplot.loglog(bins[10:], freq(xordata)[10:])
pyplot.loglog(bins[10:], freq(lfsrdata)[10:])
#pyplot.loglog(bins[10:], freq(filter(xordata, f_s, f_c, False))[10:])
#pyplot.loglog(bins[10:], freq(filter(xordata, f_s, f_c, True))[10:])
pyplot.show()