(This is still a prototype, and everything is subject to change; initial notes are in bleep-ultrasound-barcodes.ipynb.)
Bleep is a very simple modem for ultrasonic communication of cryptocurrency transactions of a few dozen bytes over distances of a few hundred millimeters; it can be used on any computer with a microphone, speaker, or ideally both, with reasonable response in the 17–20 kHz frequency range. Because of its minimal computational requirements, even small microcontrollers have enough power to encode and decode it.
Bleep uses continuous-phase frequency shift keying (CPFSK) with frequencies of 17640 Hz and 19110 Hz and a baud rate of 1470 baud. On top of this baud rate, it runs Manchester encoding, so it only gets 735 bits per second, which is further reduced by very simple ECC. (I haven't figured out what ECC to use yet.)
%pylab inline
import seaborn
style.use('seaborn-poster')
rcParams['lines.linewidth'] = 0.5
/usr/lib/python3.5/importlib/_bootstrap.py:222: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88 return f(*args, **kwds) /usr/lib/python3/dist-packages/matplotlib/__init__.py:874: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter. warnings.warn(self.msg_depr % (key, alt_key))
Populating the interactive namespace from numpy and matplotlib
44100/array([30, 17640, 19110, 1470])
array([1470. , 2.5 , 2.30769231, 30. ])
17640+1470
19110
2*3*5*7*7
1470
array([17640, 19110])/1470
array([12., 13.])
# That means they make 12 and 13 cycles in a 1470Hz (30 sample) period?
44100/array([17640, 19110])
array([2.5 , 2.30769231])
12 * 2.5, 13 * 2.30769231
(30.0, 30.000000030000002)
Yes. I think that means that over 30 samples they are perfectly orthogonal, but not over shorter periods of time.
fs = 44100
t = arange(60)/fs
c0 = sin(17640 * 2 * pi * t)
c1 = sin(19110 * 2 * pi * t)
correlations = (c0*c1).cumsum()
subplot(411); xticks(arange(0, 60, 2)); stem(c0)
subplot(412); xticks(arange(0, 60, 2)); stem(c1)
subplot(413); xticks(arange(0, 60, 2)); stem(correlations)
correlations
array([ 0.00000000e+00, 2.39073800e-01, 9.45846529e-01, 1.85035503e+00, 2.43492033e+00, 2.43492033e+00, 2.08942883e+00, 1.89169306e+00, 2.08942883e+00, 2.43492033e+00, 2.43492033e+00, 1.85035503e+00, 9.45846529e-01, 2.39073800e-01, -1.69309011e-15, -1.69309011e-15, -2.39073800e-01, -9.45846529e-01, -1.85035503e+00, -2.43492033e+00, -2.43492033e+00, -2.08942883e+00, -1.89169306e+00, -2.08942883e+00, -2.43492033e+00, -2.43492033e+00, -1.85035503e+00, -9.45846529e-01, -2.39073800e-01, 7.54951657e-15, 7.54951657e-15, 2.39073800e-01, 9.45846529e-01, 1.85035503e+00, 2.43492033e+00, 2.43492033e+00, 2.08942883e+00, 1.89169306e+00, 2.08942883e+00, 2.43492033e+00, 2.43492033e+00, 1.85035503e+00, 9.45846529e-01, 2.39073800e-01, -1.25732758e-14, -1.25732758e-14, -2.39073800e-01, -9.45846529e-01, -1.85035503e+00, -2.43492033e+00, -2.43492033e+00, -2.08942883e+00, -1.89169306e+00, -2.08942883e+00, -2.43492033e+00, -2.43492033e+00, -1.85035503e+00, -9.45846529e-01, -2.39073800e-01, 4.44089210e-15])
correlations[0], correlations[14], correlations[15], correlations[16], correlations[29], correlations[30]
(0.0, -1.6930901125533637e-15, -1.6930901125533712e-15, -0.23907380036690337, 7.549516567451064e-15, 7.549516567451094e-15)
c0[30], c1[30]
(-2.9391523179536475e-15, -1.028950903538412e-14)
c0[5:] - c0[:-5]
array([-4.89858720e-16, 3.33066907e-16, -1.11022302e-16, -2.22044605e-16, 3.33066907e-16, -4.89858720e-16, -2.44249065e-15, -2.22044605e-16, -1.11022302e-16, 4.44089210e-16, -4.89858720e-16, 3.21964677e-15, -1.11022302e-16, -2.22044605e-16, 4.44089210e-16, -4.89858720e-16, 4.44089210e-16, 1.99840144e-15, -2.33146835e-15, 3.33066907e-16, -4.89858720e-16, 3.33066907e-16, -2.33146835e-15, 2.10942375e-15, 4.44089210e-16, -4.89858720e-16, 4.44089210e-16, -1.11022302e-16, -2.22044605e-16, 3.33066907e-16, -4.89858720e-16, 4.44089210e-16, -1.11022302e-16, -1.11022302e-16, 4.44089210e-16, -4.89858720e-16, 3.33066907e-16, -2.22044605e-16, -1.11022302e-16, -1.11022302e-14, 1.37209960e-14, 1.18793864e-14, -1.11022302e-16, -2.22044605e-16, 1.18793864e-14, -1.47007134e-14, -2.25375274e-14, -2.22044605e-16, -4.55191440e-15, 4.44089210e-16, -1.47007134e-14, 1.18793864e-14, 4.32986980e-15, 4.32986980e-15, -1.11022302e-14])
c1[15:] + c1[:-15]
array([ 5.14475452e-15, 2.60902411e-15, 2.22044605e-15, -5.55111512e-16, -1.11022302e-16, -1.11022302e-16, 1.66533454e-15, -1.91513472e-15, -5.05151476e-15, 1.22124533e-15, -7.77156117e-16, 2.22044605e-16, -5.55111512e-16, 1.33226763e-15, -1.77635684e-15, -5.14475452e-15, 1.12132525e-14, 6.10622664e-15, -2.88657986e-15, 5.55111512e-16, -2.66453526e-15, 4.21884749e-15, -5.02375919e-15, -1.94289029e-15, 1.55431223e-15, -9.99200722e-16, 2.22044605e-16, 6.66133815e-16, -1.22124533e-15, 1.77635684e-15, -1.61715276e-14, 1.83186799e-15, 8.21565038e-15, -3.88578059e-15, 1.66533454e-15, -9.99200722e-16, 1.31006317e-14, -1.91513472e-15, -1.19904087e-14, -1.55431223e-15, -6.10622664e-15, -2.22044605e-16, 3.77475828e-15, 1.33226763e-15, -1.47659662e-14])
The 19110Hz signal requires a table of 8 separate values to specify it, which is the largest number that a period-30 even or odd function could need, but because 17640 Hz repeats exactly after 10 cycles, it only needs 3 separate values.
set(abs(c1.round(12))), set(abs(c0.round(12)))
({0.0, 0.207911690818, 0.406736643076, 0.587785252292, 0.743144825477, 0.866025403784, 0.951056516295, 0.994521895368}, {0.0, 0.587785252292, 0.951056516295})
By using a bit length of 30 samples (1470 baud) instead of 15, which would also have provided orthogonality (that's minimal-shift keying), the discontinuity in the derivative is smaller. By starting the wave at 0 (sin, not cos) we avoid a discontinuity when we're starting or ending a transmission.
mark = c1[:30]
space = c0[:30]
subplot(411); stem(mark); stem(space, markerfmt='r.')
mark, space
(array([ 0.00000000e+00, 4.06736643e-01, -7.43144825e-01, 9.51056516e-01, -9.94521895e-01, 8.66025404e-01, -5.87785252e-01, 2.07911691e-01, 2.07911691e-01, -5.87785252e-01, 8.66025404e-01, -9.94521895e-01, 9.51056516e-01, -7.43144825e-01, 4.06736643e-01, 5.14475452e-15, -4.06736643e-01, 7.43144825e-01, -9.51056516e-01, 9.94521895e-01, -8.66025404e-01, 5.87785252e-01, -2.07911691e-01, -2.07911691e-01, 5.87785252e-01, -8.66025404e-01, 9.94521895e-01, -9.51056516e-01, 7.43144825e-01, -4.06736643e-01]), array([ 0.00000000e+00, 5.87785252e-01, -9.51056516e-01, 9.51056516e-01, -5.87785252e-01, -4.89858720e-16, 5.87785252e-01, -9.51056516e-01, 9.51056516e-01, -5.87785252e-01, -9.79717439e-16, 5.87785252e-01, -9.51056516e-01, 9.51056516e-01, -5.87785252e-01, -1.46957616e-15, 5.87785252e-01, -9.51056516e-01, 9.51056516e-01, -5.87785252e-01, -1.95943488e-15, 5.87785252e-01, -9.51056516e-01, 9.51056516e-01, -5.87785252e-01, -2.44929360e-15, 5.87785252e-01, -9.51056516e-01, 9.51056516e-01, -5.87785252e-01]))
one = concatenate((mark, space))
zero = concatenate((space, mark))
subplot(611); plot(one)
subplot(612); plot(zero)
one
array([ 0.00000000e+00, 4.06736643e-01, -7.43144825e-01, 9.51056516e-01, -9.94521895e-01, 8.66025404e-01, -5.87785252e-01, 2.07911691e-01, 2.07911691e-01, -5.87785252e-01, 8.66025404e-01, -9.94521895e-01, 9.51056516e-01, -7.43144825e-01, 4.06736643e-01, 5.14475452e-15, -4.06736643e-01, 7.43144825e-01, -9.51056516e-01, 9.94521895e-01, -8.66025404e-01, 5.87785252e-01, -2.07911691e-01, -2.07911691e-01, 5.87785252e-01, -8.66025404e-01, 9.94521895e-01, -9.51056516e-01, 7.43144825e-01, -4.06736643e-01, 0.00000000e+00, 5.87785252e-01, -9.51056516e-01, 9.51056516e-01, -5.87785252e-01, -4.89858720e-16, 5.87785252e-01, -9.51056516e-01, 9.51056516e-01, -5.87785252e-01, -9.79717439e-16, 5.87785252e-01, -9.51056516e-01, 9.51056516e-01, -5.87785252e-01, -1.46957616e-15, 5.87785252e-01, -9.51056516e-01, 9.51056516e-01, -5.87785252e-01, -1.95943488e-15, 5.87785252e-01, -9.51056516e-01, 9.51056516e-01, -5.87785252e-01, -2.44929360e-15, 5.87785252e-01, -9.51056516e-01, 9.51056516e-01, -5.87785252e-01])
message = 'Come here, Mr. Johnson.'
bits = [1 if (1 << i) & ord(c) else 0 for c in message for i in range(8)]
subplot(811); plot(bits)
print(bits)
[1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0]
def box_filter(t, sig):
cs = sig.cumsum()
return concatenate((sig[:t], cs[t:] - cs[:-t]))
def forward_comb(t, w, sig):
return concatenate((sig[:t], sig[t:] + w * sig[:-t]))
sig = concatenate([one if bit else zero for bit in bits])
st = arange(len(sig)) / fs # signal time
subplot(811); plot(st, sig)
subplot(812); plot(st, forward_comb(5, -1, sig))
subplot(813); plot(st, forward_comb(15, 1, sig))
subplot(814); xlim(0, .01); plot(st, forward_comb(15, 1, sig))
subplot(212).set_yscale('log'); xlim(0, fs/2); ylim(.001, 10000); plot(abs(fft.fft(concatenate((sig, zeros(fs - len(sig)))))))
sig
array([ 0. , 0.40673664, -0.74314483, ..., -0.95105652, 0.74314483, -0.40673664])
efsig = forward_comb(1, -1, forward_comb(1, -1, forward_comb(1, -1, forward_comb(1, -1, sig)))) # emission-filtered sig
subplot(811); plot(st, efsig)
subplot(812); plot(st, forward_comb(5, -1, efsig))
subplot(813); plot(st, forward_comb(15, 1, efsig))
subplot(814); xlim(0, .01); plot(st, forward_comb(15, 1, efsig))
subplot(212).set_yscale('log'); xlim(0, fs/2); ylim(.001, 10000); plot(abs(fft.fft(concatenate((efsig, zeros(fs - len(sig)))))))
efsig
array([ 0. , 0.40673664, -2.3700914 , ..., -12.68434922, 14.56638913, -13.92976803])
subplot(811); xlim(0, .01); plot(st, forward_comb(5, -1, efsig))
subplot(812); xlim(0, .01); plot(st, forward_comb(15, 1, efsig))
subplot(813); xlim(0, .01); plot(st, efsig * (-1)**arange(len(efsig)))
subplot(814); xlim(0, .01); plot(st, efsig)
subplot(413).set_yscale('log'); ylim(.1, 100000); xlim(0, fs/2); plot(abs(fft.fft(resize(efsig * (-1)**arange(len(efsig)), fs))))
[<matplotlib.lines.Line2D at 0x7fc3a24fa0b8>]
subplot(811); plot(st, abs(box_filter(30, resize(mark, len(efsig))*efsig)))
subplot(812); plot(st, abs(box_filter(30, resize(space, len(efsig))*efsig)))
subplot(813); xlim(0, .01); plot(st, abs(box_filter(5, box_filter(30, resize(mark, len(efsig))*efsig))))
subplot(814); xlim(0, .01); plot(st, abs(box_filter(30, resize(space, len(efsig))*efsig)))
subplot(815); xlim(0, .01); plot(st, resize(mark, len(efsig))*efsig)
subplot(816); xlim(0, .01); plot(st, resize(space, len(efsig))*efsig)
subplot(414).set_yscale('log'); ylim(.1, 100000); xlim(0, fs/2); plot(abs(fft.fft(resize(efsig * resize(space, len(efsig)), fs))))
st
array([0.00000000e+00, 2.26757370e-05, 4.53514739e-05, ..., 2.50272109e-01, 2.50294785e-01, 2.50317460e-01])
That's not working because the emission filter is introducing a phase shift.
subplot(811); stem((resize(space, len(efsig))*efsig)[:60])
subplot(812); plot((resize(space, len(efsig))*efsig).cumsum()[:120])
subplot(813); plot(forward_comb(1, -1, (resize(space, len(efsig))*efsig).cumsum()[:1000:30]))
[<matplotlib.lines.Line2D at 0x7fc3a1f3ee10>]
The emission-filter kernel, intended to high-pass the signal to eliminate audible components, is a cascade of four backward-differencing operations, so its impulse response is the fourth row of Pascal's Triangle with alternating signs (a fairly precise approximation of a Gabor wavelet at the Nyquist frequency, amusingly), and its frequency response is just a straight high-pass with 24dB per octave (80dB per decade) and no cutoff frequency. This amplifies our carrier signals by about 20dB relative to the limit of my personal hearing, which is about 10kHz, but by less than 10dB relative to 15kHz, which some people can hear. But our CPFSK signal itself has already dropped off by about 20dB by the time you get down to 15kHz.
impulse = zeros(fs)
impulse[0] = 1
efkern = forward_comb(1, -1, forward_comb(1, -1, forward_comb(1, -1, forward_comb(1, -1, impulse))))
ef_fr = fft.fft(efkern)
subplot(211); stem(efkern[:8])
subplot(212).set_yscale('log'); xlim(0, fs/2); ylim(.0001, 100); plot(abs(ef_fr))
efkern[:8], [(freq, abs(ef_fr[freq]), ef_fr[freq]) for freq in [17640, 19110, 15000, 10000]]
(array([ 1., -4., 6., -4., 1., 0., 0., 0.]), [(17640, 13.090169943749475, (4.045084971874737+12.449491424413901j)), (19110, 14.646624873858524, (9.800504982955689+10.884563485716455j)), (15000, 9.443923661612029, (-4.006411580916092+8.551979909392735j)), (10000, 2.9205496465420064, (-2.7968594508924522-0.8409443797720568j))])
Still, maybe it would be better to use an even more aggressive emission filter. Just going up to a cascade of 8 differentiators doesn't help a lot. It widens the Gaussian window by a factor of $\sqrt 2$.
def bd(sig): return forward_comb(1, -1, sig)
efkern2 = bd(bd(bd(bd( bd(bd(bd(bd(impulse)))) ))))
ef_fr2 = fft.fft(efkern2)
subplot(211); stem(efkern2[:16])
subplot(212).set_yscale('log'); xlim(0, fs/2); ylim(.0001, 1000); plot(abs(ef_fr2))
efkern2[:16], [(freq, abs(ef_fr2[freq]), ef_fr2[freq]) for freq in [17640, 19110, 15000, 10000]]
(array([ 1., -8., 28., -56., 70., -56., 28., -8., 1., 0., 0., 0., 0., 0., 0., 0.]), [(17640, 171.3525491562421, (-138.62712429686843+100.71850133676017j)), (19110, 214.5236201955312, (-22.423824353652627+213.3484373581233j)), (15000, 89.18769412635558, (-57.085026614958416-68.52550269750563j)), (10000, 8.529610237916655, (7.11523533817622+4.704006472480768j))])
I don't know. I have to think about that.
How about feedforward comb filters for filtering out random noise, as opposed to just distinguishing the two carriers? Let's see what the carrier sample sequences' autocorrelation functions look like.
def cyclic_correlation(sig1, sig2):
return array([sum(sig1 * concatenate((sig2[i:], sig2[:i]))) for i in range(len(sig2))])/len(sig2)
mark_acf = cyclic_correlation(mark, mark)
space_acf = cyclic_correlation(space, space)
subplot(811); stem(mark_acf)
subplot(812); stem(mark)
subplot(813); stem(space_acf)
subplot(814); stem(space)
<Container object of 3 artists>
Hmm, and I guess another thing to look at is maybe what the correlations look like with the backward differences, since backward differences give us half a cycle of phase shift, which could be useful.
mark_dacf = cyclic_correlation(mark, bd(mark))
space_dacf = cyclic_correlation(space, bd(space))
subplot(811); stem(mark_dacf)
subplot(812); stem(space_dacf)
<Container object of 3 artists>
This is fucking bullshit, man. What I want is to know which feedforward comb filters give us the best response at the fucking carrier frequencies. These aren't fucking comb filters. This is a fucking comb filter.
subplot(811); stem(forward_comb(3, -1, space))
<Container object of 3 artists>
Except it's not cyclic; it's using zeroes before the beginning.
def cyclic_comb(n, w, sig):
return sig + w * concatenate((sig[n:], sig[:n]))
subplot(811); stem(cyclic_comb(3, -1, space))
<Container object of 3 artists>
Now let's look at the power gains of different unity-gain feedforward comb filters.
def comb_power_gains(w, sig):
return array([(cyclic_comb(i, w, sig)**2).mean() for i in range(len(space))])
subplot(811); stem(comb_power_gains(1, space))
subplot(812); stem(comb_power_gains(-1, space))
subplot(813); stem(comb_power_gains(1, mark))
subplot(814); stem(comb_power_gains(-1, mark))
<Container object of 3 artists>
So it looks like we can amplify the space signal by adding it to itself at lags of 2, 3, or fortunately 5, or subtracting it at lags of 1 and 4; and we can amplify the mark signal by adding it to itself at gains of 2, 5, 7, and 9, or by subtracting it at gains of 1, 6, and 8. What if we try five of those? What does the impulse response look like?
mark_kern = forward_comb(2, 1, forward_comb(5, 1, forward_comb(7, 1, forward_comb(1, -1, forward_comb(6, -1, impulse)))))
space_kern = forward_comb(2, 1, forward_comb(3, 1, forward_comb(5, 1, forward_comb(1, -1, forward_comb(4, -1, impulse)))))
mark_k_fs = fft.fft(mark_kern)
space_k_fs = fft.fft(space_kern)
subplot(811); stem(mark_kern[:30])
subplot(812); stem(mark)
subplot(813); stem(space_kern[:30])
subplot(814); stem(bd(space))
subplot(413).set_yscale('log'); ylim(1, 100); xlim(0, fs/2); plot(abs(mark_k_fs))
subplot(414).set_yscale('log'); ylim(1, 100); xlim(0, fs/2); plot(abs(space_k_fs))
space_kern[:30], mark_kern[:30], abs(mark_k_fs[:int(fs/2)]).argmax(), abs(space_k_fs[:int(fs/2)]).argmax()
(array([ 1., -1., 1., 0., -2., 3., -3., 1., 1., -3., 3., -2., 0., 1., -1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]), array([ 1., -1., 1., -1., 0., 1., -2., 3., -3., 2., -1., -1., 2., -3., 3., -2., 1., 0., -1., 1., -1., 1., 0., 0., 0., 0., 0., 0., 0., 0.]), 18602, 17260)
Well, so the frequency peaks are in roughly the right places, but the frequency selectivity is kind of crap. In fact it's giving you less than 3dB of frequency separation, which is sort of okay in the sense that we're relying on a separate comb filter for frequency separation, as shown below, but it also sort of implies that random noise over a fairly wide frequency range is going to cause interference.
subplot(211); xlim(17000, 20000); xticks([17000, 17500, 17640, 18000, 18500, 19000, 19110, 20000]);
plot(abs(mark_k_fs)/abs(mark_k_fs).max(), label="mark"); plot(abs(space_k_fs)/abs(space_k_fs).max(), label="space")
legend()
subplot(413); xlim(17000, 20000); xticks([17000, 17500, 17640, 18000, 18500, 19000, 19110, 20000]);
plot(abs(fft.fft(forward_comb(5, -1, impulse)))); plot(abs(fft.fft(forward_comb(15, 1, impulse))))
subplot(414); xlim(0, fs/2); plot(abs(fft.fft(forward_comb(5, -1, impulse)))); plot(abs(fft.fft(forward_comb(15, 1, impulse))))
mark_k_fs
array([ 0.00000000e+00+0.00000000e+00j, -9.74368443e-07+1.45765287e-09j, -3.89745800e-06+1.16612019e-08j, ..., -8.76922140e-06-3.93564383e-08j, -3.89745800e-06-1.16612026e-08j, -9.74368441e-07-1.45765187e-09j])
That's to be expected, since our kernels here have a very short time response. The solution, as I described in sparse-bandpass-filters.ipynb, is to use unity-gain feedback comb filters, in fixed point, carefully windowed so as to have a finite impulse response despite the feedback. In this case we should be able to use filters with a Q of at least 8, since we have 12 or 13 cycles within a given bit time, and that will get us a bandpass of about 2500 Hz bandwidth between 3dB points... but that may not be much better than what we're getting here actually.
Block codes would seem to complicate the software interface somewhat; either you require the transmission of an entire block at once, or you buffer data, making the user wonder why it hasn't arrived, and requiring a separate flush operation. So convolutional codes might be preferable. Systematic codes are probably preferable to non-systematic codes, because it means that a hacked-together implementation can start by ignoring the ECC bits. A really simple convolutional error-correction code that might be adequate without further work is a diagonal parity code: you append two or three parity bits to each byte. The first parity bit on byte $n$ is the (complemented) XOR of the current byte $\sim\bigoplus_{i=0}^7b_{n,i}$ indicates whether the byte has been received correctly, while the others are convolutional parity bits over the last M bytes. For example, you could have a diagonal $\bigoplus_{i=0}^7b_{n-i-1,i}$ bit, taking one bit from each of the previous 8 bytes; this would enable the reconstruction not only of any erroneous bit but also of any entirely lost byte without necessarily imposing any latency in the error-free case. If you added a $\bigoplus_{i=0}^7b_{n-2i-2,i}$ bit, drawing from the previous 17 bytes, or a $\bigoplus_{i=0}^7 b_{n-i-1,3i\%8}$ bit drawing from the previous 8 bytes, these would also be orthogonal to the per-byte parity and to the previous diagonal parity, in the sense that they have only one bit of overlap.
These schemes are not only trivial to encode but also pretty simple to correct errors with; you can actually do it by hand.
My biggest concern here is that a single parity bit is not strong enough to detect all two-bit errors, and indeed it has a 50% chance of incorrectly accepting a random byte. But any other single-bit scheme either is redundant with the parity bit or ignores some bits, at least sometimes. Some kind of better two-bit checksum that can detect all two-bit errors (i.e. a (10, 8) double-error-detection code with a minimal Hamming distance of 3) would be nice, but I'm pretty sure that's impossible as a linear code. You could of course use a (15, 11) Hamming code shortened to (12, 8), but that seems unreasonable. My reason for wanting a strong checksum on each byte is so that error-free bytes can be passed through immediately, or at least quickly.
This scheme occurred to me, and I wanted to see how well it works.
%%time
impulse = array([1.0])
impulse.resize(fs)
def rhp(n, sig):
"Recursive high-pass filter."
buf = sig.copy()
for i in range(n, len(buf)):
o = buf[i-n+1] - buf[i-n]
buf[i] += o/2 - o/8
return buf
# 8 or 38 gets a sharpish peak at 19110
# 31 or 61 has sharpish peaks just to the right of 17640 and 19110
rhp8 = rhp(8, impulse)
rhp8_f = fft.fft(rhp8)
CPU times: user 240 ms, sys: 0 ns, total: 240 ms Wall time: 240 ms
subplot(511); xlim(0, 200); plot(rhp8, '.')
subplot(512).set_yscale('log'); xlim(0, 200); ylim(.0001, 1); plot(rhp8, '.')
subplot(513); xlim(0, 90); plot(rhp8, '.')
subplot(514).set_yscale('log'); xlim(0, fs/2); plot(abs(rhp8_f))
subplot(515).set_yscale('log'); xlim(17000, 20000)
xticks([17000, 17500, 17640, 18000, 18500, 19000, 19110, 20000]); plot(abs(rhp8_f))
abs(rhp8_f[[17639, 17640, 17641, 19109, 19110, 19111]])
array([0.77604918, 0.77613132, 0.77621373, 2.39433322, 2.39490372, 2.39544802])
Above we can see that, although it does separate the signals significantly, its long exponential decay is a problem in that it smears the bits out in time quite a lot. To some extent this is just a consequence of the high Q we can see in the plot above. This is not really the best test, because the reason I was thinking this thing might be a good idea is not primarily to separate the two carrier signals from each other but from noise at other frequencies, and we can't really see how much of an improvement we're getting here.
subplot(811); plot(st, efsig)
subplot(812); plot(st, rhp(8, efsig))
subplot(813); xlim(0, .01); plot(st, rhp(8, efsig))
rhp(8, efsig)
array([ 0. , 0.40673664, -2.3700914 , ..., -35.46065241, 42.22405716, -41.67971736])
Another thing I'm thinking about is demodulation using modified square waves, sometimes fraudulently called "modified sine waves" by UPS vendors. The idea is that the waveform [1, 0, -1, 0] has the same min, max, mean, and RMS (standard deviation) as a sinewave, and indeed if sampled at only those points is identical to a sinewave, and so is an adequate substitute for some purposes. These properties remain unchanged if we resample it to a higher sample rate using nearest-neighbor resampling instead of sinc resampling.
msw0 = array([1, 0, -1, 0])
sw = cos(arange(len(msw0)) * 2 * pi / len(msw0))
subplot(411); plot(sw, 'o-'); plot(msw0, '.-');
[(w.min(), w.max(), w.mean(), w.std()) for w in [msw0, sw]]
[(-1, 1, 0.0, 0.7071067811865476), (-1.0, 1.0, -4.592425496802574e-17, 0.7071067811865476)]
So what does the spectrum of, say, a 3kHz modified square wave look like?
def msw(n):
x = arange(n)
return where(x < n/8, 1, where(x < 3*n/8, 0, where(x < 5*n/8, -1, where(x < 7*n/8, 0, 1))))
msw3k = msw(fs / 3000)
subplot(411); stem(msw3k)
subplot(412).set_yscale('log'); xlim(0, fs/2); xticks(arange(0, fs/2, 3000)); plot(abs(fft.fft(resize(msw3k, fs))))
subplot(413).set_yscale('log'); xlim(0, fs/2); ylim(100, 100000); plot(abs(fft.fft(resize(msw3k, fs))))
msw3k
array([ 1, 1, 0, 0, 0, 0, -1, -1, -1, -1, 0, 0, 0, 1, 1])
So, the third harmonic is attenuated, but only by about 20dB, and higher harmonics are just about as bad; and the filters are extremely sharp because they continue for the entire second. But the bigger thing to note here is that the filter frequencies are slightly off, because 44100/15 ≠ 3000, but rather 2940. And because they are so sharp, they pick up basically nothing at 3000. But if we window them to a narrowish time period, it's fine. Here we use a Hamming window to avoid worrying about whether the window we're using is causing unnecessary problems, before worrying about its implementation efficiency.
msw3k_4 = resize(msw3k, 60).copy() * hamming(60)
msw3k_4.resize(fs)
subplot(411).set_yscale('log'); xlim(0, fs/2); ylim(.01, 20); xticks(arange(0, fs/2, 3000)); plot(abs(fft.fft(msw3k_4)))
msw3k_4
array([0.08 , 0.08260599, 0. , ..., 0. , 0. , 0. ])
Hmm, here's an interesting thought: $y(n) = x(n) - x(n-5)$ cancels the 17640Hz wave because it repeats with a period of 2½, and $y(n) = \sum_{i=0}^{i<5} x(n-i)$ is a smoothing thing which also cancels it, which can be calculated with $y(n) = x(n) - x(n-5)$ on a prefix sum. I was thinking that the convolution of [1, 1] with [-1, 1] was [-1, 0, 1] and so you could calculate $y(n) = x(n) - x(n-10)$ as a faster equivalent to the composition of $y(n) = x(n) + x(n-5)$ and $y(n) = x(n) - x(n-5)$. But where does the $+$ come from?
Oh, it's because once we invert every other sample, the 17640Hz wave becomes a 4410Hz wave, with a period of 10 samples, so the way to cancel it is exactly with $y(n) = x(n) + x(n-5)$. So this allows you to cancel that wave and also get some smoothing with $y(n) = x(n) - x(n-10)$, which you can do on a 5-decimated sum-table signal as $y(n) = x(n) - x(n-2)$. Unfortunately 5-decimation is too extreme to allow us to also recover the 17640Hz wave with arbitrary phase, although an aliased version of the 19110Hz wave (downconverted to 2940Hz) should survive, so perhaps we can do the detection of the 17640Hz wave on a separately decimated version from the same cascade of integrators and perhaps differentiators.
This is maybe getting a bit ahead of my explanation, but what I mean is, we can run the signal first through some differentiators to filter out everything that isn't ultrasonic, then swap every other sample and integrate two or three times, and then use two separate decimated versions of that sum-table signal to attempt to filter out the original carriers. The 19110Hz wave can be recovered from a 5-fold decimated sum table as described above, plus a $y(n) = x(n) - x(n-3)$ comb or two, while the 17640Hz wave can be recovered from, say, a 3-fold decimated sum table, using $y(n) = x(n) - x(n-5)$ to cancel out the 19110Hz signal and then maybe a $y(n) = x(n) - x(n-1)$ differentiator or two.
You might think that the differentiation and the integration at the front end would cancel out, and normally they would, but the sample-swapping is not time-invariant (though linear), so the usual commutative properties of LTI signals don't hold here. In fact, I think that perhaps differentiation before sample-swapping may be the same thing as integration after it, which leads me to worry about overflow due to noise that's exactly at Nyquist. Hmm, no, can't be — differentiation is still FIR, integration is still IIR, just tamed by the following comb filters.
def get19110(sig):
out = zeros(len(sig) // 5)
i3out = [0.0, 0.0]
d1 = d2 = d3 = d4 = i1 = i2 = i3 = c1 = c2 = c3 = 0.0
for n, xn in enumerate(sig):
d1, xn = xn, xn - d1
d2, xn = xn, xn - d2
d3, xn = xn, xn - d3
d4, xn = xn, xn - d4
i1 += xn if n % 2 else -xn
i2 += i1
i3 += i2
if not n % 5:
c1 = i3 - i3out[-2]
i3out.append(i3)
c2, xn = c1, c1 - c2
c3, out[n // 5] = xn, xn - c3
return out
ef19110 = get19110(efsig)
subplot(811); plot(ef19110)
subplot(812); plot(ef19110[:200], '.-')
subplot(813); stem(ef19110[:40])
subplot(814); plot(efsig[:2000])
ef19110
array([ 0. , 417.88477782, 10161.4511827 , ..., 14283.72154999, -126.53174591, -13426.44116592])
That seems pretty reasonable; it's interleaved I and Q channels, plus their negatives. It's anyone's guess what its noise immunity will look like, and because the sample-flipping isn't LTI, we can't just look at the impulse response. Well, I guess we can add some noise to see. It looks kind of okay:
random.seed(0)
efnoise = random.random(len(efsig)) * 10 + efsig
efn19110 = get19110(efnoise)
efi = efn19110[::4] - efn19110[2::4]
efq = efn19110[1::4] - efn19110[3::4]
subplot(811); plot(efn19110)
subplot(812); plot(efn19110[:200], '.-')
subplot(813); stem(efn19110[:40])
subplot(814); plot(efnoise[:2000])
subplot(413).set_yscale('log'); xlim(0, fs/2); ylim(10, 1000000); plot(abs(fft.fft(resize(efnoise, fs))))
subplot(817); stem(efi[:50])
subplot(818); stem(efq[:50])
efn19110, (efsig**2).sum(), ((efnoise-efsig)**2).sum()
(array([-5.48813504e+00, 2.06788808e+01, 1.01513412e+04, ..., 1.13774188e+04, -3.42441587e+03, -1.51195642e+04]), 1065297.149901005, 363166.78655722574)