clear;clc;close all
%this code is adjusted to decode 'LDPC_WebSDR_France.wav'
%% Load input
[y,Fs] = audioread('LDPC_WebSDR_France.wav');
y = resample(y,3,2); %resampling is necesary as WebSDR records at Fs=8kHz
Fs = 12000;
duration = 1920;
deltaF = 6.25;
M = 8;
Fc1 = 1000; %expected carrier frequency
%% Carrier frequncy shift compensation
%due to transmision and remote receiver properties, shift in frequency may
%occur, we find value of shifted carrier frequency by finding bandwidth with
%max power, this calculation is the most time occuping in the code
A = 0;
for i=1:30
    Fc = Fc1+i*deltaF/2;
    [yPower,fPower]=pspectrum(y,Fs,'FrequencyResolution',20,...
        'FrequencyLimits',[Fc-deltaF Fc+(M+1)*deltaF]);
    if A < sum(yPower)
        A = sum(yPower);
        index = i;
    end
end
Fc = Fc1+index*deltaF/2+deltaF;
clear A Fc1 index
%% Demodulate and detect symbols
numSym = floor(length(y)/duration);
detected = zeros(1,numSym);
yFilter = bandpass(y,[Fc-(M+1)*deltaF Fc+(M+1)*deltaF],Fs); %filtration

for i=1:numSym
    [yPower,fPower]=pspectrum(yFilter((i-1)*duration+1:i*duration),Fs,...
        'FrequencyResolution',20,'FrequencyLimits',[Fc-deltaF Fc+M*deltaF]);
    detected(i) = round( (fPower(find(yPower==max(yPower)))-Fc)/deltaF );
end

% figure;stem((detected));grid on;title('detected');
% figure;pspectrum(yFilter,Fs);
%% Remove synchronization
desync = [detected(34:62) detected(70:98)]; %4.wav
% figure;stem(desync);grid on;title('desync');
%% 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