multi_rayleigh_vita_213encoded_all_wdiv2.m

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

%function pb=multi_rayleigh_vita_213encoded_all(Eb_to_Nj_in_dB,Eb_to_No_in_dB,BPH,number_of_states,Q)
%VITERBI      This program simulates the Viterbi sequnce detection of the
%             differential frequency hopping system with diversity of level 2
%
%             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
%
%             这个程序是正确的,带(n,k,K)=(1,2,3)的卷积纠错码,编码后在后面只补上L个0,即符号流总长度是2N+L而不是2(N+L),
%                        
%             这个程序是带(1,2,3)的卷积纠错码并且加了基于比特的分组(block)交织
%             

%********系统参数********%
Eb_to_Nj_in_dB=[12,15,18];
Eb_to_No_in_dB=13.35;
BPH=2;
number_of_states=16;
Q=4; % 干扰音的个数
N=1000; % 每次符号流长度
times=500; % 重复做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进制符号错误个数计数器
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的扇出系数
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级是其状态位,与网格图中的状态一一对应(注意并不是与跳频频点一一对应)

% 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进制的,而是fanout进制的)随机符号流,最后补上L个0符号,使移位寄存器的状态清零
source=randint(1,N,fanout);
%***********************信源部分************************%

% 将信息符号流转化成二进制信息比特流
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

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

% *******************差错控制编码部分(n,k,K)=(1,2,3)卷积码**********************%
KK=3; %卷积码的约束长度K=3
trel=poly2trellis(KK,[7,5]); % 卷积码的网格图trel
dsource_coded=convenc(dsource,trel,0);
% *****************************************************************%


% *********************交织部分(块交织)(基于比特的交织)****************%
% 仅适用于BPH=2且N=1000时
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
%************************************************%


% % *******************随机交织************************%
% [interleave_outcome,alpha]=interleave(dsource_coded); % dsource_coded 长度为(n/k)*N*BPH=2*N*BPH
% 
% % ***************************************************%
% 编码和交织后将二进制序列转换成十进制序列,准备输入到G函数进行映射
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是编码输出十进制符号序列,长度为2N个符号(补2L个0符号之前)
end

source_coded=[source_coded_ba0,zeros(1,L)];% G函数映射前补上L个0符号

% ******************* 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_coded);%*******************
Eb=1;
Es=Eb*BPH*(1/2);
diversi=2; % 2重分集
Ec=Es/diversi; % 分集后每个分集码片的能量
% 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的均方根
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表示网格图的时间走势
%     f(i)=output(D+1,source_coded(i)+1); % f(i)是i时刻的分支转移输出,即i时刻的跳频频率号,频率号范围是[0,number_of_out-1]而不是[1,number_of_out]
%     theta=2*pi*rand;% 干扰音与跳频信号的相对相位
% %   J=randint(1,Q,number_of_out); % J 矩阵中存放Q个干扰音所在的频率号,干扰音所在频率号范围也是[0,number_of_out-1]而不是[1,number_of_out]
%     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)=raylrnd(sqrt(Es))+sgma*randn;% 改了,原来是sqrt(Es)+sgma*randn,信号幅度服从瑞利分布,瑞利分布的参数是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));% 干扰音的幅度服从瑞利分布,瑞利分布的参数是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 *************** %

% *******************信道和非相干解调部分:加多音干扰和噪声,然后非相干解调,得到解调输出即维特比译码器的输入****************** %
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