%function pb=multi_rayleigh_vita_213encoded_all(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,K)=(1,2,3)的卷積糾錯(cuò)碼,編碼后在后面只補(bǔ)上L個(gè)0,即符號(hào)流總長度是2N+L而不是2(N+L),
%                        
%             這個(gè)程序是帶(1,2,3)的卷積糾錯(cuò)碼并且加了基于比特的分組(block)交織
%             

%********系統(tǒng)參數(shù)********%
Eb_to_Nj_in_dB=[0,3,6];
Eb_to_No_in_dB=13.35;
BPH=2;
number_of_states=16;
Q=4; % 干擾音的個(gè)數(shù)
N=1000; % 每次符號(hào)流長度
times=10; % 重復(fù)做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進(jìn)制符號(hào)錯(cuò)誤個(gè)數(shù)計(jì)數(shù)器
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的扇出系數(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級(jí)是其狀態(tài)位,與網(wǎng)格圖中的狀態(tài)一一對(duì)應(yīng)(注意并不是與跳頻頻點(diǎn)一一對(duì)應(yīng))

% num_of_err=zeros(1,times);
% 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進(jìn)制的,而是fanout進(jìn)制的)隨機(jī)符號(hào)流,最后補(bǔ)上L個(gè)0符號(hào),使移位寄存器的狀態(tài)清零
source=randint(1,N,fanout);
%***********************信源部分************************%

% 將信息符號(hào)流轉(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);
    end
else
    dsource=source(1:N);%*****
end

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

% *******************差錯(cuò)控制編碼部分(n,k,K)=(1,2,3)卷積碼**********************%
KK=3; %卷積碼的約束長度K=3
trel=poly2trellis(KK,[7,5]); % 卷積碼的網(wǎng)格圖trel
dsource_coded=convenc(dsource,trel,0);
% *****************************************************************%


% *********************交織部分(塊交織)(基于比特的交織)****************%
% 僅適用于BPH=2且N=1000時(shí)
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
%************************************************%


% % *******************隨機(jī)交織************************%
% [interleave_outcome,alpha]=interleave(dsource_coded); % dsource_coded 長度為(n/k)*N*BPH=2*N*BPH
% 
% % ***************************************************%
% 編碼和交織后將二進(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);
    end
else
    source_coded_ba0=interleave_outcome(1:2*N);% source_coded_ba0 是行矢量,source_coded_ba0是編碼輸出十進(jìn)制符號(hào)序列,長度為2N個(gè)符號(hào)(補(bǔ)2L個(gè)0符號(hào)之前)
end

source_coded=[source_coded_ba0,zeros(1,L)];% G函數(shù)映射前補(bǔ)上L個(gè)0符號(hào)

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

% 先定義三個(gè)關(guān)鍵矩陣"nextstates" "output" "input"
nextstate=zeros(number_of_states,fanout); % nextstate矩陣:行代表網(wǎng)格圖中的各狀態(tài)(一一對(duì)應(yīng)),列與輸入移位寄存器的信息符號(hào)一一對(duì)應(yīng),
                                          % 矩陣中存儲(chǔ)的內(nèi)容是與當(dāng)前狀態(tài)和輸入符號(hào)對(duì)應(yīng)的下一狀態(tài)號(hào)(即存儲(chǔ)網(wǎng)格圖的狀態(tài)轉(zhuǎn)移規(guī)則)  
output=zeros(number_of_states,fanout);    % output矩陣:行代表網(wǎng)格圖中的各狀態(tài)(一一對(duì)應(yīng)),列與輸入移位寄存器的信息符號(hào)一一對(duì)應(yīng),
                                          % 矩陣中存儲(chǔ)的內(nèi)容是與當(dāng)前狀態(tài)和輸入符號(hào)對(duì)應(yīng)的網(wǎng)格圖分支轉(zhuǎn)移輸出(分支轉(zhuǎn)移輸出是跳頻頻率號(hào))
input=zeros(number_of_states,number_of_states);
number_of_out=number_of_states*fanout;% 跳頻頻點(diǎn)數(shù)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);%*******************
Eb=1;
Es=Eb*BPH*(1/2);
% Ej0=(Eb*number_of_out*Q)/(Eb_to_Nj);% 每個(gè)多音干擾的能量Ej0
Ej0=(Es*number_of_out)/(BPH*Q*Eb_to_Nj); % 每跳時(shí)間內(nèi)每個(gè)多音干擾的能量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 記錄網(wǎng)格圖的當(dāng)前狀態(tài),這里初始狀態(tài)是0狀態(tài)

