%function pb=multi_rayleigh_vita_84encoded_all(Eb_to_Nj_in_dB,Eb_to_No_in_dB,BPH,number_of_states,Q)
%VITERBI      This program simulates the Viterbi sequnce detection of the
%             differential frequency hopping system with diversity of level 2
%            
%             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
%             這個程序是正確的,帶(8,4)糾錯碼,編碼后在后面只補上L個0,即符號流總長度是2N+L而不是2(N+L),
%             而multi_awgn_vita_encode.m是帶(8,4)糾錯碼但符號流總長度是2(N+L),這兩個程序都對
%           
%             這個程序是在帶(8,4)糾錯碼的multi_awgn_vita_encodev.m基礎上加了交織,并且是基于比特的分組(block)交織
%             

%********系統參數********%
% Eb_to_Nj_in_dB=11;
Eb_to_Nj_in_dB=[9,12,15];
Eb_to_No_in_dB=13.35;
BPH=2;
number_of_states=16;
Q=4; % 干擾音的個數
N=1000; % 每次符號流長度
times=500; % 重復做500次

PB_all=zeros(length(Eb_to_Nj_in_dB),1);
vita_symbol_err=zeros(length(Eb_to_Nj_in_dB),times);% 維特比譯碼后, (8,4) 譯碼之前2^BPH=4進制符號錯誤個數計數器
num_of_err=zeros(length(Eb_to_Nj_in_dB),times);
%************************%
WTbarpp=waitbar(0,'Outer SNRs Loop:Please wait...');
for pp=1:length(Eb_to_Nj_in_dB)


fanout=2^BPH; % DFH的扇出系數
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級是其狀態位,與網格圖中的狀態一一對應(注意并不是與跳頻頻點一一對應)

% pb=zeros(1,times);
WTbar=waitbar(0,'SNR inside loop:please wait...');
for rep=1:times

% source=[randint(1,N,fanout),zeros(1,L)]; % 信息源:(注意不是10進制的,而是fanout進制的)隨機符號流,最后補上L個0符號,使移位寄存器的狀態清零
source=randint(1,N,fanout);
%***********************信源部分************************%

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

%******************************************************%


% *******************差錯控制編碼部分(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)擴展漢明碼的生成矩陣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)擴展漢明碼的一致校驗矩陣H
% trt=syndtable(H); % (8,4,4)擴展漢明碼的伴隨式譯碼表
dsource_coded=encode(dsource,8,4,'linear',G)';
% source_coded0=zeros(1,2*N);
% *****************************************************************%

% *********************交織部分(塊交織)(基于比特的交織)****************%
% 僅適用于BPH=2且N=1000時
interleave_outcome=zeros(1,2*BPH*N);
A1=zeros(50,80);
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
%************************************************%


% % *******************隨機交織************************%
% [interleave_outcome,alpha]=interleave(dsource_coded); % dsource_coded 長度為(n/k)*N*BPH=2*N*BPH
% 
% % ***************************************************%
% 編碼和交織后將二進制序列轉換成十進制序列,準備輸入到G函數進行映射
if (BPH~=1)
    for i=1:2*N
        register=interleave_outcome((i-1)*BPH+1:i*BPH);
        source_coded_ba0(i)=change2deci(register,2); % 其實是把BPH個二進制比特轉化為一個2^BPH進制的符號,這里BPH=2,就是把2個二進制比特轉化為1個四進制符號(是對的)
    end
else
    source_coded_ba0=interleave_outcome(1:2*N);% source_coded_ba0 是行矢量,source_coded_ba0是編碼輸出十進制符號序列,長度為2N個符號(補2L個0符號之前)
end

source_coded=[source_coded_ba0,zeros(1,L)];% G函數映射前補上L個0符號,長度為2N+L個fanout進制符號

% ******************* G 函數實現部分 ************************ %

% 先定義三個關鍵矩陣"nextstates" "output" "input"
nextstate=zeros(number_of_states,fanout); % nextstate矩陣:行代表網格圖中的各狀態(一一對應),列與輸入移位寄存器的信息符號一一對應,
                                          % 矩陣中存儲的內容是與當前狀態和輸入符號對應的下一狀態號(即存儲網格圖的狀態轉移規則)  
output=zeros(number_of_states,fanout);    % output矩陣:行代表網格圖中的各狀態(一一對應),列與輸入移位寄存器的信息符號一一對應,
                                          % 矩陣中存儲的內容是與當前狀態和輸入符號對應的網格圖分支轉移輸出(分支轉移輸出是跳頻頻率號)
