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

data = [0 0 0 1 0 1 1 0 1 0 0 1 1 0 1 0 1 0 0 1 1 1 0 0 1 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) == 1
            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 = sin(y) + sin(2*pi*1500*t); % cont. interference signal
%y = awgn(y, 20);               % white noise
%plot(t, y)                     % time domain plot

%% rtty demodulator
offset_ideal = 2 * pi * 0.5*shift/fs

bb = y.*exp(-1j*2*pi*fmid.*t);
Hf = fdesign.lowpass('N,Fc',50, 0.5*fmid/fs);
Hd1 = design(Hf,'window','window',@hamming,'SystemObject',true);
bbf = filter(Hd1.Numerator,1,bb);

i = real(bbf);
q = imag(bbf);
i = i./abs(bbf);
q = q./abs(bbf);
demod = i(2:end).*q(1:end-1) - i(1:end-1).*q(2:end);
plot(demod)