% *******************信道和非相干解調(diào)部分:加多音干擾和噪聲,然后非相干解調(diào)****************** %
for i=1:depth_of_trellis % i表示網(wǎng)格圖的時(shí)間走勢(shì)
    f(i)=output(D+1,source_coded(i)+1); % f(i)是i時(shí)刻的分支轉(zhuǎn)移輸出,即i時(shí)刻的跳頻頻率號(hào),頻率號(hào)范圍是[0,number_of_out-1]而不是[1,number_of_out]
    theta=2*pi*rand;% 干擾音與跳頻信號(hào)的相對(duì)相位
%   J=randint(1,Q,number_of_out); % J 矩陣中存放Q個(gè)干擾音所在的頻率號(hào),干擾音所在頻率號(hào)范圍也是[0,number_of_out-1]而不是[1,number_of_out]
    J=gen_multijammer(Q,number_of_out);% J 矩陣中存放Q個(gè)干擾音所在的頻率號(hào),干擾音所在頻率號(hào)范圍也是[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,信號(hào)幅度服從瑞利分布,瑞利分布的參數(shù)是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));% 干擾音的幅度服從瑞利分布,瑞利分布的參數(shù)是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(theta);%改了,原來是rc(j+1)+sqrt(Ej0)*cos(thyta)*2222222222222222222222222222222222222222222222222222222222*
                rs(j+1)=rs(j+1)+jam_rayleigh*sin(theta);%改了,原來是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 Modulation Modular *************** %

% ***********************維特比譯碼模塊********************************* %
state_metric=zeros(number_of_states,2);
survivor_state=zeros(number_of_states,depth_of_trellis+1);
for i=1:depth_of_trellis-L
    flag=zeros(1,number_of_states);
    if i<=L+1
        step=2^((L+1-i)*BPH);
    else
        step=1;
    end
    for j=0:step:number_of_states-1
        for m=0:fanout-1
            branch_metric=demod_input(output(j+1,m+1)+1,i);
            if((state_metric(nextstate(j+1,m+1)+1,2)<state_metric(j+1,1)...
                    +branch_metric)|flag(nextstate(j+1,m+1)+1)==0)
                state_metric(nextstate(j+1,m+1)+1,2)=state_metric(j+1,1)+branch_metric;
                survivor_state(nextstate(j+1,m+1)+1,i+1)=j;
                flag(nextstate(j+1,m+1)+1)=1;
            end
        end
    end
    state_metric=state_metric(:,2:-1:1);
end
for i=depth_of_trellis-L+1:depth_of_trellis
    flag=zeros(1,number_of_states);
    last_stop=number_of_states/(2^((i-depth_of_trellis+L-1)*BPH));
    for j=0:last_stop-1
        branch_metric=demod_input(output(j+1,m+1)+1,i);
        if((state_metric(nextstate(j+1,1)+1,2)<state_metric(j+1,1)...
            +branch_metric)|flag(nextstate(j+1,1)+1)==0)
            state_metric(nextstate(j+1,1)+1,2)=state_metric(j+1,1)+branch_metric;
            survivor_state(nextstate(j+1,1)+1,i+1)=j;
            flag(nextstate(j+1,1)+1)=1;
        end
    end
    state_metric=state_metric(:,2:-1:1);
end
state_sequence=zeros(1,depth_of_trellis+1);
for i=1:depth_of_trellis
    state_sequence(1,depth_of_trellis-i+1)=survivor_state((state_sequence(1,depth_of_trellis+2-i)...
        +1),depth_of_trellis-i+2);
