%function pb=multi_awgn_sym100sym(Eb_to_Nj_in_dB,Eb_to_No_in_dB,BPH,number_of_states,Q)
%Symbol-by-Symbol      This procedure simulates the symbol-by-symbol 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 band


%********系統(tǒng)參數(shù)********%
Eb_to_Nj_in_dB=[6,7,8,9];
Eb_to_No_in_dB=13.35;
BPH=2;
number_of_states=16;
Q=4; % 干擾音的個(gè)數(shù)
%************************%

N=1000; % 每次符號(hào)流長度
times=50; % 重復(fù)做500次
fanout=2^BPH; % DFH的扇出系數(shù)

num_of_err=zeros(length(Eb_to_Nj_in_dB),times);
Ps=zeros(length(Eb_to_Nj_in_dB),1);
Pb=zeros(length(Eb_to_Nj_in_dB),1);

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

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))

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

% ********************信源模塊****************************** %

source=[randint(1,N,fanout),zeros(1,L)]; % 信息源:隨機(jī)符號(hào)流,最后補(bǔ)上L個(gè)0符號(hào),使移位寄存器的狀態(tài)清零

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

% ******************* 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);% depth_of_trellis長度為N+L,其中L=2
Eb=1;
Es=BPH*Eb;
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);% demod_input矩陣為number_of_out=64行,depth_of_trellis=N+L列,行表示頻率號(hào),列表示時(shí)間走勢
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í)間走勢,depth_of_trellis長度為N+L,其中L=2
    f(i)=output(D+1,source(i)+1); % f(i)是i時(shí)刻的分支轉(zhuǎn)移輸出,即i時(shí)刻的跳頻頻率號(hào),f(i)是跳頻頻率號(hào)序列
    thyta=2*pi*rand;% 干擾音與跳頻信號(hào)的相對(duì)相位,每一跳信號(hào)音干擾音的相位都應(yīng)該不同
%   J=randint(1,Q,number_of_out); % J 矩陣中存放Q個(gè)干擾音所在的頻率號(hào)
    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)=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);% 非相干解調(diào)輸出存放于demod_input矩陣
    end
    D=nextstate(D+1,source(i)+1);
end
% *********************上面是信道和非相干解調(diào)模塊*********************************%

% ********************逐符號(hào)(Symbol-by-Symbol)譯碼模塊************************* %
[Max_energy freq_number]=max(demod_input,[],1);% 找出每跳時(shí)間內(nèi)能量最大的那個(gè)頻率號(hào),存放于freq_number矢量中

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

for i=1:N  %計(jì)算逐符號(hào)譯碼之后的符號(hào)錯(cuò)誤率,發(fā)送頻率序列與接收頻率序列比較
    if ((freq_number(i)-1)~=f(i))
        num_of_err(pp,rep)=num_of_err(pp,rep)+1;
    end
end

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

Ps=sum(num_of_err,2)/(N*times) % 計(jì)算符號(hào)錯(cuò)誤率,N 是一遍的符號(hào)流長度,總共進(jìn)行times遍,共N*times個(gè)符號(hào)
Pb=Ps*fanout/(2*(fanout-1))