Use 2 frequencies to achieve 8kbps

common.py

@@ -1,3 +1,4 @@
+import numpy as np
 import hashlib
 import struct
 import logging
@@ -12,7 +13,7 @@ Tc = 1.0 / Fc
 Tsym = 1e-3
 Nsym = int(Tsym / Fs)
 
-F0 = Fc - 1e3
+F0 = Fc
 F1 = Fc + 1e3
 Nsym = int(Tsym / Ts)
 baud = int(1/Tsym)
@@ -47,6 +48,21 @@ def to_bytes(bits):
     byte = sum(b << i for i, b in enumerate(bits))
     return chr(byte)
 
+def load(fname):
+    x = np.fromfile(open(fname, 'rb'), dtype='int16')
+    x = x / scaling
+    t = np.arange(len(x)) / Fs
+    return t, x
+
+class Signal(object):
+    def __init__(self, fd):
+        self.fd = fd
+    def send(self, sym, n=1):
+        sym = sym.imag * scaling
+        sym = sym.astype('int16')
+        for i in range(n):
+            sym.tofile(self.fd)
+
 if __name__ == '__main__':
 
     import pylab

recv.py

@@ -1,27 +1,22 @@
 import numpy as np
-from numpy import linalg
 import pylab
 
 import sys
 import struct
 import logging
+import itertools
 logging.basicConfig(level=0, format='%(message)s')
 log = logging.getLogger(__name__)
 
 import sigproc
+import show
 from common import *
 
-NFFT = 256
 COHERENCE_THRESHOLD = 0.9
 
 CARRIER_DURATION = 300
 CARRIER_THRESHOLD = int(0.9 * CARRIER_DURATION)
 
-def load(fname):
-    x = np.fromfile(open(fname, 'rb'), dtype='int16') / scaling
-    t = np.arange(len(x)) / Fs
-    return t, x
-
 def norm(x):
     return np.sqrt(np.dot(x.conj(), x).real)
 
@@ -73,32 +68,59 @@ def find_start(x, start):
     log.info('Carrier starts at {:.3f} ms'.format(start * Tsym * 1e3 / Nsym))
     return start
 
-def equalize(symbols):
+def extract_symbols(x, freq, offset=0):
+    Hc = exp_iwt(-freq, Nsym) / (0.5*Nsym)
+    func = lambda y: np.dot(Hc, y)
+    for _, symbol in iterate(x, Nsym, advance=Nsym, func=func):
+        yield symbol
+
+def demodulate(x, freq, filt):
+    S = extract_symbols(x, freq)
+    S = np.array(list(filt.apply(S)))
+    for bits in sigproc.modulator.decode(S):  # list of bit tuples
+        yield bits
+
+def equalize(x, freqs):
+    prefix = [1]*300 + [0]*100
+    symbols = list(itertools.islice(extract_symbols(x, Fc), len(prefix)))
     bits = np.round(np.abs(symbols))
     bits = np.array(bits, dtype=int)
-    prefix = [1]*300 + [0]*100 + ([1]*10 + [0]*10)*20 + [0]*100
-    n = len(prefix)
-    if all(bits[:n] == prefix):
+    if all(bits[:len(prefix)] != prefix):
+        return None
 
-        S = symbols[:n]
+    log.info( 'Prefix OK')
+    x = x[len(prefix)*Nsym:]
+    filters = {}
+    for freq in freqs:
+        training = ([1]*10 + [0]*10)*20 + [0]*100
+        S = list(itertools.islice(extract_symbols(x, freq), len(training)))
 
-        A = np.array([ S[1:], S[:-1], prefix[:-1] ]).T
-        b = prefix[1:]
+        filt = sigproc.Filter.train(S, training)
+        filters[freq] = filt
 
-        b0, b1, a1 = linalg.lstsq(A, b)[0]
-        y = np.array(list(sigproc.lfilter([b0, b1], [1, -a1], symbols)))
-        constellation(y)
+        y = np.array(list(filt.apply(S))).real
+        #constellation(y)
 
-        prefix_bits = y[:n] > 0.5
-        noise = y[:n] - prefix_bits
-        assert(all(prefix_bits == np.array(prefix)))
-        log.info( 'Prefix OK')
+        train_result = y > 0.5
+        assert(all(train_result == np.array(training)))
+
+        noise = y - train_result
         Pnoise = power(noise)
