refactor configuration in sigproc.MODEM object

calib.py

@@ -1,12 +1,13 @@
 #!/usr/bin/env python
 import numpy as np
 import common
+import config
 import sigproc
 import wave
 
 Tsample = 1
-t = np.arange(int(Tsample * common.Fs)) * common.Ts
-sig = np.exp(2j * np.pi * common.Fc * t)
+t = np.arange(int(Tsample * config.Fs)) * config.Ts
+sig = np.exp(2j * np.pi * config.Fc * t)
 sig_dump = common.dumps(sig)
 
 

colorhash.py

@@ -4,7 +4,12 @@ import binascii
 
 lines = sys.stdin.read().strip().split('\n')
 for line in lines:
-    head, tail = line.split(' ', 1)
+    try:
+        head, tail = line.split(' ', 1)
+    except ValueError:
+        print line
+        continue
+
     bars = ''
     try:
         data = binascii.unhexlify(head)

common.py

@@ -5,22 +5,6 @@ import numpy as np
 import logging
 log = logging.getLogger(__name__)
 
-Fs = 32e3
-Ts = 1.0 / Fs
-
-frequencies = (1 + np.arange(9)) * 1e3
-carrier_index = 0
-Fc = frequencies[carrier_index]
-Tc = 1.0 / Fc
-
-symbols = np.array([complex(x, y)
-                   for x in np.linspace(-1, 1, 8)
-                   for y in np.linspace(-1, 1, 8)]) / np.sqrt(2)
-
-Tsym = 1e-3
-Nsym = int(Tsym / Ts)
-baud = int(1/Tsym)
-
 scaling = 32000.0  # out of 2**15
 SATURATION_THRESHOLD = 1.0
 
@@ -52,18 +36,14 @@ def check_saturation(x):
         raise SaturationError(peak)
 
 
-def load(fileobj, time=False):
-    return loads(fileobj.read(), time=time)
+def load(fileobj):
+    return loads(fileobj.read())
 
 
-def loads(data, time=False):
+def loads(data):
     x = np.fromstring(data, dtype='int16')
     x = x / scaling
-    if time:
-        t = np.arange(len(x)) / Fs
-        return t, x
-    else:
-        return x
+    return x
 
 
 def dumps(sym, n=1):
@@ -129,14 +109,3 @@ def icapture(iterable, result):
 
 def take(iterable, n):
     return np.array(list(itertools.islice(iterable, n)))
-
-if __name__ == '__main__':
-
-    import pylab
-    t = pylab.linspace(0, Tsym, 1e3)
-    x = pylab.sin(2 * pylab.pi * Fc * t)
-    pylab.plot(t / Tsym, x)
-    t = pylab.linspace(0, Tsym, Nsym + 1)
-    x = pylab.sin(2 * pylab.pi * Fc * t)
-    pylab.plot(t / Tsym, x, '.k')
-    pylab.show()

config.py

@@ -0,0 +1,17 @@
+import numpy as np
+
+Fs = 32e3
+Ts = 1.0 / Fs
+
+frequencies = (1 + np.arange(9)) * 1e3
+carrier_index = 0
+Fc = frequencies[carrier_index]
+Tc = 1.0 / Fc
+
+symbols = np.array([complex(x, y)
+                   for x in np.linspace(-1, 1, 8)
+                   for y in np.linspace(-1, 1, 8)]) / np.sqrt(2)
+
+Tsym = 1e-3
+Nsym = int(Tsym / Ts)
+baud = int(1/Tsym)

recv.py

@@ -2,6 +2,7 @@
 import numpy as np
 import logging
 import itertools
+import functools
 import collections
 import time
 import sys
@@ -13,13 +14,17 @@ import stream
 import sigproc
 import loop
 import train
-from common import *
+import common
+import config
+
+modem = sigproc.MODEM(config)
+
 
 if os.environ.get('PYLAB') == '1':
     import pylab
     import show
-    WIDTH = np.floor(np.sqrt(len(frequencies)))
-    HEIGHT = np.ceil(len(frequencies) / float(WIDTH))
+    WIDTH = np.floor(np.sqrt(len(modem.freqs)))
+    HEIGHT = np.ceil(len(modem.freqs) / float(WIDTH))
 else:
     pylab = None
 
@@ -32,30 +37,30 @@ SEARCH_WINDOW = 10  # symbols
 
 def report_carrier(bufs, begin):
     x = np.concatenate(tuple(bufs)[-CARRIER_THRESHOLD:-1])
