multi_rayleigh_sym100sym.m

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

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


%********系统参数********%
Eb_to_Nj_in_dB=[21,22];
Eb_to_No_in_dB=13.35;
BPH=2;
number_of_states=16;
Q=4; % 干扰音的个数
%************************%

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

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级是其状态位,与网格图中的状态一一对应(注意并不是与跳频频点一一对应)

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

% ********************信源模块****************************** %

source=[randint(1,N,fanout),zeros(1,L)]; % 信息源:随机符号流,最后补上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);% depth_of_trellis长度为N+L,其中L=2
Eb=1;
Es=BPH*Eb;
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);% demod_input矩阵为number_of_out=64行,depth_of_trellis=N+L列,行表示频率号,列表示时间走势
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表示网格图的时间走势,depth_of_trellis长度为N+L,其中L=2
    f(i)=output(D+1,source(i)+1); % f(i)是i时刻的分支转移输出,即i时刻的跳频频率号,f(i)是跳频频率号序列
    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))
            rc(j+1)=sqrt(Es)+sgma*randn;% 改了,原来是sqrt(Es)+sgma*randn,现在信号幅度服从瑞利分布,瑞利分布的参数是sqrt(Es/2)
            rs(j+1)=sgma*randn;
        else
            rc(j+1)=sgma*randn;
            rs(j+1)=sgma*randn;
        end
    end
    jam_rayleigh=raylrnd(sqrt(Ej0/2));% 干扰音的幅度服从瑞利分布,瑞利分布的参数是sqrt(Ej0/2)
    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)
                rs(j+1)=rs(j+1)+jam_rayleigh*sin(theta);%改了,原来是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矩阵
    end
    D=nextstate(D+1,source(i)+1);
end
% ********************* 上面是信道和非相干解调模块 *****************%



% ****************** 逐符号(Symbol-by-Symbol)译码模块 ********************* %
[Max_energy freq_number]=max(demod_input,[],1);% 找出每跳时间内能量最大的那个频率号,存放于freq_number矢量中

% *************** End of Symbol-by-Symbol Decoding Modular **************** %

for i=1:N  %计算逐符号译码之后的符号错误率,发送频率序列与接收频率序列比较
    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循环对应的end
close(WTbar)
waitbar(pp/length(Eb_to_Nj_in_dB),WTbarpp)
end % 与最最外面for pp=1:length(Eb_to_Nj_in_dB)信干比循环相对应的end
close(WTbarpp)

Ps=sum(num_of_err,2)/(N*times) % 计算符号错误率,N 是一遍的符号流长度,总共进行times遍,共N*times个符号
Pb=Ps*fanout/(2*(fanout-1))