-        log.debug('Noise sigma={:.4f}, SNR={:.1f} dB'.format( Pnoise**0.5, 10*np.log10(1/Pnoise) ))
+        log.debug('{:10.1f}Hz: Noise sigma={:.4f}, SNR={:.1f} dB'.format( freq, Pnoise**0.5, 10*np.log10(1/Pnoise) ))
+
+        x = x[len(training)*Nsym:]
+
+    results = []
+    for freq in freqs:
+        results.append( demodulate(x, freq, filters[freq]) )
+
+    bitstream = []
+    for block in itertools.izip(*results):
+        for bits in block:
+            bitstream.extend(bits)
+
+    return bitstream
 
-        data_bits = sigproc.qpsk.decode(y[n:])
-        data_bits = list(data_bits)
-        return data_bits
 
 def constellation(y):
     theta = np.linspace(0, 2*np.pi, 1000)
@@ -107,8 +129,8 @@ def constellation(y):
     pylab.subplot(121)
     pylab.plot(y.real, y.imag, '.')
     pylab.plot(np.cos(theta), np.sin(theta), ':')
-    keys = np.array(sigproc.qpsk._enc.values())
-    pylab.plot(keys.real, keys.imag, 'o')
+    points = np.array(sigproc.modulator.points)
+    pylab.plot(points.real, points.imag, 'o')
     pylab.grid('on')
     pylab.axis('equal')
     pylab.axis(np.array([-1, 1, -1, 1]) * 1.1)
@@ -138,16 +160,8 @@ def main(t, x):
 
     start = find_start(x, begin)
     x = x[start:] / amp
-    t = t[start:]
 
-    Hc = exp_iwt(-Fc, Nsym) / (0.5*Nsym)
-    func = lambda y: np.dot(Hc, y)
-    symbols = []
-    for _, coeff in iterate(x, Nsym, advance=Nsym, func=func):
-        symbols.append(coeff)
-    symbols = np.array(symbols)
-
-    data_bits = equalize(symbols)
+    data_bits = equalize(x, [F0, F1])
     if data_bits is None:
         log.info('Cannot demodulate symbols!')
     else:

send.py

@@ -7,19 +7,23 @@ import sys
 import time
 import signal
 import logging
+import itertools
+
 logging.basicConfig(level=0, format='%(message)s')
 log = logging.getLogger(__name__)
 
 import sigproc
 from common import *
 
+dev_null = open('/dev/null')
+
 def play(fd):
-    args = ['aplay', fd.name, '-f', 'S16_LE', '-c', '1', '-r', str(int(Fs))]
+    args = ['aplay', fd.name, '-q', '-f', 'S16_LE', '-c', '1', '-r', str(int(Fs))]
     ret = sp.call(args=args)
     assert ret == 0
 
 def record(fname):
-    args = ['arecord', fname, '-f', 'S16_LE', '-c', '1', '-r', str(int(Fs))]
+    args = ['arecord', fname, '-q', '-f', 'S16_LE', '-c', '1', '-r', str(int(Fs))]
     p = sp.Popen(args=args)
     p.stop = lambda: os.kill(r.pid, signal.SIGINT)
     return p
@@ -30,43 +34,47 @@ log.info('MODEM Fc={}KHz, {} BAUD'.format(Fc/1e3, baud))
 
 class Symbol(object):
     t       = np.arange(0, Nsym) * Ts
-#    c0      = np.sin(2 * np.pi * F0 * t)
-#    c1      = np.sin(2 * np.pi * F1 * t)
-    carrier = np.exp(2j * np.pi * Fc * t)
-
-class Signal(object):
-    def __init__(self, fd):
-        self.fd = fd
-    def send(self, sym, n=1):
-        sym = sym.imag * scaling
-        sym = sym.astype('int16')
-        for i in range(n):
-            sym.tofile(fd)
+    c0      = np.exp(2j * np.pi * F0 * t)
+    c1      = np.exp(2j * np.pi * F1 * t)
 
 sym = Symbol()
 
 data = open('data.send', 'r').read()
-data = pack(data)
 
-bits = list(to_bits(data))
+def start(sig, c):
+    sig.send(c*0, n=100)
+    sig.send(c*1, n=300)
+    sig.send(c*0, n=100)
+
+def train(sig, c):
+    for i in range(20):
+        sig.send(c*1, n=10)
+        sig.send(c*0, n=10)
+    sig.send(c*0, n=100)
 
 if __name__ == '__main__':
 
     with open('tx.int16', 'wb') as fd:
         sig = Signal(fd)
