Use 4 carriers wit QAM16 to achieve 16kbps.

common.py

@@ -7,32 +7,26 @@ log = logging.getLogger(__name__)
 Fs = 32e3
 Ts = 1.0 / Fs
 
-Fc = 9e3
+frequencies = (np.arange(4) + 8) * 1e3
+carrier_index = len(frequencies)/2
+Fc = frequencies[carrier_index]
 Tc = 1.0 / Fc
 
 Tsym = 1e-3
-Nsym = int(Tsym / Fs)
-
-F0 = Fc
-F1 = Fc + 1e3
 Nsym = int(Tsym / Ts)
 baud = int(1/Tsym)
 
-scaling = 10000.0
+scaling = 32000.0  # out of 2**15
+SATURATION_THRESHOLD = 1.0
 
 LENGTH_FORMAT = '<I'
 
 def pack(data):
     log.info('Sending {} bytes'.format(len(data)))
-    data = struct.pack(LENGTH_FORMAT, len(data)) + data
     return data
 
 def unpack(data):
-    length_size = struct.calcsize(LENGTH_FORMAT)
+    log.info('Received {} bytes'.format(len(data)))
-    length, data = data[:length_size], data[length_size:]
-    length, = struct.unpack(LENGTH_FORMAT, length)
-    data = data[:length]
-    log.info('Decoded {} bytes, leftover: {}'.format(len(data), len(data)-length))
     return data
 
 def to_bits(chars):

errors.py

@@ -0,0 +1,17 @@
+import common
+import sys
+
+tx, rx = sys.argv[1:]
+tx = open(tx).read()
+rx = open(rx).read()
+
+L = min(len(tx), len(rx))
+rx = list(common.to_bits(rx[:L]))
+tx = list(common.to_bits(tx[:L]))
+indices = [index for index, (r, t) in enumerate(zip(rx, tx)) if r != t]
+
+if indices:
+    total = L*8
+    errors = len(indices)
+    print('{}/{} bit error rate: {:.3f}%'.format(errors, total, (100.0 * errors) / total))
+    sys.exit(1)

recv.py

@@ -77,6 +77,7 @@ def extract_symbols(x, freq, offset=0):
 def demodulate(x, freq, filt):
     S = extract_symbols(x, freq)
     S = np.array(list(filt.apply(S)))
+    #constellation(S)
     for bits in sigproc.modulator.decode(S):  # list of bit tuples
         yield bits
 
@@ -98,11 +99,13 @@ def equalize(x, freqs):
         filt = sigproc.Filter.train(S, training)
         filters[freq] = filt
 
-        y = np.array(list(filt.apply(S))).real
+        S = list(filt.apply(S))
-        #constellation(y)
+        y = np.array(S).real
 
         train_result = y > 0.5
-        assert(all(train_result == np.array(training)))
+        if not all(train_result == np.array(training)):
+            pylab.plot(y, '-', training, '-')
+            return None
 
         noise = y - train_result
         Pnoise = power(noise)
@@ -112,7 +115,7 @@ def equalize(x, freqs):
 
     results = []
     for freq in freqs:
-        results.append( demodulate(x, freq, filters[freq]) )
+        results.append( demodulate(x * len(freqs), freq, filters[freq]) )
 
     bitstream = []
     for block in itertools.izip(*results):
@@ -159,15 +162,18 @@ def main(t, x):
     ))
 
     start = find_start(x, begin)
-    x = x[start:] / amp
+    x = x[start:]
+    peak = np.max(np.abs(x))
+    if peak > SATURATION_THRESHOLD:
+        raise ValueError('Saturation detected: {:.3f}'.format(peak))
 
-    data_bits = equalize(x, [F0, F1])
+    data_bits = equalize(x / amp, frequencies)
     if data_bits is None:
         log.info('Cannot demodulate symbols!')
     else:
         data = iterate(data_bits, bufsize=8, advance=8, func=to_bytes)
         data = ''.join(c for _, c in data)
-        log.info( 'Demodulated {} payload bydes'.format(len(data)) )
+        log.info( 'Demodulated {} payload bytes'.format(len(data)) )
         data = unpack(data)
         with file('data.recv', 'wb') as f:
             f.write(data)

send.py

@@ -29,13 +29,9 @@ def record(fname):
     return p
 
 
-log.info('MODEM Fc={}KHz, {} BAUD'.format(Fc/1e3, baud))
-
-
 class Symbol(object):
     t       = np.arange(0, Nsym) * Ts
-    c0      = np.exp(2j * np.pi * F0 * t)
+    carrier = [ np.exp(2j * np.pi * F * t) for F in frequencies ]
-    c1      = np.exp(2j * np.pi * F1 * t)
 
 sym = Symbol()
 
@@ -52,24 +48,30 @@ def train(sig, c):
         sig.send(c*0, n=10)
     sig.send(c*0, n=100)
 
+def modulate(sig, bits):
+    symbols_iter = sigproc.modulator.encode(list(bits))
+    symbols_iter = itertools.chain(symbols_iter, itertools.repeat(0))
+    carriers = np.array(sym.carrier) / len(sym.carrier)
+    while True:
+        symbols = itertools.islice(symbols_iter, len(sym.carrier))
+        symbols = np.array(list(symbols))
+        sig.send(np.dot(symbols, carriers))
+        if all(symbols == 0):
+            break
+
 if __name__ == '__main__':
 
+    bps = baud * sigproc.modulator.bits_per_symbol * len(sym.carrier)
+    log.info('Running MODEM @ {:.1f} kbps'.format(bps / 1e3))
+
     with open('tx.int16', 'wb') as fd:
         sig = Signal(fd)
-        start(sig, sym.c0)
+        start(sig, sym.carrier[carrier_index])
-        train(sig, sym.c0)
+        for c in sym.carrier:
-        train(sig, sym.c1)
+            train(sig, c)
-        carriers = [sym.c0, sym.c1]
 
         bits = to_bits(pack(data))
-        symbols_iter = sigproc.modulator.encode(list(bits))
+        modulate(sig, 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')

show.py

@@ -8,11 +8,15 @@ def spectrogram(t, x, Fs, NFFT=256):
     Pxx, freqs, bins, im = pylab.specgram(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)
+    for fname in sys.argv[1:]:
-    spectrogram(t, x, common.Fs)
+        t, x = common.load(fname)
+        pylab.figure()
+        pylab.title(fname)
+        spectrogram(t, x, common.Fs)
+
+    pylab.show()
+

sigproc.py

@@ -33,12 +33,11 @@ class Filter(object):
 
 
 class QAM(object):
-    def __init__(self, bits_per_symbol):
+    def __init__(self, bits_per_symbol, radii):
-        self._enc = {tuple(): 0.0} # End-Of-Stream symbol
+        self._enc = {}
         index = 0
-        amps = [0.6, 1.0]
+        N = (2 ** bits_per_symbol) / len(radii)
-        N = (2 ** bits_per_symbol) / len(amps)
+        for a in radii:
-        for a in amps:
             for i in range(N):
                 k = tuple(int(index & (1 << j) != 0) for j in range(bits_per_symbol))
                 v = np.exp(2j * i * np.pi / N)
@@ -61,7 +60,7 @@ class QAM(object):
             index = np.argmin(np.abs(s - keys))
             yield self._dec[ keys[index] ]
 
-modulator = QAM(bits_per_symbol=4)
+modulator = QAM(bits_per_symbol=4, radii=[0.6, 1.0])
 
 def test():
     q = QAM(bits_per_symbol=2)

test.sh

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