%% rtty modulator
clear all
fs = 8000;     % audio sample rate
dt = 1/fs;      % time per sample
fmid = 1000;
shift = 170;    % frequency deviation
baud = 45;      % bit per second
samp_per_bit = round(fs/baud);    % samples per bit (depending on fs)

data = [0 0 0 0 0 0 ...
    1 0 1 1 0 1 0 0 ...
    1 1 0 1 0 0 0 0 ...
    1 1 1 0 0 1 0 0 ...
    0 0 0 ...
    1 0 1 1 0 1 0 0 ...
    1 1 0 1 0 1 0 0 ...
    1 1 1 0 1 0 0 0 ...
    0 0 0 0 0 0 ...
    1 0 1 1 0 1 0 0 ...
    1 1 0 1 0 0 0 0 ...
    1 1 1 0 0 1 0 0 ...
    0 0 0 0 0 0 0];

f1 = fmid - shift/2;
f2 = fmid + shift/2;
f1_inc = 2*pi*f1/fs;
f2_inc = 2*pi*f2/fs;

t=(0:dt:length(data)/baud)';
y=zeros(length(data)*samp_per_bit,1);

for m=1:length(data)
    for n=1:samp_per_bit
        if data(m) == 0
            y(samp_per_bit*(m-1)+n+1) = y(samp_per_bit*(m-1)+n) + f1_inc;
        else
            y(samp_per_bit*(m-1)+n+1) = y(samp_per_bit*(m-1)+n) + f2_inc;
        end
    end
end
y = sin(y);
y = [zeros(100,1); y];
%y = sin(y) + sin(2*pi*1500*t); % cont. interference signal
%y = awgn(y, 20);               % white noise
%plot(t, y)                     % time domain plot

%% read rtty from wav
clear all
%rtty.wav = inverted
%rtty2.wav = not inverted
[y,fs] = wavread('rtty2.wav');
dt = 1/fs;
fmid = 1000;
shift = 170;
baud = 45;
inverted = 0;
y = awgn(y, 10);

%% rtty demodulator
offset_ideal = 2 * pi * 0.5*shift/fs;   % ideal DC value after FM demod for shift
nco_inc = 2 * pi * (fmid+100)/fs;
nco_pac = 0;
fil_len = 50;
fil_n = fil_len + 1;
RC = 0.01; % fg = 1 / 2 * pi * RC = 
alph = dt/(RC + dt);
IQf = fdesign.lowpass('N,Fc',fil_len, 0.5*fmid/fs);
IQd1 = design(IQf,'window','window',@hamming,'SystemObject',true);

bb = zeros(length(y),1);
bbf = zeros(length(y),1);
demod = zeros(length(y),1);
err = zeros(length(y),1);
errf = zeros(length(y),1);
samp_inst = zeros(length(y),1);
phi = zeros(2,1);

d = 0;
in_char = 0;
bit_cnt = 1;
bit_wait = 1;
byte_cnt = 1;
decode = 0;

for n=2:length(y)
    % iq mixing
    nco_pac = nco_pac + nco_inc;
    bb(n) = y(n) * exp(-1j*nco_pac);
    bbfd = bbf(n);
    if n > fil_n
        fil_tmp = filter(IQd1.Numerator,1,bb(n-fil_len:n));
        bbf(n) = fil_tmp(end);
    end
    phi(1) = phi(2);
    phi(2) = angle(bbf(n));
    demod(n) = diff(unwrap(phi));
    
    % error signal
    err(n) = sign(demod(n)) * (offset_ideal - abs(demod(n)));
    if isnan(err(n))
        err(n) = 0;
    end
    errf(n) = err(n);
    if n > fil_n
        rcf = errf(n-1) * (1-alph) + err(n) * alph;
        errf(n) = rcf;
    end
    nco_inc = nco_inc - errf(n)/100;
    
    % symbol recovery
    if inverted == 0
        mark = 1;
        space = 0;
    else
        mark = 0;
        space = 1;
    end
    dd = d;
    d = demod(n) > 0;
    if in_char == 0
        if d == space && d ~= dd % mark to space transition
            in_char = 1;
            bit_cnt = 0;
            bit_wait = 0.5 * fs/baud;
            samp_inst(n) = 1;
        end
    else
        bit_wait = bit_wait - 1;
        if bit_wait < 0
            if bit_cnt == 0
                if d == mark % start bit error
                    in_char = 0;
                    samp_inst(n) = -1;
                else 
                    bit_wait = fs/baud;
                    samp_inst(n) = 1;
                end
            elseif bit_cnt > 5
                byte_cnt = byte_cnt + 1;
                in_char = 0;
            else
                bit_wait = fs/baud;
                if inverted == 1
                    decode(byte_cnt, bit_cnt) = d;
                else 
                    decode(byte_cnt, bit_cnt) = ~d;
                end
                samp_inst(n) = 1;
            end
            bit_cnt = bit_cnt + 1;
        end     
    end
end

%plot(demod)
%hold on
%plot(errf,'r')
% lt = l
% lf = f
% cr = c
% nr = n
% sp = s
baudotchar = ['l' 'V' 'X' 'M' 'n' 'G' 'B' 'O' 'Q' 'P' 'Y' 'H' 'W' 'L' ...
    'Z' 'T' 'K' 'C' 'F' 'N' 'J' 'R' 'D' 'c' 'U' 'I' 'S' ' ' 'A' 'f' 'E' '*' ];
baudotnum = [ 'l' '=' '/' '.' 'n' '*' '?' '9' '1' '0' '6' '*' '2' ')' '+' ...
    '5' '(' ':' '*' ',' '*' '4' '*' 'c' '7' '8' '\' ' ' '-' 'f' '3' '*'];
let = 1;
for chr=1:length(decode)
    code = decode(chr,:);
    str = '';
    for i=1:5
       str = [num2str(code(i)) str];
    end
    if let == 1
        text(chr) = baudotchar(int8(bin2dec(str)) + 1);
    else
        text(chr) = baudotnum(int8(bin2dec(str)) + 1);
    end
    if text(chr) == 'l'
        let = 1;
    end
    if text(chr) == 'n'
        let = 0;
    end
end

text
