%function pb=multi_awgn_vita(Eb_to_Nj_in_dB,Eb_to_No_in_dB,BPH,number_of_states,Q)
%VITERBI      This procedure simulates the Viterbi sequnce decoding of the
%             differential frequency hopping system
%
%             Eb_to_Nj_in_dB is the signal-to-jamming ratio given in dB
%
%             Eb_to_No_in_dB is the signal-to-noise ratio given in dB
%
%             BPH is the number of bits transmitted by one hop
%
%             number_of_states is the number of states in the DFH
%             trellis,corresponding to the right L stages of the DFH encoding 
%             shift register 
%
%             Q is the number of jamming tones in the DFH bandwidth
%
%             這個(gè)程序是正確的,帶(N,K,D)=(15,7,9)的 RS 糾錯(cuò)碼,編碼后在后面只補(bǔ)上L個(gè)0,即符號流總長度是2N+L而不是2(N+L),
%             
%           
%             
%             

%********系統(tǒng)參數(shù)********%
% Eb_to_Nj_in_dB=11;
Eb_to_Nj_in_dB=[14];
Eb_to_No_in_dB=15;
BPH=2;
number_of_states=16;
Q=1; % 干擾音的個(gè)數(shù)

nRS=15; % (15,7)RS code
kRS=7;
M=4; % 伽羅華域GF(2^m)中m=4,即GF(16)域
% N=7*3000; % 每次符號流長度,N個(gè)16進(jìn)制的信息符號,即4N個(gè)比特,注意N一定要取7的倍數(shù)才行
N=kRS*3000;
times=20; % 重復(fù)做500次

PB_all=zeros(length(Eb_to_Nj_in_dB),1);
PS_all=zeros(length(Eb_to_Nj_in_dB),1);

vita_symbol_err=zeros(length(Eb_to_Nj_in_dB),times);% 維特比譯碼后, RS 譯碼之前2^BPH=4進(jìn)制符號錯(cuò)誤個(gè)數(shù)計(jì)數(shù)器
num_of_err=zeros(length(Eb_to_Nj_in_dB),times);% 最終 RS 碼譯碼輸出的2^M=16進(jìn)制符號錯(cuò)誤個(gè)數(shù)計(jì)數(shù)器


% 4位二進(jìn)制序列與GF(16)域中1個(gè)16進(jìn)制GF符號對應(yīng)轉(zhuǎn)換關(guān)系表
binary_table=[0 0 0 0;
              0 0 0 1;
              0 0 1 0;
              0 0 1 1;
              0 1 0 0;
              0 1 0 1;
              0 1 1 0;
              0 1 1 1;
              1 0 0 0;
              1 0 0 1;
              1 0 1 0;
              1 0 1 1;
              1 1 0 0;
              1 1 0 1;
              1 1 1 0;
              1 1 1 1];
% 2位四進(jìn)制序列與GF(16)域中1個(gè)16進(jìn)制GF符號對應(yīng)轉(zhuǎn)換關(guān)系表
quarter_table=[0 0;
               0 1;
               0 2;
               0 3;
               1 0;
               1 1;
               1 2;
               1 3;
               2 0;
               2 1;
               2 2;
               2 3;
               3 0;
               3 1;
               3 2;
               3 3];
decimal_table=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15]';
GF16_table=gf(decimal_table,M);% GF(16)域中GF符號表


%************************%
WTbarpp=waitbar(0,'Outer SNRs Loop:Please wait');
for pp=1:length(Eb_to_Nj_in_dB)


fanout=2^BPH; % DFH的扇出系數(shù)
Eb_to_Nj=10^(Eb_to_Nj_in_dB(pp)/10); % 比值形式的Eb/Nj
Eb_to_No=10^(Eb_to_No_in_dB/10); % 比值形式的Eb/No
L=floor(log(number_of_states)/log(fanout)); % 編碼移位寄存器的長度為L+1,最右邊L級是其狀態(tài)位,與網(wǎng)格圖中的狀態(tài)一一對應(yīng)(注意并不是與跳頻頻點(diǎn)一一對應(yīng))

% pb=zeros(1,times);
WTbar=waitbar(0,'SNR inside loop:please wait for a moment...');
for rep=1:times
% source=[randint(1,N,fanout),zeros(1,L)]; % 信息源:(注意不是10進(jìn)制的,而是fanout進(jìn)制的)隨機(jī)符號流,最后補(bǔ)上L個(gè)0符號,使移位寄存器的狀態(tài)清零