-    Hc = sigproc.exp_iwt(-Fc, len(x))
+    Hc = sigproc.exp_iwt(-config.Fc, len(x))
     Zc = np.dot(Hc, x) / (0.5*len(x))
     amp = abs(Zc)
     log.info('Carrier detected at ~%.1f ms @ %.1f kHz:'
              ' coherence=%.3f%%, amplitude=%.3f',
-             begin * Tsym * 1e3 / Nsym, Fc / 1e3,
-             np.abs(sigproc.coherence(x, Fc)) * 100, amp)
+             begin * config.Tsym * 1e3 / config.Nsym, config.Fc / 1e3,
+             np.abs(sigproc.coherence(x, config.Fc)) * 100, amp)
 
 
 def detect(samples, freq):
 
     counter = 0
-    bufs = collections.deque([], maxlen=baud)  # 1 second of symbols
-    for offset, buf in iterate(samples, Nsym):
+    bufs = collections.deque([], maxlen=config.baud)  # 1 second of symbols
+    for offset, buf in common.iterate(samples, config.Nsym):
         bufs.append(buf)
 
-        coeff = sigproc.coherence(buf, Fc)
+        coeff = sigproc.coherence(buf, config.Fc)
         if abs(coeff) > COHERENCE_THRESHOLD:
             counter += 1
         else:
             counter = 0
 
         if counter == CARRIER_THRESHOLD:
-            length = (CARRIER_THRESHOLD - 1) * Nsym
+            length = (CARRIER_THRESHOLD - 1) * config.Nsym
             begin = offset - length
             report_carrier(bufs, begin=begin)
             break
@@ -65,28 +70,30 @@ def detect(samples, freq):
     log.debug('Buffered %d ms of audio', len(bufs))
 
     to_append = SEARCH_WINDOW + (CARRIER_DURATION - CARRIER_THRESHOLD)
-    for _, buf in itertools.islice(iterate(samples, Nsym), to_append):
+    bufs_iterator = common.iterate(samples, config.Nsym)
+    for _, buf in itertools.islice(bufs_iterator, to_append):
         bufs.append(buf)
 
     bufs = tuple(bufs)[-CARRIER_DURATION-2*SEARCH_WINDOW:]
     buf = np.concatenate(bufs)
 
-    offset = find_start(buf, length=Nsym*CARRIER_DURATION)
-    start = begin - Nsym * SEARCH_WINDOW + offset
-    log.info('Carrier starts at {:.3f} ms'.format(start * Tsym * 1e3 / Nsym))
+    offset = find_start(buf, length=config.Nsym*CARRIER_DURATION)
+    start = begin - config.Nsym * SEARCH_WINDOW + offset
+    log.info('Carrier starts at %.3f ms',
+             start * config.Tsym * 1e3 / config.Nsym)
 
     return itertools.chain(buf[offset:], samples)
 
 
 def find_start(buf, length):
-    Hc = sigproc.exp_iwt(Fc, len(buf))
+    Hc = sigproc.exp_iwt(config.Fc, len(buf))
     P = np.abs(Hc.conj() * buf) ** 2
     cumsumP = P.cumsum()
     return np.argmax(cumsumP[length:] - cumsumP[:-length])
 
 
 def receive_prefix(symbols):
-    S = take(symbols, len(train.prefix))[:, carrier_index]
+    S = common.take(symbols, len(train.prefix))[:, config.carrier_index]
     sliced = np.round(S)
 
     nonzeros = np.array(train.prefix, dtype=bool)
@@ -99,10 +106,10 @@ def receive_prefix(symbols):
 
     pilot_tone = S[nonzeros]
 
-    freq_err, mean_phase = sigproc.drift(pilot_tone) / (Tsym * Fc)
+    freq_err, mean_phase = sigproc.drift(pilot_tone) / (config.Tsym * config.Fc)
     expected_phase, = set(np.angle(sliced[nonzeros]) / (2 * np.pi))
 
-    sampling_err = (mean_phase - expected_phase) * Nsym
+    sampling_err = (mean_phase - expected_phase) * config.Nsym
     log.info('Frequency error: %.2f ppm', freq_err * 1e6)
     log.info('Sampling error: %.2f samples', sampling_err)
     return freq_err, sampling_err
@@ -115,7 +122,7 @@ def train_receiver(symbols, freqs):
     if pylab:
         pylab.figure()
 
-    symbols = take(symbols, len(training) * len(freqs))
+    symbols = common.take(symbols, len(training) * len(freqs))
     for i, freq in enumerate(freqs):
         size = len(training)
         offset = i * size
@@ -154,17 +161,17 @@ def demodulate(symbols, filters, freqs, sampler):
     def error_handler(received, decoded, freq):
         errors.setdefault(freq, []).append(received / decoded)
 
-    generators = split(symbols, n=len(freqs))
+    generators = common.split(symbols, n=len(freqs))
     for freq, S in zip(freqs, generators):
         S = filters[freq](S)
 
         if pylab:
             equalized = []
-            S = icapture(S, result=equalized)
+            S = common.icapture(S, result=equalized)
             symbol_list.append(equalized)
 
         freq_handler = functools.partial(error_handler, freq=freq)