end
decoder_output=zeros(1,BPH*(depth_of_trellis-L));% length(decoder_output)=BPH*2N
for i=1:depth_of_trellis-L  % 在G函數(shù)維特比譯碼的時(shí)候就已經(jīng)考慮到將原來補(bǔ)上的L個(gè)0符號(hào)給去掉了,G函數(shù)維特比譯碼后的輸出是2N個(gè)十進(jìn)制符號(hào)(轉(zhuǎn)換成BPH*2N個(gè)比特)
    dec_output_deci=input(state_sequence(1,i)+1,state_sequence(1,i+1)+1);% 輸出的十進(jìn)制符號(hào)寄存器dec_output_deci(其實(shí)是fanout進(jìn)制的符號(hào))
                                                                         % 因?yàn)橄旅嬉痪鋎ec_output_bin=deci2change(dec_output_deci,BPH,2)顯然是
                                                                         % 把dec_output_deci轉(zhuǎn)化為 BPH 位二進(jìn)制的,所以dec_output_deci
                                                                         % 應(yīng)該是2^BPH=fanout進(jìn)制而不是十進(jìn)制的    
    vitdec_output(i)=dec_output_deci; % 存放維特比譯碼輸出的fanout進(jìn)制的符號(hào)序列
    if(BPH~=1)
        dec_output_bin=deci2change(dec_output_deci,BPH,2);
        decoder_output((i-1)*BPH+1:i*BPH)=dec_output_bin; %decoder_output是維特比譯碼輸出的二進(jìn)制序列,length(decoder_output)=BPH*2N
    else
        decoder_output(1,i)=dec_output_deci;
    end
end
% ****************** End of the Viterbi Do

% % ***********************隨機(jī)解交織*********************%
% deint_output_bin=deinterleave(decoder_output,alpha); % BPH*2N個(gè)比特進(jìn)行解交織,2N個(gè)符號(hào)即BPH*2N個(gè)比特
% 
% % ****************************************************%

% *********************解交織部分(塊交織)(基于比特的解交織)****************%
% 僅適用于BPH=2且N=1000時(shí)
deint_output_bin=zeros(1,2*BPH*N);% 解交織輸出deint_output_bin
B1=zeros(80,50);
for i=1:2*BPH*N
    B1(i)=decoder_output(i);
end
B=B1.';
for i=1:2*BPH*N
    deint_output_bin(i)=B(i);% deint_output_bin為解交織后輸出的二進(jìn)制比特序列,行矢量(也就是糾錯(cuò)碼譯碼的輸入序列)
end
% ************************************************************************%


% ***********(n,k,K)=(1,2,3)卷積碼譯碼部分***********%
deco_output=vitdec(deint_output_bin,trel,15,'trunc','hard');

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

for i=1:2*N % % 計(jì)算維特比譯碼后,(8,4)譯碼之前的fanout=4進(jìn)制符號(hào)錯(cuò)誤數(shù),共有2N個(gè)符號(hào)
    if(vitdec_output(i)~=source_coded(i))
        vita_symbol_err(pp,rep)=vita_symbol_err(pp,rep)+1;% 這時(shí)vita_symbol_err是維特比譯碼后,(8,4)譯碼之前的4進(jìn)制符號(hào)錯(cuò)誤數(shù)
    end
end

for i=1:N*BPH
    if(dsource(i)~=deco_output(i))
       num_of_err(pp,rep)=num_of_err(pp,rep)+1;
    end
end

waitbar(rep/times,WTbar)
end % 與最外層"多少遍rep" for循環(huán)對(duì)應(yīng)的end
close(WTbar)
waitbar(pp/length(Eb_to_Nj_in_dB),WTbarpp)
end % 對(duì)應(yīng)于最最外面pp信干比循環(huán)的end
close(WTbarpp)



PB_all=sum(num_of_err,2)/(N*times*BPH) % N 是一遍的符號(hào)流長度,總共進(jìn)行times遍,共N*times個(gè)符號(hào)

symerr_num_biterr=zeros(2*length(Eb_to_Nj_in_dB),times);% symerr_num_biterr矩陣中奇數(shù)行是某一信干比下維特比譯碼后,RS譯碼之前的2^BPH=4進(jìn)制符號(hào)錯(cuò)誤數(shù)
                                                        % 偶數(shù)行是對(duì)應(yīng)的最終RS譯碼后的2^M=16進(jìn)制符號(hào)錯(cuò)誤數(shù)
for i=1:length(Eb_to_Nj_in_dB)
    symerr_num_biterr(2*(i-1)+1,:)=vita_symbol_err(i,:);
    symerr_num_biterr(2*i,:)=num_of_err(i,:);
end