clear all; close all; clc

%% Input parametrs
robustens_mode = 2; % Режим помехоустойчивости 1; 2; 3; 4;
spectrum_mode = 5; % Размещение спектра 4.5кГц - 0; 5кГц - 1; 9кГц - 2; 10кГц - 3; 18кГц - 4; 20кГц - 5;

level_protection = 0; % Уровень защиты для QAM16 - 0; QAM4 - 1;

QAM_Pmax = 2; % 4-QAM Pmax = 1; 16-QAM Pmax = 2;

%% Initialization parametrs

% Колличество ячеек
N_mux_list = [167, 190, 359, 405, 754, 846; ...
130, 150, 282, 322, 588, 662; ...
NaN, NaN, NaN, 288, NaN, 607; ...
NaN, NaN, NaN, 152, NaN, 332];

% Кодовые скорости канала
Rate = cell(2,2,1);
Rate{1,1} = [1 3];
Rate{1,2} = [2 3];
Rate{2,1} = [1 2];
Rate{2,2} = [NaN];

% Схема выкалывания канала
puncpat = cell(2,2,1);
puncpat{1,1} = [1 1 1 0 0 0];
puncpat{1,2} = [1 1 0 0 0 0 1 0 0 0 0 0];
puncpat{2,1} = [1 1 0 0 0 0];
puncpat{2,2} = [NaN];

% Перфорирование остаточных 36 бит
puncpat_Rp = cell(12,1);
puncpat_Rp{1,1} = [1 1 0 0 0 0 1 1 0 0 0 0 1 1 0 0 0 0 1 1 0 0 0 0 1 1 0 0 0 0 1 1 0 0 0 0];
puncpat_Rp{2,1} = [1 1 1 0 0 0 1 1 0 0 0 0 1 1 0 0 0 0 1 1 0 0 0 0 1 1 0 0 0 0 1 1 0 0 0 0];
puncpat_Rp{3,1} = [1 1 1 0 0 0 1 1 0 0 0 0 1 1 0 0 0 0 1 1 1 0 0 0 1 1 0 0 0 0 1 1 0 0 0 0];
puncpat_Rp{4,1} = [1 1 1 0 0 0 1 1 1 0 0 0 1 1 0 0 0 0 1 1 1 0 0 0 1 1 0 0 0 0 1 1 0 0 0 0];
puncpat_Rp{5,1} = [1 1 1 0 0 0 1 1 1 0 0 0 1 1 0 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 0 0 0 0];
puncpat_Rp{6,1} = [1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 0 0 0 0];
puncpat_Rp{7,1} = [1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0];
puncpat_Rp{8,1} = [1 1 1 1 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0];
puncpat_Rp{9,1} = [1 1 1 1 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 1 0 0 1 1 1 0 0 0 1 1 1 0 0 0];
puncpat_Rp{10,1} = [1 1 1 1 0 0 1 1 1 1 0 0 1 1 1 0 0 0 1 1 1 1 0 0 1 1 1 0 0 0 1 1 1 0 0 0];
puncpat_Rp{11,1} = [1 1 1 1 0 0 1 1 1 1 0 0 1 1 1 0 0 0 1 1 1 1 0 0 1 1 1 0 0 0 1 1 1 1 0 0];
puncpat_Rp{12,1} = [1 1 1 1 0 0 1 1 1 1 0 0 1 1 1 1 0 0 1 1 1 1 0 0 1 1 1 0 0 0 1 1 1 1 0 0];

Rp = zeros(1,QAM_Pmax);
for k = 1:QAM_Pmax
Rp(k) = (2*N_mux_list(robustens_mode,spectrum_mode+1)-12) - Rate{level_protection+1,k}(2)...
*floor((2*N_mux_list(robustens_mode,spectrum_mode+1) - 12)./Rate{level_protection+1,k}(2)) + 1;
end

% Кодовая скорость материнского кода
mother_code = 1/6;

% Задание таблицы переходов
trellis = poly2trellis(7,[133, 171, 145, 133,171,145]);

Mbase = [2 4];
M=2.^Mbase(QAM_Pmax);

a = [1/sqrt(2), 1/sqrt(10)]; % Коэффициент нормализации QAM

% Модуляция QAM16 и QAM4
if QAM_Pmax == 1
SMQAM = [1+1*i, -1+1*i, -1-1*i, 1-1*i].*a(QAM_Pmax);

elseif QAM_Pmax == 2
SMQAM = [3+3*i, 3-1*i, 3+1*i, 3-3*i,...
-1+3*i, -1-1*i, -1+1*i, -1-3*i,...
1+3*i, 1-1*i, 1+1*i, 1-3*i,...
-3+3*i, -3-1*i, -3+1*i, -3-3*i].*a(QAM_Pmax);
end

