%function pb=multi_awgn_vita(Eb_to_Nj_in_dB,Eb_to_No_in_dB,BPH,number_of_states,Q)
%VITERBI      This script m-file simulates the Viterbi sequnce decoding of the
%             differential frequency hopping system with diversity of 3 levels
%
%             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


%********系統(tǒng)參數(shù)********%
Eb_to_Nj_in_dB=9;
Eb_to_No_in_dB=13.35;
BPH=2;
number_of_states=16;
Q=8; % 干擾音的個數(shù)
%************************%




N=1000; % 每次符號流長度
times=1; % 重復(fù)做500次
fanout=2^BPH; % DFH的扇出系數(shù)
Eb_to_Nj=10^(Eb_to_Nj_in_dB/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))

vita_symbol_err=zeros(1,times);% 維特比譯碼后符號錯誤計數(shù)器
num_of_err=zeros(1,times);
% pb=zeros(1,times);
WTbar=waitbar(0,'Please wait...');% waitbar等待條
for rep=1:times

source=[randint(1,N,fanout),zeros(1,L)]; % 信息源:隨機(jī)符號流,最后補(bǔ)上L個0符號,使移位寄存器的狀態(tài)清零,很顯然source長度為N+L個fanout進(jìn)制的符號

%***********************信源部分************************%

% 將信息符號流轉(zhuǎn)化成二進(jìn)制信息比特流,不包括最后補(bǔ)上的L個0符號(dsource只是在最后計算誤碼率的時候用)
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

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



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

% 先定義三個關(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
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);
% Es=1;
% Eb=Es/BPH;
Eb=1;
Es=BPH*Eb;
diversi=3; % 3重分集
Ec=Es/diversi; % 分集后每個分集碼片的能量
% Ej0=(Eb*number_of_out*Q)/(Eb_to_Nj);% 每個多音干擾的能量Ej0
Ej0=(Es*number_of_out)/(BPH*Q*Eb_to_Nj); % 每跳時間內(nèi)每個多音干擾的能量Ej0               “這個能量是否要隨著每跳信號的能量而變還待定”
sgma=sqrt(Eb/(2*Eb_to_No));% AWGN的均方根sgma
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)格圖的時間走勢
%     f(i)=output(D+1,source(i)+1); % f(i)是i時刻的分支轉(zhuǎn)移輸出,即i時刻的跳頻頻率號
%     thyta=2*pi*rand;% 干擾音與跳頻信號的相對相位
% %   J=randint(1,Q,number_of_out); % J 矩陣中存放Q個干擾音所在的頻率號
%     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)=sqrt(Es)+sgma*randn;
%             rs(j+1)=sgma*randn;
%         else
%             rc(j+1)=sgma*randn;
%             rs(j+1)=sgma*randn;
%         end
%     end
%     for k=1:Q
%         for j=0:number_of_out-1
%             if (j==J(k))
%                 rc(j+1)=rc(j+1)+sqrt(Ej0)*cos(thyta);
%                 rs(j+1)=rs(j+1)+sqrt(Ej0)*sin(thyta);
%             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(i)+1);
% end
% % ***************************************************************************************** %

% *******************信道和非相干解調(diào)部分:加多音干擾和噪聲,然后非相干解調(diào),得到解調(diào)輸出即維特比譯碼器的輸入****************** %
for i=1:depth_of_trellis % i表示網(wǎng)格圖的時間走勢
    f(i)=output(D+1,source(i)+1); % f(i)是i時刻的分支轉(zhuǎn)移輸出,即i時刻輸出的跳頻頻率號
    
    rc1=zeros(1,number_of_out);% 每次搞完一跳后都把rc數(shù)組清零,準(zhǔn)備存放下一跳的相關(guān)解調(diào)器輸出數(shù)據(jù),第1個分集chip的rc
    rc2=zeros(1,number_of_out);% 第2個分集chip的rc,當(dāng)然也可以把rc定義為分集重數(shù)diversi那么多行,number_of_out那么多列的一個矩陣，可能會更簡練一些
    rc3=zeros(1,number_of_out);% 第3個分集chip的rc
    rs1=zeros(1,number_of_out);% 每次搞完一跳后都把rs數(shù)組清零,準(zhǔn)備存放下一跳的相關(guān)解調(diào)器輸出數(shù)據(jù),第1個分集chip的rs
    rs2=zeros(1,number_of_out);% 第2個分集chip的rs
    rs3=zeros(1,number_of_out);% 第3個分集chip的rs
    
    % 這里我沒有必要把整個頻率號都生成出來再重復(fù)diversity那么多次再交織,我只需把加干擾的那部分程序運(yùn)行diversity次其結(jié)果與前面采用交織是一樣的
    % 給chip1加干擾和噪聲
    theta=2*pi*rand;% 干擾音與跳頻信號的相對相位
