clear;clc;close all;
%this code is adjusted to decode 'JT9_local.wav'
%% Load input
[y,Fs] = audioread('JT9_local.wav');
duration = 6912;
deltaF = 1.736;
Fc = 1000+1*deltaF; %manual carrier frequency compensation
%% Demodulate and detect symbols
numSym = 87; %number of detected symbols expanded
detected = zeros(1,numSym);
yFilter=bandpass(y,[Fc-deltaF Fc+9*deltaF],Fs);

for i=1:numSym
    [yPower,fPower]=pspectrum(yFilter((i-1)*duration+1:i*duration),Fs,...
        'FrequencyResolution',10,'FrequencyLimits',[Fc-deltaF/2 Fc+8.5*deltaF]);
    detected(i) = round((fPower(find(yPower==max(yPower)))-Fc)/deltaF);
end
detected(1:85) = detected(2:86); %manual correction of delayed start of transmission
% figure;stem(detected);title(Fc);
clear yPower fPower
%% Remove sync symbols
syncVectorS = '1100100001000001000000100000000010100000000000000011001000010000010000001000000000101';
syncVector = zeros(1,85);
for i=1:85
    syncVector(i) = str2double(syncVectorS(i));
end

j=1;
desync = zeros(1,69);
for i=1:85
    if syncVector(i) == 0
        desync(j) = detected(i);
        j = j+1;
    end
end
clear syncVectorS syncVector
%% Decode Gray code
deGray = gray2bin(desync-1,'fsk',8);
%% Convert from symbols to bits
A = [];
for i = 1:69
    A = [A (dec2bin(deGray(i),3))];
end
bits = zeros(1,207);
for i=1:length(A)
    bits(i) = str2double(A(i));
end
clear A
%% Deinterleaving
%"mapping" contains new indexes for deinterleaving "bits"
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
     symbol(i) = bits(mapping(i)+1);
 end

clear n j0 k
%% Fano decode initialization
%this Fano algorithm implementation was created by Jonathon Cheah
%source: https://www.mathworks.com/matlabcentral/fileexchange/39221-wspr-txrx
%initial conditions and memory allocation adjusted for our requirements

% set the initial variable configuration
MAXBITS=103;
nstate=zeros(MAXBITS-1,1);
gamma=zeros(MAXBITS-1,1);
metrics=zeros(4,MAXBITS-1);
tm=zeros(2,MAXBITS-1);

% configure FEC sequence
npoly1=2^32-221228207; % 0xF2D05351, defined in JT9 specification
npoly2=2^32-463389625; % 0xE4613C47, defined in JT9 specification

% configure the Fano decoder
partab=textread('partab.dat','%u');
mettab=textread('mettab.dat','%d');
mettab=reshape(mettab,256,2);
%% Fano decode
symbol = double(symbol);
nsym=206;
nbits=103;
ndelta=72;
maxcycles=20000;
ntail=nbits-31;
i4a=0;
i4b=0;
for np=1:nbits
    j=2*np;
    i4a=-symbol(j-1);
    i4b=-symbol(j);
    if (i4a<0)
        i4a=i4a+256;
    end
    if (i4b<0)
        i4b=i4b+256;
    end
    metrics(1,np) = mettab(i4a+1,1) + mettab(i4b+1,1);
    metrics(2,np) = mettab(i4a+1,1) + mettab(i4b+1,2);
    metrics(3,np) = mettab(i4a+1,2) + mettab(i4b+1,1);
    metrics(4,np) = mettab(i4a+1,2) + mettab(i4b+1,2);
end
np=1;
nstate(np)=0;

% Compute and sort fano branch metrics
n=bitand(nstate(np),npoly1);
n=bitxor(n,bitshift(n,-16));
m=bitand(bitxor(n,bitshift(n,-8)),255);
lsym=partab(m+1);
n=bitand(nstate(np),npoly2);
n=bitxor(n,bitshift(n,-16));
m=bitand(bitxor(n,bitshift(n,-8)),255);
lsym=lsym+lsym+partab(m+1);
m0=metrics(lsym+1,np);
m1=metrics(bitxor(3,lsym)+1,np);
% tm(1,np)=m0 if 0-branch is better
%         =m1 if 1-branch is better
if m0>m1
    tm(1,np) = m0;
    tm(2,np) = m1;
else
    tm(1,np) = m1;
    tm(2,np) = m0;
    nstate(np) = mod(nstate(np)+1,2^32);
end

ii(np)=0;
gamma(np)=0;
nt=0;

