%     D=nextstate(D+1,source_coded(i)+1); 
% end
% % ********************************************************************************************************** %            

% *******************信道和非相干解調部分:加多音干擾和噪聲,然后非相干解調,得到解調輸出即維特比譯碼器的輸入****************** %
for i=1:depth_of_trellis % i表示網格圖的時間走勢
    f(i)=output(D+1,source_coded(i)+1); % f(i)是i時刻的分支轉移輸出,即i時刻輸出的跳頻頻率號
    
    rc1=zeros(1,number_of_out);% 每次搞完一跳后都把rc數組清零,準備存放下一跳的相關解調器輸出數據,第1個分集chip的rc
    rc2=zeros(1,number_of_out);% 第2個分集chip的rc,當然也可以把rc定義為分集重數diversi那么多行,number_of_out那么多列的一個矩陣，可能會更簡練一些
    rs1=zeros(1,number_of_out);% 每次搞完一跳后都把rs數組清零,準備存放下一跳的相關解調器輸出數據,第1個分集chip的rs
    rs2=zeros(1,number_of_out);% 第2個分集chip的rs
    
    % 這里我沒有必要把整個頻率號都生成出來再重復diversity那么多次再交織,我只需把加干擾的那部分程序運行diversity次其結果與前面采用交織是一樣的
    % 給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; % 很顯然,分集后這里應該是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的平方律解調輸出判決統計量(能量)
    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; % 很顯然,分集后這里應該是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的平方律解調輸出判決統計量(能量)
    end
    
    for j=0:number_of_out-1
        demod_input(j+1,i)=R_chip1(j+1)+R_chip2(j+1);
    end
    D=nextstate(D+1,source_coded(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=zeros(1,BPH*(depth_of_trellis-L));% length(decoder_output)=BPH*2N
for i=1:depth_of_trellis-L  % 在G函數維特比譯碼的時候就已經考慮到將原來補上的L個0符號給去掉了,G函數維特比譯碼后的輸出是2N個十進制符號(轉換成BPH*2N個比特)
    dec_output_deci=input(state_sequence(1,i)+1,state_sequence(1,i+1)+1);% 輸出的十進制符號寄存器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; %decoder_output是維特比譯碼輸出的二進制序列,length(decoder_output)=BPH*2N
    else
        decoder_output(1,i)=dec_output_deci;
    end
end

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


% % ***********************隨機解交織*********************%
% deint_output_bin=deinterleave(decoder_output,alpha); % BPH*2N個比特進行解交織,2N個符號即BPH*2N個比特
% 
% % ****************************************************%


% *********************解交織部分(塊交織)(基于比特的解交織)****************%
% 僅適用于BPH=2且N=1000時
deint_output_bin=zeros(1,2*BPH*N);% 解交織輸出deint_output_bin
B1=zeros(200,200);
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為解交織后輸出的二進制比特序列,行矢量(也就是糾錯碼譯碼的輸入序列)
end
% ************************************************************************%


% % ***********(8,4)譯碼部分***********%
% deco_input=zeros(BPH*2*N,1); % 列矢量
% 
% deco_input(1:BPH*2*N)=deint_output_bin(1:BPH*2*N);% 把行矢量轉化為列矢量(其實完全沒有這個必要,decode函數完全可以接受行矢量,效果與列矢量是一樣的)
% deco_output=decode(deco_input,8,4,'linear',G)';% (8,4)譯碼器的譯碼輸出(長為BPH*N)個比特,行矢量)
% 
% % ***********************************%

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

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

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循環對應的end
close(WTbar)
waitbar(pp/length(Eb_to_Nj_in_dB),WTbarpp)
end % 對應于最最外面pp信干比循環的end
close(WTbarpp)

PB_all=sum(num_of_err,2)/(N*times*BPH) % N 是一遍的符號流長度,總共進行times遍,共N*times個符號