Make a (software) modem!

Question

Objective

Design a modulator/demodulator pair to accurately transmit data as quickly as possible over simulated plain old telephone service (POTS).

Steps

1. Generate some random (/dev/random or the like) data that will take 3-4 seconds to transmit
2. Modulate the data with your modulator to produce an audio file
3. Pass the audio file through the POTS simulator. If you don't have Python/Scipy you can upload a file with the form, or do a JSON API request.
4. Demodulate the audio file back to binary data
5. Validate that the input and output are equal-ish* (limit 1 out of every 1000 bits can be corrupted)
6. Score is number of bits transmitted divided by the length of the audio file (bits/second)

Rules

• Input file must be 3-4 seconds, 44.1 kHz, mono.
• Run the simulator with an SNR of 30 dB (it's default)
• The demodulator must reconstruct the transmitted data with a bit error rate of no more than 10^-3 (1 per thousand bits).
• No digital compression is allowed (i.e. zipping the data. It's outside the scope of the challenge.)
• No attempting to shove data into frequencies above 4 kHz. (My filters aren't perfect, but they're reasonably POTS-like with a relatively small number of taps.)
• If your modem protocol requires a short preamble (no more than 1 second) to synchronize/calibrate the receiver, it isn't penalized.
• If possible, please host the audio file somewhere accessible so we can listen to a cacophony of beeps and boops.

Example

Here's an example notebook that demonstrates the modulation/demodulation with simple "on-off keying" (audio samples included!).

It would score 100 (bits/second). Note that it's transmitting with a much worse 5 dB SNR.

Answer 1

MATLAB, 1960 bps

Here is my updated attempt:

fs = 44100; %44.1kHz audio rate
fc = 2450;  %2.45kHz carrier - nice fraction of fs!
fsym = fc/5; %symbol rate

tmax = 4; %about 4 seconds worth

preamblesyms = 6;

t = 1/fs:1/fs:(tmax+preamblesyms/fsym);

symbols = preamblesyms+fsym*tmax;
symbollength = length(t)/symbols;
bits = symbols*3;
bitstream = [zeros(1,preamblesyms*3),rand(1,bits-preamblesyms*3)>0.5]; %Add a little preamble of 18 bits
data = bin2dec(char(reshape(bitstream,3,symbols)'+'0'))';

greycode = [0 1 3 2 6 7 5 4];

%Encode the symbols using QAM8 - we use effectively grey code so that
%adjacent symbols in the constellation have only one bit difference
%(minimises error rate)
encoded = zeros(2,symbols);
encoded(1,data==1) = 1/sqrt(2);
encoded(1,data==3) = 1;
encoded(1,data==2) = 1/sqrt(2);
encoded(1,data==7) = -1/sqrt(2);
encoded(1,data==5) = -1;
encoded(1,data==4) = -1/sqrt(2);
encoded(2,data==0) = 1;
encoded(2,data==1) = 1/sqrt(2);
encoded(2,data==2) = -1/sqrt(2);
encoded(2,data==6) = -1;
encoded(2,data==7) = -1/sqrt(2);
encoded(2,data==4) = 1/sqrt(2);

%Modulate onto carrier
carrier = [sin(2*pi*fc*t);cos(2*pi*fc*t)];
signal = reshape(repmat(encoded(1,:)',1,symbollength)',1,[]);
signal(2,:) = reshape(repmat(encoded(2,:)',1,symbollength)',1,[]);
modulated = sum(signal.*carrier)';

%Write out an audio file
audiowrite('audio.wav',modulated,fs);

%Wait for the user to run through the POTS simulator
input('');

%Read in the filtered data
filtered=audioread('audio.pots-filtered.wav')';

%Recover the two carrier signals
preamblecos = filtered(symbollength+1:symbollength*2);
preamblesin = filtered(symbollength+1+round(symbollength*3/4):symbollength*2+round(symbollength*3/4));

%Replicated the recovered carriers for all symbols
carrierfiltered = [repmat(preamblesin,1,symbols);repmat(preamblecos,1,symbols)];

%Generate a demodulation filter (pass up to 0.66*fc, stop at 1.33*fc
%(really we just need to kill everything around 2*fc where the alias ends up)
d=fdesign.lowpass('Fp,Fst,Ap,Ast',0.05,0.1,0.5,60);
Hd = design(d,'equiripple');

%Demodulate the incoming stream
demodulated = carrierfiltered .* [filtered;filtered];
demodulated(1,:)=filtfilt(Hd.Numerator,1,demodulated(1,:));
demodulated(2,:)=filtfilt(Hd.Numerator,1,demodulated(2,:));

%Split signal up into bit periods
recovereddemodulated=[];
recovereddemodulated(1,:,:) = reshape(demodulated(1,:),symbollength,symbols);
recovereddemodulated(2,:,:) = reshape(demodulated(2,:),symbollength,symbols);

%Extract the average level for each bit period. Only look at the second
%half to account for slow rise times in the signal due to filtering
recoveredsignal=mean(recovereddemodulated(1,round(symbollength/2):symbollength,:));
recoveredsignal(2,:)=mean(recovereddemodulated(2,round(symbollength/2):symbollength,:));

%Convert the recovered signal into a complex number.
recoveredsignal=recoveredsignal(2,:) + 1j*recoveredsignal(1,:);

%Determine the magnitude and angle of the symbol. The phase is normalised
%to pi/4 as that is the angle between the symbols. Rounding this to the
%nearest integer will tell us which of the 8 phases it is closest to
recoveredphase = round(angle(recoveredsignal)/(pi/4));
recoveredphase = mod(recoveredphase+8,8)+1; %Remap to an index in the grey code vector.

%Determine the symbol in the QAM8 constellation
recoveredencoded=greycode(recoveredphase);
recoveredencoded(1:preamblesyms)=0; %Assume the preamble is correct for comparison

%Turn it back in to a bit stream
bitstreamRecovered = reshape(dec2bin(recoveredencoded)'-'0',1,[]);

%And check if they are all correct...
if(all(bitstream==bitstreamRecovered))
    disp(['Woop, ' num2str(fsym*4) 'bps']);
else
    error('Its corrupt Jim.');
end

Since my first attempt, I have played around a bit. There is now a small preamble at the beginning (18 bit periods, but could be shorter) which contains just a cosine wave. I extract this and replicated it to create correctly phased sine and cosine carriers for demodulation - as it is a very short preamble, I haven't counted it in the bit rate as per your instructions.

Also since the first attempt I am now using a QAM8 constellation to achieve 3 bits per symbol instead of 2. This effectively doubles the transfer rate. So with a ~2.4kHz carrier I am now achieving 1960bps.

I've also improved the symbol detection so that the averaging doesn't get affected by slow rise times caused by the filtering - basically only the second half of each bit period is averaged to remove the impact of rise times.

Still nowhere near the 40kbps theoretical channel bandwidth from the Shannon-Hartley theory (assuming the 30dB SNR)

Just for those that like horrible sounds, this is the new entry:

---

And in case anyone is interested, this is the previous 960bps entry