% Расчет колличества бит
L_2 = 0;
for k = 1:QAM_Pmax
Mp_2(k) = Rate{level_protection+1,k}(1)*floor((2*N_mux_list(robustens_mode,spectrum_mode+1
)...
- 12)./Rate{level_protection+1,k}(2)); % Формула 77
end
L_2 = sum(Mp_2);

% Задаем схему выкалывания
punc = cell(QAM_Pmax,1);
for n = 1:QAM_Pmax
Num_puncpat = ((Mp_2(n)+ 6)/mother_code)./length(puncpat{level_protection+1,n});
m = [];
for k = 1:ceil(Num_puncpat)
m = [m ,puncpat{level_protection+1,n}];
end
punc{n,:} = [m(1:Num_puncpat*length(puncpat{level_protection+1,n}) - 36), puncpat_Rp{Rp(n),1}];
end

% QAM модулятор
H_QAM = modem.genqammod('Constellation', SMQAM, 'InputType', 'bit');

% Пременные
x_in = 2*N_mux_list(robustens_mode,spectrum_mode+1);
pBuff_out = zeros(x_in, 1);
decoded_out = [];
Input_Data_coder = [];
N_simbol = 1000;

%% Программа кодера декодера
tblen = 35;
quant = 3;
AmplQuant = 0.4;
delta = 2*AmplQuant/(2^quant-2);
for EbNo = 1:20
for l = 1:N_simbol
%% Формирование последовательности бит
Input_Data=randi([0 1],L_2,1);
Input_Data_coder = [Input_Data_coder; Input_Data];

%% Сверточный кодер
init_state = 0;
if QAM_Pmax == 2
[In_Data_1] = convenc([Input_Data(Mp_2(1)+1:Mp_2(1)+Mp_2(2)); [0 0 0 0 0 0]'],trellis,punc{2, 1}, init_state);
[In_Data_0] = convenc([Input_Data(1:Mp_2(1)); [0 0 0 0 0 0]'],trellis,punc{1, 1}, init_state);
else
[In_Data_0] = convenc([Input_Data(1:Mp_2(1)); [0 0 0 0 0 0]'],trellis,punc{1, 1}, init_state);
end

%% Мультиплексор
if QAM_Pmax == 2
pBuff_out = [In_Data_0; In_Data_1];
else
pBuff_out = In_Data_0;
end

% Отображение на созвездие QAM
Mod_qam=modulate(H_QAM,pBuff_out);

%% Канал передачи AWGN
k=log2(M);
Qam_rx=awgn(Mod_qam,EbNo+10*log10(k*1/2),'measured','dB');

%% Расчет ОСШ канала
EsNo = EbNo*k/2;
SNR = sum(real(Qam_rx).^2 + imag(Qam_rx).^2)/(EsNo*M);

%% Обратное отображение на созвездие QAM
HD_QAM = modem.genqamdemod('Constellation', SMQAM, 'OutputType', 'bit','DecisionType', 'approximate llr','NoiseVariance', SNR);
Demod_QAM=demodulate(HD_QAM,Qam_rx);

%% Демультиплексор
if QAM_Pmax == 2
In_Data_0 = Demod_QAM(1:x_in);
In_Data_1 = Demod_QAM(x_in+1:2*x_in);
else
In_Data_0 = Demod_QAM(1:x_in);
end

%% Декодирование
if QAM_Pmax == 2
[x,qcode_1] = quantiz(In_Data_1(:, 1),(-AmplQuant:delta:AmplQuant), (2^quant-1:-1:0));
[x,qcode_0] = quantiz(In_Data_0(:, 1),(-AmplQuant:delta:AmplQuant), (2^quant-1:-1:0));
[decoded_1] = vitdec(qcode_1,trellis,tblen,'trunc','soft', quant, punc{2, 1});
[decoded_0] = vitdec(qcode_0,trellis,tblen,'trunc','soft', quant, punc{1, 1});
else
[x,qcode_0] = quantiz(In_Data_0(:, 1),(-AmplQuant:delta:AmplQuant), (2^quant-1:-1:0));
[decoded_0] = vitdec(qcode_0,trellis,tblen,'trunc','soft', quant, punc{1, 1});
end

%% Мультиплексор
if QAM_Pmax == 2
decoded = [decoded_0(1:Mp_2(1)), decoded_1(1:Mp_2(2))];
else
decoded = decoded_0(1:Mp_2(1));
end

decoded_out = [decoded_out, decoded];
end
%% Подсчет количества ошибок
[num,BER]=biterr(Input_Data_coder(1:end),decoded_out(1:end)');
BERcod(EbNo) = BER;
fprintf('ОСШ: %f BER: %f \n', EbNo, BER);
Input_Data_coder = [];
decoded_out = [];
end

semilogy(BERcod);
grid on;