multi_rayleigh_vita.asv

短波信道抗多音干扰的性能分析及其仿真

%function pb=multi_rayleigh_vita(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


%********系统参数********%
Eb_to_Nj_in_dB=[1];
Eb_to_No_in_dB=15;
BPH=2;
number_of_states=16;
Q=8; % 干扰音的个数
%************************%

N=100000; % 每次符号流长度
times=50; % 重复做50次
fanout=2^BPH; % DFH的扇出系数

vita_symbol_err=zeros(length(Eb_to_Nj_in_dB),times);% 维特比译码后符号错误计数器
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...');
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级是其状态位,与网格图中的状态一一对应(注意并不是与跳频频点一一对应)

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

source=[randint(1,N,fanout),zeros(1,L)]; % 信息源:随机符号流,最后补上L个0符号,使移位寄存器的状态清零

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

% 将信息符号流转化成二进制信息比特流,不包括最后补上的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 函数实现部分 ************************ %

% 先定义三个关键矩阵"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);
% Es=1;
% Eb=Es/BPH;
Eb=1;
Es=BPH*Eb;
% 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的均方根
sgma=sqrt(Eb/(2*Eb_to_No));% AWGN的均方根
% 原来这里是thyta=2*pi*rand;% 干扰音与跳频信号的相对相位
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表示网格图的时间走势
    rc=zeros(1,number_of_out);
    rs=zeros(1,number_of_out);
%     theta_s=2*pi*rand;
%     theta_j=2*pi*rand;% 干扰音与跳频信号的相对相位,每一跳信号与干扰之间的相对相位都应该不同
    f(i)=output(D+1,source(i)+1); % f(i)是i时刻的分支转移输出,即i时刻的跳频频率号
    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))
            theta_s=2*pi*rand;% //random phase of the signal tone
            rc(j+1)=raylrnd(sqrt(Es/2))*cos(theta_s)+sgma*randn;
            rs(j+1)=raylrnd(sqrt(Es/2))*sin(theta_s)+sgma*randn;
        else
            rc(j+1)=sgma*randn;
            rs(j+1)=sgma*randn;
        end
    end
    jam_rayleigh=raylrnd(sqrt(Ej0/2));% 干扰音的幅度服从瑞利分布,瑞利分布的参数是sqrt(Ej0)
    for k=1:Q
        theta_j=2*pi*rand;% //random phase of each jamming tone
        for j=0:number_of_out-1
            if (j==J(k))
                rc(j+1)=rc(j+1)+jam_rayleigh*cos(theta_j);%改了,原来是rc(j+1)+sqrt(Ej0)*cos(thyta)
                rs(j+1)=rs(j+1)+jam_rayleigh*sin(theta_j);%改了,原来是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);
%         demod_input(j+1,i)=rc(j+1)^2+rs(j+1)^2;
    end
    D=nextstate(D+1,source(i)+1);
end
% **************** End of Channel and Noncoherent demodulation modular ***************** %

% ****************** Viterbi decoding 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_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
% ******************** End of Viterbi Decoding Modular ****************** %

for i=1:N  %计算维特比译码之后的符号错误数
    if (decoder_output_symbol(i)~=source(i))
        vita_symbol_err(pp,rep)=vita_symbol_err(pp,rep)+1;
    end
end

for i=1:N*BPH % 计算维特比译码之后的比特错误数
    if(dsource(i)~=decoder_output(i))
       num_of_err(pp,rep)=num_of_err(pp,rep)+1;
   end
end

waitbar(rep/times,WTbar)
end % 与最外层"多少遍" for rep=1:times循环对应的end
close(WTbar)
waitbar(pp/length(Eb_to_Nj_in_dB),WTbarpp)
end % 与最最外层信干比循环 pp=1:length(Eb_to_Nj_in_dB)对应的end
close(WTbarpp)




Pb=sum(num_of_err,2)/(N*times*BPH) % N 是一遍的符号流长度,总共进行times遍,共N*times个符号
symerr_num_biterr=zeros(2*length(Eb_to_Nj_in_dB),times);% symerr_num_biterr矩阵中奇数行是某一信干比下维特比译码后的2^BPH=4进制符号错误数
                                                        % 偶数行是对应的比特错误数
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