clear;clc;
%this code is adjusted to decode 'LDPC_local.wav'
%% Load input
[y,Fs] = audioread('LDPC_local.wav');
duration = 1920;
deltaF = 6.25;
Fc = 1000;
M = 8;
%% Demodulate and detect symbols
detected = zeros(1,94);
yFilter=bandpass(y,[Fc-(M/2+1)*deltaF Fc+(M/2+1)*deltaF],Fs); %filtration

for i=1:90
    [yPower,fPower]=pspectrum(yFilter((i-1)*duration+1:i*duration),Fs,...
        'FrequencyResolution',20,'FrequencyLimits',[Fc-M*deltaF Fc+M*deltaF]);
    detected(i) = round( (fPower(find(yPower==max(yPower)))-Fc)/deltaF );
end
clear yPower fPower
%% Remove synchronization
desync = [detected(12:40) detected(44:72)]; %matched manually
%% Decode Gray code
deGray = gray2bin(desync,'fsk',8);
%% Convert from symbols to bits
A = [];
for i = 1:58
    A = [A (dec2bin(deGray(i),3))];
end
bits = zeros(1,174);
for i=1:length(A)
    bits(i) = str2double(A(i));
end
clear A
%% Decode
%LDPC code (174,91) is specified by parity matrix which is stored separately
load('LDPCparity.mat');
ldpcDecoder = comm.LDPCDecoder(sparse(transpose(parity)));
decodedBits = mod(1+(ldpcDecoder(transpose(bits))),2);
%% 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