%***************信源部分(信源是長度為N的16進(jìn)制符號流)*********************%

source_sequence=randint(1,N,2^M);
source_matrix=reshape(source_sequence,kRS,N/kRS); % 特別注意reshape函數(shù)是從矩陣A中按列的順序取出元素,返回一個(gè)M*N的矩陣,但是注意是把元素按先填第1列,然后填第2列……的順序填到新矩陣中去的，不是按行填寫
source=source_matrix';               % reshape之后得到的是一個(gè)7行N/7(kRS=7)列的矩陣,轉(zhuǎn)置后得到的是一個(gè)N/7行7列的矩陣source,矩陣source的第一行存放的是source_sequence
                                     % 序列的第1到7個(gè)符號,第二行存放的是source_sequence序列的第8到14個(gè)符號,依此類推    
%*************************************************************************%

% %***********************信源部分************************%
% 
% % 將信息符號流轉(zhuǎn)化成二進(jìn)制信息比特流
% dsource=zeros(1,N*BPH); 
% if(BPH~=1)
%     for i=1:N
%         dsource((i-1)*BPH+1:i*BPH)=deci2change(source(i),BPH,2); %%%%%%%%%%%%其實(shí)這里是把 2^BPH 進(jìn)制的符號source(i)轉(zhuǎn)換成BPH個(gè)二進(jìn)制比特,前面產(chǎn)生的符號流source=randint(1,N,fanout)
%                                                                  %產(chǎn)生的符號流原來本身也就是fanout=2^BPH進(jìn)制的,這正好相吻合,所以這里并不是把十進(jìn)制的符號source(i)轉(zhuǎn)化為
%                                                                  %BPH位二進(jìn)制比特,而是把2^BPH進(jìn)制的符號source(i)(顯然source=randint(1,N,fanout)表明source(i)
%                                                                  %的取值范圍被限定在0到fanout-1之間)轉(zhuǎn)化為BPH位二進(jìn)制比特,
%     end
% else
%     dsource=source(1:N);%*****
% end
% 
% %******************************************************%


% % *******************差錯(cuò)控制編碼部分(8,4)碼**********************%
% % G=[1 0 0 0 1 0 1 1;
% %    0 1 0 0 1 1 1 0;
% %    0 0 1 0 1 1 0 1;
% %    0 0 0 1 0 1 1 1];
% G=[1 0 0 0 1 0 1 1;0 1 0 0 1 1 0 1;0 0 1 0 0 1 1 1;0 0 0 1 1 1 1 0]; % (8,4,4)擴(kuò)展?jié)h明碼的生成矩陣G
% % H=[1 1 1 1 1 1 1 1;1 1 0 1 1 0 0 0;0 1 1 1 0 1 0 0;1 0 1 1 0 0 1 0]; % (8,4,4)擴(kuò)展?jié)h明碼的一致校驗(yàn)矩陣H
% % trt=syndtable(H); % (8,4,4)擴(kuò)展?jié)h明碼的伴隨式譯碼表
% dsource_coded=encode(dsource,8,4,'linear',G)';
% % source_coded0=zeros(1,2*N);
% % *****************************************************************%

%****************差錯(cuò)控制編碼部分(N,K,Dmin)=(15,7,9)的 RS 碼*****************%

msg_GF=gf(source,M); %把信息符號轉(zhuǎn)化為GF域的符號
% source_coded_GF=rsenc(msg_GF,15,7); % RS編碼后,source_coded_GF是一個(gè)N/7行15列的矩陣,矩陣每一行是一個(gè) RS 碼的碼字
source_coded_GF=rsenc(msg_GF,nRS,kRS);

%****************************************************************************%
% % *********************交織部分(塊交織)(基于比特的交織)****************%
% % 僅適用于BPH=2且N=1000時(shí)
% interleave_outcome=zeros(1,2*BPH*N);
% A1=zeros(200,200);
% for i=1:2*BPH*N
%     A1(i)=dsource_coded(i);
% end
% A=A1.';
% for i=1:2*BPH*N
%     interleave_outcome(i)=A(i);
% end
% %************************************************%