-        bits = sigproc.modulator.decode(S, freq_handler)  # list of bit tuples
+        bits = modem.qam.decode(S, freq_handler)  # list of bit tuples
         streams.append(bits)  # stream per frequency
 
     stats['symbol_list'] = symbol_list
@@ -177,15 +184,16 @@ def demodulate(symbols, filters, freqs, sampler):
             stats['rx_bits'] = stats['rx_bits'] + len(bits)
             yield bits
 
-        if i and i % baud == 0:
+        if i and i % config.baud == 0:
             mean_err = np.array([e for v in errors.values() for e in v])
             correction = np.mean(np.angle(mean_err)) / (2*np.pi)
+            duration = time.time() - stats['rx_start']
             log.debug('%10.1f kB, realtime: %6.2f%%, sampling error: %+.3f%%',
                       stats['rx_bits'] / 8e3,
-                      (time.time() - stats['rx_start']) * 100.0 / (i*Tsym),
+                      duration * 100.0 / (i*config.Tsym),
                       correction * 1e2)
             errors.clear()
-            sampler.freq -= 0.01 * correction / Fc
+            sampler.freq -= 0.01 * correction / config.Fc
             sampler.offset -= correction
 
 
@@ -220,19 +228,19 @@ def decode(bits_iterator):
 
 def main(args):
 
-    log.info('Running MODEM @ {:.1f} kbps'.format(sigproc.modem_bps / 1e3))
+    log.info('Running MODEM @ {:.1f} kbps'.format(modem.modem_bps / 1e3))
     start = time.time()
 
     fd = sys.stdin
     signal = stream.iread(fd)
-    skipped = take(signal, args.skip)
-    log.debug('Skipping first %.3f seconds', len(skipped) / float(baud))
+    skipped = common.take(signal, args.skip)
+    log.debug('Skipping first %.3f seconds', len(skipped) / float(modem.baud))
 
-    stream.check = check_saturation
+    stream.check = common.check_saturation
 
     size = 0
-    signal = detect(signal, Fc)
-    bits = receive(signal, frequencies)
+    signal = detect(signal, config.Fc)
+    bits = receive(signal, modem.freqs)
     try:
         for chunk in decode(bits):
             sys.stdout.write(chunk)
@@ -241,7 +249,7 @@ def main(args):
         log.exception('Decoding failed')
 
     duration = time.time() - stats['rx_start']
-    audio_time = stats['rx_bits'] / float(sigproc.modem_bps)
+    audio_time = stats['rx_bits'] / float(modem.modem_bps)
     log.info('Demodulated %.3f kB @ %.3f seconds (%.1f%% realtime)',
              stats['rx_bits'] / 8e3, duration, 100 * duration / audio_time)
 
@@ -252,7 +260,7 @@ def main(args):
     if pylab:
         pylab.figure()
         symbol_list = np.array(stats['symbol_list'])
-        for i, freq in enumerate(frequencies):
+        for i, freq in enumerate(modem.freqs):
             pylab.subplot(HEIGHT, WIDTH, i+1)
             show.constellation(symbol_list[i], '$F_c = {} Hz$'.format(freq))
 

send.py

@@ -7,16 +7,19 @@ import time
 
 log = logging.getLogger(__name__)
 
-import sigproc
 import train
 import wave
 
-from common import *
+import common
+import config
+import sigproc
+
+modem = sigproc.MODEM(config)
 
 
 class Symbol(object):
-    t = np.arange(0, Nsym) * Ts
-    carrier = [np.exp(2j * np.pi * F * t) for F in frequencies]
+    t = np.arange(0, config.Nsym) * config.Ts
+    carrier = [np.exp(2j * np.pi * F * t) for F in modem.freqs]
 
 sym = Symbol()
 
@@ -26,7 +29,7 @@ class Writer(object):
         self.last = time.time()
 
     def write(self, fd, sym, n=1):
-        fd.write(dumps(sym, n))
+        fd.write(common.dumps(sym, n))
         if time.time() > self.last + 1:
             log.debug('%10.3f seconds of data audio',
                       fd.tell() / wave.bytes_per_second)
@@ -46,7 +49,7 @@ def training(fd, c):
 
 
 def modulate(fd, bits):
-    symbols_iter = sigproc.modulator.encode(bits)
+    symbols_iter = modem.qam.encode(bits)
     symbols_iter = itertools.chain(symbols_iter, itertools.repeat(0))
     carriers = np.array(sym.carrier) / len(sym.carrier)
     while True:
@@ -77,14 +80,14 @@ class Reader(object):
 
 def main(args):
     import ecc
-    log.info('Running MODEM @ {:.1f} kbps'.format(sigproc.modem_bps / 1e3))
+    log.info('Running MODEM @ {:.1f} kbps'.format(modem.modem_bps / 1e3))
 
     fd = sys.stdout
 
     # padding audio with silence
