clear;clc;
%% Declaration
% msg = 'A1B2 .?C3D4E5'; %message1
% msg = 'OK1FLZ ABCDEF'; %message2
msg = 'OK1FLZ JN79IA'; %message3
msgLength = length(msg);
if msgLength < 13
    for i=1:(13-msgLength), msg(msgLength+i) = ' '; end
    msgLength = length(msg);
elseif msgLength > 13, disp("msg too long"); return; end
alphabet = ['0','1','2','3','4','5','6','7','8','9','A','B','C','D','E',...
    'F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T','U','V',...
    'W','X','Y','Z',' ','+','-','.','/','?'];
alphabetLength = length(alphabet);
%% Source coding
%for every character in message we use its index in "alphabet" to encode it
msgCoded = zeros(1,msgLength);
for i=1:msgLength
    for j=1:alphabetLength
        if msg(i)==alphabet(j)
            msgCoded(i)=j-1;
        end
    end
end
clear i j alphabetLength alphabet;
%% Forming long integers
N1 = (((msgCoded(1)*42+msgCoded(2))*42+msgCoded(3))*42+msgCoded(4))*42+msgCoded(5);
N2 = (((msgCoded(6)*42+msgCoded(7))*42+msgCoded(8))*42+msgCoded(9))*42+msgCoded(10);
N3 = (msgCoded(11)*42+msgCoded(12))*42+msgCoded(13);
%% Rearanging long integers
N1 = bitset(bitshift(N1,1),1,mod(round(N3/32768),2));
N2 = bitset(bitshift(N2,1),1,mod(round(N3/65536),2));
N3 = mod(N3,32768) + 32768;
%adding 32768=2^15 is for setting bit 15 to 1, that marks message type
%as plain text for WSJT-X software
%% Bit sequence
%forming 103 bit long sequence for convolutional encoding
bitSeq = [dec2bin(N1,28),dec2bin(N2,28),dec2bin(N3,16),dec2bin(0,31)];
bitSeqLength = length(bitSeq);
%% Convolution encoding
reg0 = uint32(0);
reg1 = uint32(0);
H0 = 0xF2D05351; %defined in JT9 specification
H1 = 0xE4613C47; %defined in JT9 specification
bitOut = zeros(1,206);
bitCount = 1; %output bits counter
%for every iteration of cycle one input bit is shifted to 32bit registers
%reg0,1 then by AND with H0 and XOR through resulting 32bits even parity bit
%is obtained and append to output, same process with H1, thus for one input
%bit two output bits are obtained
for i=1:bitSeqLength
        reg0 = bitshift(reg0,1);
        reg0 = bitset(reg0,1,str2double(bitSeq(i)));
        bitAnd1 = bitand(reg0,H0);
        parity1 = 0;
        for k=1:32
            parity1 = bitxor(parity1,bitand(bitAnd1,1));
            bitAnd1 = bitshift(bitAnd1,-1);
        end
        bitOut(bitCount) = parity1;
        bitCount = bitCount + 1;
        
        reg1 = bitshift(reg1,1);
        reg1 = bitset(reg1,1,str2double(bitSeq(i)));
        bitAnd2 = bitand(reg1,H1);
        parity2 = 0;
        for k=1:32
            parity2 = bitxor(parity2,bitand(bitAnd2,1));
            bitAnd2 = bitshift(bitAnd2,-1);
        end
        bitOut(bitCount) = parity2;
        bitCount = bitCount + 1;
end
clear parity1 parity2 bitCount bitAnd1 bitAnd2 reg0 reg1 H0 H1;
%% Interleaving
%"mapping" contains new indexes for interleaving "bitOut"
bitDest(1,1:207) = uint8(0);
mapping(1,1:206) = uint8(0);
k = 1;
for i=0:254
    n = 0;
	for j=0:7, n = bitshift(n,1) + bitand(bitshift(i,-j),1); end
	if n < 206, mapping(k) = n; k=k+1; end
	if k > 206, break, end
end
for i=1:206
    bitDest(mapping(i)+1) = bitOut(i);
end
clear j0 n k i;
%% Convert to symbols
symbols = zeros(1,69);
j = 1;
for i=1:69
    symbols(i) = bitDest(j)*4+bitDest(j+1)*2+bitDest(j+2);
    j = j+3;
end
clear i j;
%% Gray Coding
symbolsGray = bin2gray(symbols,'fsk',8);
%% Sync
syncVector = '1100100001000001000000100000000010100000000000000011001000010000010000001000000000101';
symbolsOut = zeros(1,85);
j = 1;
for i=1:length(syncVector)
    if syncVector(i)=='1'
        symbolsOut(i) = 0;
    else
        if syncVector(i)=='0'
            symbolsOut(i) = symbolsGray(j)+1; %adding 1 to shift data values from 0-8 to 1-9
            j = j+1;
        end
    end
end
clear i j
%% Modulation 9-FSK
%for every channel symbol in sequence we assign waveform with
%corresponding frequency
Fs = 12000;
duration = 6912;
M = 9;
freqSep = Fs/duration;
fc = 1000;

t = linspace(0,1/freqSep,duration);
y = zeros(duration,85);

for i=1:85
    y(:,i) = 0.5*cos(2*pi*(fc+symbolsOut(i)*freqSep).*t);
end
y=transpose(y(:));
%% Pulse shaping
%pulse shaping may be used, however it was not implemented for experimental
%transmission, for use uncomment lines 134-142 and comment line 130
% T = 1/freqSep;
% shapeLength = round(T/8*Fs);
% A = 0.5*(1-cos(8*pi.*t(1:shapeLength)/T));
% yA = y;
% for i=1:85
%     yA(1:shapeLength,i) = y(1:shapeLength,i).*transpose(A(1,1:shapeLength));
%     yA(duration-shapeLength+1:duration,i) = y(duration-shapeLength+1:duration,i).*transpose(A(1,shapeLength:-1:1));
% end
% y=transpose(yA(:));
%% Save output
%audio output for radiostation is generated
%60s length is due to protocol detection in WSJT-X
out = zeros(1,Fs*60);
out(1,1:length(y)) = y(:);
% out = awgn(out,-15);
audiowrite('out.wav',out,Fs)