% Start fano with the best branch
for i=1:nbits*maxcycles-1
    ngamma=gamma(np) + tm(ii(np)+1,np);
    if ngamma>nt
        if gamma(np)<(nt+ndelta)
            nt=nt + ndelta * floor((ngamma-nt)/ndelta);
        end
        
        gamma(np+1)=ngamma;
        nstate(np+1)=mod(bitshift(nstate(np),1),2^32);
        np=np+1;
        
        % fano decoding done.
        if np == nbits
            break
        end
        
        % fano processing.
        n=mod(bitand(nstate(np),npoly1),2^32);
        n=mod(bitxor(n,bitshift(n,-16)),2^32);
        lsym=partab(bitand(bitxor(n,bitshift(n,-8)),255)+1);
        n=mod(bitand(nstate(np),npoly2),2^32);
        n=mod(bitxor(n,bitshift(n,-16)),2^32);
        lsym=lsym+lsym+partab(bitand(bitxor(n,bitshift(n,-8)),255)+1);
        
        % fano at the tail
        if np >ntail+1
            tm(1,np)=metrics(lsym+1,np);
        else
            m0=metrics(lsym+1,np);
            m1=metrics(bitxor(3,lsym)+1,np);
            if m0>m1
                tm(1,np)=m0;
                tm(2,np)=m1;
            else
                tm(1,np)=m1;
                tm(2,np)=m0;
                nstate(np)=mod(nstate(np) + 1,2^32);
            end
        end
        % Start at the best branch
        ii(np)=0  ;
        continue
    end
    npp=0;
    
    while(true)
        noback=0;
        if np==1
            noback=1;
        end
        if np>1
            if gamma(np-1)<nt
                noback=1;
            end
        end
        if(noback)
            nt=nt-ndelta;
            if ii(np)~=0
                ii(np)=0;
                nstate(np)=mod(bitxor(nstate(np),1),32);
            end
            break
        end
        
        % Back up search
        np=np-1;
        if(np<ntail+1 && ii(np)~=1)
            %Search next best branch
            ii(np)=ii(np)+1;
            nstate(np)=mod(bitxor(nstate(np),1),2^32);
            break
        end
    end
end
metric=gamma(np);

decodedBits = mod(nstate(:,1),2);
clear nstate gamma metric maxcycles nt ii np m0 m1 nsym ntail tm n nbits
clear ndelta ngamma npoly1 npoly2 m lsym j i4a i4b MAXBITS metrics mettab
clear partab
%% Unpacking of decoded bit sequence
N1 = bin2dec(join(string(decodedBits(1:28))));
N2 = bin2dec(join(string(decodedBits(29:56))));
N3 = bin2dec(join(string(decodedBits(57:72))));

N3 = (N3 - 32768); %substraction of 32768=2^15 is to remove plain text flag bit
N2 = bitshift(N2,-1);
N1 = bitshift(N1,-1);

msgCoded = zeros(1,13);
msgCoded(13) = mod(N3,42);
msgCoded(12) = mod((N3-msgCoded(13))/42,42);
msgCoded(11) = ((N3-msgCoded(13))/42-msgCoded(12))/42;

msgCoded(10) = mod(N2,42);
msgCoded(9) = mod((N2-msgCoded(10))/42,42);
msgCoded(8) = mod(((N2-msgCoded(10))/42-msgCoded(9))/42,42);
msgCoded(7) = mod((((N2-msgCoded(10))/42-msgCoded(9))/42-msgCoded(8))/42,42);
msgCoded(6) = mod(((((N2-msgCoded(10))/42-msgCoded(9))/42-msgCoded(8))/42-msgCoded(7))/42,42);

msgCoded(5) = mod(N1,42);
msgCoded(4) = mod((N1-msgCoded(5))/42,42);
msgCoded(3) = mod(((N1-msgCoded(5))/42-msgCoded(4))/42,42);
msgCoded(2) = mod((((N1-msgCoded(5))/42-msgCoded(4))/42-msgCoded(3))/42,42);
msgCoded(1) = mod(((((N1-msgCoded(5))/42-msgCoded(4))/42-msgCoded(3))/42-msgCoded(2))/42,42);
clear N1 N2 N3
%% Source decoding
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',' ','+','-','.','/','?'];
msgLength = 13;
for i=1:msgLength
    msgDecoded(i) = alphabet(msgCoded(i)+1);
end
disp(msgDecoded); %should be equal to one of messages below
% msg = 'A1B2 .?C3D4E5'; %message1
% msg = 'OK1FLZ ABCDEF'; %message2
% msg = 'OK1FLZ JN79IA'; %message3
clear msgLength alphabet i