%   J=randint(1,Q,number_of_out); % J 矩陣中存放Q個干擾音所在的頻率號
    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))
            rc1(j+1)=sqrt(Ec)+sgma*randn; % 很顯然,分集后這里應(yīng)該是Ec,而不是原來的Es
            rs1(j+1)=sgma*randn;
        else
            rc1(j+1)=sgma*randn;
            rs1(j+1)=sgma*randn;
        end
    end
    for k=1:Q
        for j=0:number_of_out-1
            if (j==J(k))
                rc1(j+1)=rc1(j+1)+sqrt(Ej0)*cos(theta);
                rs1(j+1)=rs1(j+1)+sqrt(Ej0)*sin(theta);
            end
        end
    end
    for j=0:number_of_out-1
        R_chip1(j+1)=rc1(j+1)^2+rs1(j+1)^2;% 第1個chip的平方律解調(diào)輸出判決統(tǒng)計量(能量)
    end
    
    % 給chip2加干擾和噪聲
    theta=2*pi*rand;% 干擾音與跳頻信號的相對相位
%   J=randint(1,Q,number_of_out); % J 矩陣中存放Q個干擾音所在的頻率號
    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))
            rc2(j+1)=sqrt(Ec)+sgma*randn; % 很顯然,分集后這里應(yīng)該是Ec,而不是原來的Es
            rs2(j+1)=sgma*randn;
        else
            rc2(j+1)=sgma*randn;
            rs2(j+1)=sgma*randn;
        end
    end
    for k=1:Q
        for j=0:number_of_out-1
            if (j==J(k))
                rc2(j+1)=rc2(j+1)+sqrt(Ej0)*cos(theta);
                rs2(j+1)=rs2(j+1)+sqrt(Ej0)*sin(theta);
            end
        end
    end
    for j=0:number_of_out-1
        R_chip2(j+1)=rc2(j+1)^2+rs2(j+1)^2;% 第2個chip的平方律解調(diào)輸出判決統(tǒng)計量(能量)
    end
    
    % 給chip3加干擾和噪聲
    theta=2*pi*rand;% 干擾音與跳頻信號的相對相位
%   J=randint(1,Q,number_of_out); % J 矩陣中存放Q個干擾音所在的頻率號
    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))
            rc3(j+1)=sqrt(Ec)+sgma*randn; % 很顯然,分集后這里應(yīng)該是Ec,而不是原來的Es
            rs3(j+1)=sgma*randn;
        else
            rc3(j+1)=sgma*randn;
            rs3(j+1)=sgma*randn;
        end
    end
    for k=1:Q
        for j=0:number_of_out-1
            if (j==J(k))
                rc3(j+1)=rc3(j+1)+sqrt(Ej0)*cos(theta);
                rs3(j+1)=rs3(j+1)+sqrt(Ej0)*sin(theta);
            end
        end
    end
    for j=0:number_of_out-1
        R_chip3(j+1)=rc3(j+1)^2+rs3(j+1)^2;% 第2個chip的平方律解調(diào)輸出判決統(tǒng)計量(能量)
    end
    
    for j=0:number_of_out-1
        demod_input(j+1,i)=R_chip1(j+1)+R_chip2(j+1)+R_chip3(j+1);
    end
    D=nextstate(D+1,source(i)+1);
end
% *************************************************************************************************************************** %




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_symbol=zeros(1,depth_of_trellis-L);% length(decoder_output_symbol)=N,維特比譯碼后的輸出符號流
decoder_output=zeros(1,BPH*(depth_of_trellis-L));
for i=1:depth_of_trellis-L
    dec_output_deci=input(state_sequence(1,i)+1,state_sequence(1,i+1)+1);
    decoder_output_symbol(i)=dec_output_deci;
    if(BPH~=1)
        dec_output_bin=deci2change(dec_output_deci,BPH,2);
        decoder_output((i-1)*BPH+1:i*BPH)=dec_output_bin;
    else
        decoder_output(1,i)=dec_output_deci;
    end
end

for i=1:N  %計算維特比譯碼之后的符號錯誤率
    if (decoder_output_symbol(i)~=source(i))
        vita_symbol_err(rep)=vita_symbol_err(rep)+1;
    end
end

for i=1:N*BPH
    if(dsource(i)~=decoder_output(i))
       num_of_err(rep)=num_of_err(rep)+1;
   end
end
waitbar(rep/times,WTbar)
end % 與最外層"多少遍" for循環(huán)對應(yīng)的end
close(WTbar)

Pb=sum(num_of_err)/(N*times*BPH) % N 是一遍的符號流長度,總共進(jìn)行times遍,共N*times個符號
symerr_num_biterr=zeros(2,times);
symerr_num_biterr=[vita_symbol_err;num_of_err];