% % *******************隨機(jī)交織************************%
% [interleave_outcome,alpha]=interleave(dsource_coded); % dsource_coded 長度為(n/k)*N*BPH=2*N*BPH
% 
% % ***************************************************%

%**************交織部分(基于16進(jìn)制GF(16)域符號(即RS碼元)的交織)***************%
% codedsequence=zeros(1,15*N/7);
codedsequence=zeros(1,nRS*N/kRS);
codedsequence_GF=gf(codedsequence,M);
% for i=1:15*N/7
for i=1:nRS*N/kRS
    codedsequence_GF(i)=source_coded_GF(i); % 把RS編碼后的碼字矩陣轉(zhuǎn)化為一串碼字序列,按列讀出,相當(dāng)于nRS*N/kRS的分組交織
end
%*****************************************************************************%
% % % 編碼和交織后將二進(jìn)制序列轉(zhuǎn)換成十進(jìn)制序列,準(zhǔn)備輸入到G函數(shù)進(jìn)行映射
% % if (BPH~=1)
% %     for i=1:2*N
% %         register=interleave_outcome((i-1)*BPH+1:i*BPH);
% %         source_coded_ba0(i)=change2deci(register,2); % 其實(shí)是把BPH個(gè)二進(jìn)制比特轉(zhuǎn)化為一個(gè)2^BPH進(jìn)制的符號,這里BPH=2,就是把2個(gè)二進(jìn)制比特轉(zhuǎn)化為1個(gè)四進(jìn)制符號(是對的)
% %     end
% % else
% %     source_coded_ba0=interleave_outcome(1:2*N);% source_coded_ba0 是行矢量,source_coded_ba0是編碼輸出十進(jìn)制符號序列,長度為2N個(gè)符號(補(bǔ)2L個(gè)0符號之前)
% % end
% 
% source_coded=[source_coded_ba0,zeros(1,L)];% G函數(shù)映射前補(bǔ)上L個(gè)0符號




% 編碼和交織后將GF(16)域的16進(jìn)制符號轉(zhuǎn)化為普通域的fanout=4進(jìn)制的符號,準(zhǔn)備輸入到G函數(shù)進(jìn)行映射
% codedsequence_quart=zeros(1,(M/BPH)*15*N/7);
codedsequence_quart=zeros(1,(M/BPH)*nRS*N/kRS);
% for i=1:15*N/7
for i=1:nRS*N/kRS
    for j=0:2^M-1
        if (codedsequence_GF(i)==j)
            codedsequence_quart(2*(i-1)+1:2*i)=quarter_table(j+1,:);
        end
    end
end   % 把1500個(gè)編碼后的16進(jìn)制 RS 碼元都轉(zhuǎn)化為3000個(gè)4進(jìn)制符號,準(zhǔn)備輸入到G函數(shù)進(jìn)行映射(通過測試,結(jié)果正確,可以轉(zhuǎn)化)
source_coded=[codedsequence_quart,zeros(1,L)]; % G函數(shù)映射前補(bǔ)上L=2個(gè)0符號,總共3002個(gè)4進(jìn)制符號進(jìn)行G函數(shù)映射

% ******************* G 函數(shù)實(shí)現(xiàn)部分 ************************ %

% 先定義三個(gè)關(guān)鍵矩陣"nextstates" "output" "input"
nextstate=zeros(number_of_states,fanout); % nextstate矩陣:行代表網(wǎng)格圖中的各狀態(tài)(一一對應(yīng)),列與輸入移位寄存器的信息符號一一對應(yīng),
                                          % 矩陣中存儲的內(nèi)容是與當(dāng)前狀態(tài)和輸入符號對應(yīng)的下一狀態(tài)號(即存儲網(wǎng)格圖的狀態(tài)轉(zhuǎn)移規(guī)則)  
output=zeros(number_of_states,fanout);    % output矩陣:行代表網(wǎng)格圖中的各狀態(tài)(一一對應(yīng)),列與輸入移位寄存器的信息符號一一對應(yīng),
                                          % 矩陣中存儲的內(nèi)容是與當(dāng)前狀態(tài)和輸入符號對應(yīng)的網(wǎng)格圖分支轉(zhuǎn)移輸出(分支轉(zhuǎn)移輸出是跳頻頻率號)
input=zeros(number_of_states,number_of_states);
number_of_out=number_of_states*fanout;% 跳頻頻點(diǎn)數(shù)Nt