From 807c03a8e8585b13bd3e5466072cd6ab6998b345 Mon Sep 17 00:00:00 2001
From: Roman Zeyde <roman.zeyde@gmail.com>
Date: Fri, 13 Feb 2015 14:47:28 +0200
Subject: [PATCH] equalizer: use PRBS for equalization sequence

---
 amodem/dsp.py       | 11 +++++++++++
 amodem/equalizer.py | 14 ++++----------
 tests/test_dsp.py   | 19 +++++++++++++++++++
 3 files changed, 34 insertions(+), 10 deletions(-)

diff --git a/amodem/dsp.py b/amodem/dsp.py
index d8da7d5..ba7f88f 100644
--- a/amodem/dsp.py
+++ b/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
diff --git a/amodem/equalizer.py b/amodem/equalizer.py
index 7afa108..de43363 100644
--- a/amodem/equalizer.py
+++ b/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)
diff --git a/tests/test_dsp.py b/tests/test_dsp.py
index ea45285..9a16e58 100644
--- a/tests/test_dsp.py
+++ b/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))