input=zeros(number_of_states,number_of_states);
number_of_out=number_of_states*fanout;% 跳頻頻點數Nt
for i=0:number_of_states-1
    for j=0:fanout-1
        [next_state,out_put]=G_func1(i,j,L,fanout);
        nextstate(i+1,j+1)=next_state;
        output(i+1,j+1)=out_put;
        input(i+1,next_state+1)=j;
    end
end 
% ********************************************************* %


% ********************維特比譯碼部分**********************%

depth_of_trellis=length(source_coded);%*******************% depth_of_trellis長度為2N+L
Eb=1;
Es=Eb*BPH*(1/2);
diversi=2; % 2重分集
Ec=Es/diversi; % 分集后每個分集碼片的能量
% Ej0=(Eb*number_of_out*Q)/(Eb_to_Nj);% 每個多音干擾的能量Ej0
Ej0=(Es*number_of_out)/(BPH*Q*Eb_to_Nj); % 每跳時間內每個多音干擾的能量Ej0
sgma=sqrt(Eb/(2*Eb_to_No));% AWGN的均方根
demod_input=zeros(number_of_out,depth_of_trellis);
f=zeros(1,depth_of_trellis);
% rc=zeros(1,number_of_out);
% rs=zeros(1,number_of_out);
D=0; % D 記錄網格圖的當前狀態,這里初始狀態是0狀態

% % *******************信道和非相干解調部分:加多音干擾和噪聲,然后非相干解調****************** %
% for i=1:depth_of_trellis % i表示網格圖的時間走勢
%     f(i)=output(D+1,source_coded(i)+1); % f(i)是i時刻的分支轉移輸出,即i時刻的跳頻頻率號,頻率號范圍是[0,number_of_out-1]而不是[1,number_of_out]
%     thyta=2*pi*rand;% 干擾音與跳頻信號的相對相位
% %   J=randint(1,Q,number_of_out); % J 矩陣中存放Q個干擾音所在的頻率號,干擾音所在頻率號范圍也是[0,number_of_out-1]而不是[1,number_of_out]
%     J=gen_multijammer(Q,number_of_out);% J 矩陣中存放Q個干擾音所在的頻率號,干擾音所在頻率號范圍也是[0,number_of_out-1]而不是[1,number_of_out]
%     for j=0:number_of_out-1
%         if (j==f(i))
%             rc(j+1)=raylrnd(sqrt(Es))+sgma*randn;% 改了,原來是sqrt(Es)+sgma*randn,信號幅度服從瑞利分布,瑞利分布的參數是sqrt(Es)*222222222222222222222222*
%             rs(j+1)=sgma*randn;
%         else
%             rc(j+1)=sgma*randn;
%             rs(j+1)=sgma*randn;
%         end
%     end
%     jam_rayleigh=raylrnd(sqrt(Ej0));% 干擾音的幅度服從瑞利分布,瑞利分布的參數是sqrt(Ej0)*222222222222222222222222222222222222222222222222222222222222222*
%     for k=1:Q
%         for j=0:number_of_out-1
%             if (j==J(k))
%                 rc(j+1)=rc(j+1)+jam_rayleigh*cos(thyta);%改了,原來是rc(j+1)+sqrt(Ej0)*cos(thyta)*2222222222222222222222222222222222222222222222222222222222*
%                 rs(j+1)=rs(j+1)+jam_rayleigh*sin(thyta);%改了,原來是rs(j+1)+sqrt(Ej0)*sin(thyta)*2222222222222222222222222222222222222222222222222222222222*
%             end
%         end
%     end
%     for j=0:number_of_out-1
%         demod_input(j+1,i)=sqrt(rc(j+1)^2+rs(j+1)^2);
%     end
%     D=nextstate(D+1,source_coded(i)+1); 
% end
%             
% % **************** End of the Rayleigh Fading Channel and Noncoherent Demodulation Modular *************** %

% *******************信道和非相干解調部分:加多音干擾和噪聲,然后非相干解調,得到解調輸出即維特比譯碼器的輸入****************** %
for i=1:depth_of_trellis % i表示網格圖的時間走勢
    f(i)=output(D+1,source_coded(i)+1); % f(i)是i時刻的分支轉移輸出,即i時刻輸出的跳頻頻率號