-    write(fd, np.zeros(int(Fs * args.silence_start)))
+    write(fd, np.zeros(int(config.Fs * args.silence_start)))
 
-    start(fd, sym.carrier[carrier_index])
+    start(fd, sym.carrier[config.carrier_index])
     for c in sym.carrier:
         training(fd, c)
     training_size = fd.tell()
@@ -94,14 +97,14 @@ def main(args):
     reader = Reader(sys.stdin, 64 << 10)
     data = itertools.chain.from_iterable(reader)
     encoded = itertools.chain.from_iterable(ecc.encode(data))
-    modulate(fd, bits=to_bits(encoded))
+    modulate(fd, bits=common.to_bits(encoded))
 
     data_size = fd.tell() - training_size
     log.info('%.3f seconds of data audio, for %.3f kB of data',
              data_size / wave.bytes_per_second, reader.total / 1e3)
 
     # padding audio with silence
-    write(fd, np.zeros(int(Fs * args.silence_stop)))
+    write(fd, np.zeros(int(config.Fs * args.silence_stop)))
 
 if __name__ == '__main__':
     logging.basicConfig(level=logging.DEBUG,

show.py

@@ -27,11 +27,13 @@ def spectrogram(t, x, Fs, NFFT=256):
 if __name__ == '__main__':
     import sys
     import common
+    from config import Fs, Ts
 
     for fname in sys.argv[1:]:
-        t, x = common.load(open(fname, 'rb'), time=True)
+        x = common.load(open(fname, 'rb'))
+        t = np.arange(len(x)) * Ts
         pylab.figure()
         pylab.title(fname)
-        spectrogram(t, x, common.Fs)
+        spectrogram(t, x, Fs)
 
     pylab.show()

sigproc.py

@@ -2,6 +2,7 @@ import numpy as np
 from numpy import linalg
 
 import common
+from config import Ts, Nsym
 
 
 class Filter(object):
@@ -56,10 +57,14 @@ class QAM(object):
                 error_handler(received=s, decoded=S)
             yield self._dec[S]
 
-modulator = QAM(common.symbols)
 
-bits_per_baud = modulator.bits_per_symbol * len(common.frequencies)
-modem_bps = common.baud * bits_per_baud
+class MODEM(object):
+    def __init__(self, config):
+        self.qam = QAM(config.symbols)
+        self.baud = config.baud
+        self.freqs = config.frequencies
+        self.bits_per_baud = self.qam.bits_per_symbol * len(self.freqs)
+        self.modem_bps = self.baud * self.bits_per_baud
 
 
 def clip(x, lims):
@@ -71,7 +76,7 @@ def power(x):
 
 
 def exp_iwt(freq, n):
-    iwt = 2j * np.pi * freq * np.arange(n) * common.Ts
+    iwt = 2j * np.pi * freq * np.arange(n) * Ts
     return np.exp(iwt)
 
 
@@ -82,14 +87,18 @@ def norm(x):
 def coherence(x, freq):
     n = len(x)
     Hc = exp_iwt(-freq, n) / np.sqrt(0.5*n)
-    return np.dot(Hc, x) / norm(x)
+    norm_x = norm(x)
+    if norm_x:
+        return np.dot(Hc, x) / norm_x
+    else:
+        return 0.0
 
 
 def extract_symbols(x, freq, offset=0):
-    Hc = exp_iwt(-freq, common.Nsym) / (0.5*common.Nsym)
+    Hc = exp_iwt(-freq, Nsym) / (0.5*Nsym)
     func = lambda y: np.dot(Hc, y)
 
-    iterator = common.iterate(x, common.Nsym, func=func)
+    iterator = common.iterate(x, Nsym, func=func)
     for _, symbol in iterator:
         yield symbol
 

test.sh

@@ -4,7 +4,7 @@ set -e
 
 run() {
 	echo "SRC $HOST: $CMD" 1>&2
-	ssh $HOST $CMD
+	ssh $HOST "$CMD"
 }
 
 run_src() {
@@ -20,7 +20,7 @@ run_dst() {
 }
 
 ## generate 1Mbit of random data
-#run_src dd if=/dev/urandom of=data.send bs=125kB count=1 status=none
+run_src dd if=/dev/urandom of=data.send bs=125kB count=1 status=none
 SRC_HASH=`run_src sha256sum data.send | ./colorhash.py`
 
 # modulate data into audio file

wave.py

@@ -6,8 +6,9 @@ import logging
 
 log = logging.getLogger(__name__)
 
-from common import Fs
-Fs = int(Fs)  # sampling rate
+import config
+Fs = int(config.Fs)  # sampling rate
+
 bits_per_sample = 16
 bytes_per_sample = bits_per_sample / 8.0
 bytes_per_second = bytes_per_sample * Fs