equalizer: use PRBS for equalization sequence

amodem/dsp.py

@@ -109,3 +109,14 @@ class MODEM(object):
                 if error_handler:
                     error_handler(received=received, decoded=decoded)
                 yield bits
+
+
+def prbs(reg, poly, bits):
+    ''' Simple pseudo-random number generator. '''
+    mask = (1 << bits) - 1
+    while True:
+        yield reg & mask
+        reg = reg << 1
+        size = poly.bit_length() - 1
+        if reg >> size:
+            reg = reg ^ poly

amodem/equalizer.py

@@ -5,29 +5,23 @@ import numpy as np
 from numpy.linalg import lstsq
 
 import itertools
-import random
 
 
 class Equalizer(object):
 
-    _constellation = [1, 1j, -1, -1j]
-
     def __init__(self, config):
         self.carriers = config.carriers
         self.omegas = 2 * np.pi * np.array(config.frequencies) / config.Fs
         self.Nfreq = config.Nfreq
         self.Nsym = config.Nsym
 
-    def train_symbols(self, length, seed=0, constant_prefix=16):
-        r = random.Random(seed)
-
-        def random_symbol():
-            ''' use low-level randomness for cross-version compatibility. '''
-            return self._constellation[r.getrandbits(2)]
+    def train_symbols(self, length, constant_prefix=16):
+        r = dsp.prbs(reg=1, poly=0x1100b, bits=2)
+        constellation = [1, 1j, -1, -1j]
 
         symbols = []
         for _ in range(length):
-            symbols.append([random_symbol() for _ in range(self.Nfreq)])
+            symbols.append([constellation[next(r)] for _ in range(self.Nfreq)])
 
         symbols = np.array(symbols)
         # Constant symbols (for analog debugging)

tests/test_dsp.py

@@ -72,3 +72,22 @@ def test_overflow():
     for i in range(10000):
         s = 10*(r.normal() + 1j * r.normal())
         quantize(q, s)
+
+
+def test_prbs():
+    r = list(itertools.islice(dsp.prbs(reg=1, poly=0x7, bits=2), 4))
+    assert r == [1, 2, 3, 1]
+
+    r = list(itertools.islice(dsp.prbs(reg=1, poly=0x7, bits=1), 4))
+    assert r == [1, 0, 1, 1]
+
+    r = list(itertools.islice(dsp.prbs(reg=1, poly=0xd, bits=3), 8))
+    assert r == [1, 2, 4, 5, 7, 3, 6, 1]
+
+    r = list(itertools.islice(dsp.prbs(reg=1, poly=0xd, bits=2), 8))
+    assert r == [1, 2, 0, 1, 3, 3, 2, 1]
+
+    period = 2 ** 16 - 1
+    r = list(itertools.islice(dsp.prbs(reg=1, poly=0x1100b, bits=16), period))
+    r.sort()
+    assert r == list(range(1, 2 ** 16))