-        sig.send(sym.carrier*0, n=100)
-        sig.send(sym.carrier*1, n=300)
-        sig.send(sym.carrier*0, n=100)
-        for i in range(20):
-            sig.send(sym.carrier*1, n=10)
-            sig.send(sym.carrier*0, n=10)
+        start(sig, sym.c0)
+        train(sig, sym.c0)
+        train(sig, sym.c1)
+        carriers = [sym.c0, sym.c1]
 
-        sig.send(sym.carrier*0, n=100)
-        for s in sigproc.qpsk.encode(bits):
-            sig.send(sym.carrier * s)
+        bits = to_bits(pack(data))
+        symbols_iter = sigproc.modulator.encode(list(bits))
+        symbols_iter = itertools.chain(symbols_iter, itertools.repeat(0))
+        while True:
+            symbols = itertools.islice(symbols_iter, len(carriers))
+            symbols = np.array(list(symbols))
+            sig.send(np.dot(symbols, carriers))
+            if all(symbols == 0):
+                break
 
 
     r = record('rx.int16')
+    start = time.time()
     p = play(fd)
-    time.sleep(0.2)
+    log.debug('Took %.2f seconds', time.time() - start)
+    time.sleep(0.1)
     r.stop()

show.py

@@ -1,4 +1,6 @@
-def spectrogram(t, x):
+import pylab
+
+def spectrogram(t, x, Fs, NFFT=256):
     ax1 = pylab.subplot(211)
     pylab.plot(t, x)
 
@@ -7,3 +9,10 @@ def spectrogram(t, x):
         NFFT=NFFT, Fs=Fs, noverlap=NFFT/2, cmap=pylab.cm.gist_heat)
 
     pylab.show()
+
+if __name__ == '__main__':
+    import sys
+    import common
+    fname, = sys.argv[1:]
+    t, x = common.load(fname)
+    spectrogram(t, x, common.Fs)

sigproc.py

@@ -1,4 +1,5 @@
 import numpy as np
+from numpy import linalg
 
 def lfilter(b, a, x):
     b = np.array(b) / a[0]
@@ -14,11 +15,28 @@ def lfilter(b, a, x):
         y_ = [u] + y_[1:]
         yield u
 
-class QPSK(object):
+
+class Filter(object):
+    def __init__(self, b, a):
+        self.b = b
+        self.a = a
+
+    def apply(self, x):
+        return lfilter(self.b, self.a, x)
+
+    @classmethod
+    def train(cls, S, training):
+        A = np.array([ S[1:], S[:-1], training[:-1] ]).T
+        b = training[1:]
+        b0, b1, a1 = linalg.lstsq(A, b)[0]
+        return cls([b0, b1], [1, -a1])
+
+
+class QAM(object):
     def __init__(self, bits_per_symbol):
         self._enc = {tuple(): 0.0} # End-Of-Stream symbol
         index = 0
-        amps = [0.4, 0.6, 0.8, 1.0]
+        amps = [0.6, 1.0]
         N = (2 ** bits_per_symbol) / len(amps)
         for a in amps:
             for i in range(N):
@@ -27,6 +45,7 @@ class QPSK(object):
                 self._enc[k] = v * a
                 index += 1
         self._dec = {v: k for k, v in self._enc.items()}
+        self.points = self._enc.values()
         self.bits_per_symbol = bits_per_symbol
 
     def encode(self, bits):
@@ -40,15 +59,12 @@ class QPSK(object):
         keys = np.array(self._dec.keys())
         for s in symbols:
             index = np.argmin(np.abs(s - keys))
-            bits = self._dec[ keys[index] ]
-            for b in bits:
-                yield b
+            yield self._dec[ keys[index] ]
 
+modulator = QAM(bits_per_symbol=4)
 
-qpsk = QPSK(bits_per_symbol=6)
-
-def test_qpsk():
-    q = QPSK(bits_per_symbol=2)
+def test():
+    q = QAM(bits_per_symbol=2)
     bits = [1,1, 0,1, 0,0, 1,0]
     S = qpsk.encode(bits)
     assert list(qpsk.decode(list(S))) == bits

test.sh

@@ -3,7 +3,7 @@ set -u
 set -x
 set -e 
 
-dd if=/dev/urandom of=data.send bs=1024 count=1
+dd if=/dev/urandom of=data.send bs=1024 count=8
 python send.py
 python recv.py
-diff data.*
+diff data.* || sha256sum data.*