multi_rician_vita_rscoded.m

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

Aj=sqrt((2*Kj*Ej0)/(Kj+1));% //the amplitude of the LOS component of an individual jamming tone
sigmaS=sqrt(Es/(Ks+1));% //represent the multipath component of the DFH signal tone
sigmaJ=sqrt(Ej0/(Kj+1));% //represent the multipath component of an individual jamming tone
% sgma=sqrt(Eb/(2*Eb_to_No));% AWGN的均方根
sigmaN=sqrt(Eb/Eb_to_No);% //Attention: here the sigmaN^2 should be equal to No
demod_input=zeros(number_of_out,depth_of_trellis);
f=zeros(1,depth_of_trellis);
D=0; % D 记录网格图的当前状态,这里初始状态是0状态

% *******************信道和非相干解调部分:加多音干扰和噪声,然后非相干解调****************** %
for i=1:depth_of_trellis % i表示网格图的时间走势
    rc=zeros(1,number_of_out);
    rs=zeros(1,number_of_out);
    f(i)=output(D+1,source_coded(i)+1); % f(i)是i时刻的分支转移输出,即i时刻的跳频频率号,频率号范围是[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))
            as=raylrnd(sigmaS);% //amplitude of the multipath component of the DFH signal tone
            theta_s=2*pi*rand;% //random phase of the LOS component of the signal tone
            thetaS_mp=2*pi*rand;% //random phase of the multipath component of the signal tone
            rc(j+1)=As*cos(theta_s)+as*cos(thetaS_mp)+sigmaN*randn;
            rs(j+1)=As*sin(theta_s)+as*sin(thetaS_mp)+sigmaN*randn;
        else
            rc(j+1)=sigmaN*randn;
            rs(j+1)=sigmaN*randn;
        end
    end
    for k=1:Q
        aj=raylrnd(sigmaJ); % //amplitude of the multipath component of an individual jamming tone
        theta_j=2*pi*rand;% //random phase of the LOS component of each jamming tone
        thetaJ_mp=2*pi*rand;% //random phase of the multipath component 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);
%                 rs(j+1)=rs(j+1)+jam_rayleigh*sin(theta_j);
                rc(j+1)=rc(j+1)+Aj*cos(theta_j)+aj*cos(thetaJ_mp);% // add the jamming tone components to 'rc'
                rs(j+1)=rs(j+1)+Aj*sin(theta_j)+aj*sin(thetaJ_mp);% // add the jamming tone components to 'rs'
            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 Demodulation Modular **************** %


% ************* The Viterbi Sequence Detection 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=zeros(1,BPH*(depth_of_trellis-L));% length(decoder_output)=BPH*2N
decoder_output=zeros(1,depth_of_trellis-L);% depth_of_trellis=length(source_coded)=3002=3000+L,存放3000个fanout=4进制符号流
for i=1:depth_of_trellis-L
    dec_output_deci=input(state_sequence(1,i)+1,state_sequence(1,i+1)+1);% 输出的十进制符号寄存器dec_output_deci(其实是fanout=4进制的符号)
                                                                         % 因为下面一句dec_output_bin=deci2change(dec_output_deci,BPH,2)显然是
                                                                         % 把dec_output_deci转化为 BPH 位二进制的,所以dec_output_deci
                                                                         % 应该是2^BPH=fanout进制而不是十进制的
    decoder_output(i)=dec_output_deci; % decoder_output矩阵存放的是G函数维特比译码后得到的3000个普通域的fanout进制的符号
%     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
% ****************** End of the Viterbi Sequence Detection Modular ****************** %


% ****************** Normal Domain to Galois Domain Change Modular ****************** %
% 解交织和译码之前要先将普通域的fanout=4进制的符号转化为GF(16)域的16进制符号,准备进行基于GF(16)域的16进制符号的解交织和RS译码
decoder_output_hex=zeros(1,nRS*N/kRS);% 存放普通域的16进制符号
decoder_output_GF=zeros(1,nRS*N/kRS);% 存放GF域的16进制符号
for i=1:(depth_of_trellis-L)/2   % 即 i=1:1500
    w1=decoder_output(2*(i-1)+1); % decoder_output中有3000个fanout=4进制符号,把两位4进制符号转化为1个普通域的16进制符号
    w2=decoder_output(2*i);       % 这段"把两位4进制符号转化为1个普通域的16进制符号"的程序通过测试,正确!
    w=4*w1+1*w2;
    for j=0:2^M-1
        if (w==j)
            decoder_output_hex(i)=j;
        end
    end
end
decoder_output_GF=gf(decoder_output_hex,M);% 把普通域的16进制符号转化为GF(16)域的16进制符号,总共1500个符号,准备解交织后进行RS译码
% ****************** End of the Domains Change Modular ***************** %

% % ***********************随机解交织*********************%
% 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
% % ************************************************************************%

%**************解交织部分(基于16进制GF(16)域符号(即RS码元)的解交织)***************%
deint_output=zeros(N/kRS,nRS);
deint_output_GF=gf(deint_output,M); % 存放解交织后的GF(16)域的16进制符号矩阵(尺寸为100*15),是 RS 译码器的输入矩阵
for i=1:nRS*N/kRS
    deint_output_GF(i)=decoder_output_GF(i);% 对!必须将解调后的序列按列存放为100行15列的矩阵
end
%************** End of the De-interleaving Modular **************** %

% % ***********(8,4)译码部分***********%
% deco_input=zeros(BPH*2*N,1); % 列矢量
% 
% deco_input(1:BPH*2*N)=deint_output_bin(1:BPH*2*N);
% deco_output=decode(deco_input,8,4,'linear',G)';% (8,4)译码器的译码输出(长为BPH*N)个比特,行矢量)
% 
% % ***********************************%

%**************** RS 码的译码部分(N,K,Dmin)=(15,7,9)的 RS 码*****************%
deco_output_GF=rsdec(deint_output_GF,nRS,kRS);% RS 译码输出矩阵GF(16)域(尺寸100*7)
                                           % 注意这个deco_output_GF矩阵对应于RS编码前的source矩阵,矩阵source的第一行存放的是原始信息符号流source_sequence
                                           % 序列的第1到7个符号,第二行存放的是source_sequence序列的第8到14个符号,依此类推.而这个deco_output_GF矩阵同样如此,
                                           % 它的每个位置上的元素与前面source矩阵相同位置上的元素应该是一一对应的!
                                           % 下一步应该把这个deco_output_GF矩阵的元素按行取出构成一个1*N的序列,与原始信息符号流source_sequence进行
                                           % 比较,计算符号错误率
%****************************************************************************%

deco_output_GF1=deco_output_GF';% 转置一下变成7行100列矩阵,每列为一个信息符号组    
deco_output_GF2=reshape(deco_output_GF1,1,N); % 把译码输出reshape(reshape按列取出)成1行N=700列的GF(16)域符号矢量,准备计算符号错误个数
 
for i=1:(M/BPH)*nRS*N/kRS   % 计算维特比译码后,RS译码之前的4进制符号错误数 
    if (decoder_output(i)~=codedsequence_quart(i))
        vita_symbol_err(pp,rep)=vita_symbol_err(pp,rep)+1;% 这时vita_symbol_err是维特比译码后,RS译码之前的2^BPH=4进制符号错误数
    end
end

for i=1:N    % N=700,700个16进制符号,计算最终RS译码后的16进制符号错误数
    if(deco_output_GF2(i)~=source_sequence(i))
       num_of_err(pp,rep)=num_of_err(pp,rep)+1; % 这时num_of_err是2^M=16进制符号错误数
   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)



PS_all=sum(num_of_err,2)/(N*times) % 计算出16进制符号错误率,其中N 是一遍的符号流长度,总共进行times遍,共N*times个符号
PB_all=2^(M-1)/(2^M-1)*PS_all

symerr_num_biterr=zeros(2*length(Eb_to_Nj_in_dB),times);% symerr_num_biterr矩阵中奇数行是某一信干比下维特比译码后,RS译码之前的2^BPH=4进制符号错误数
                                                        % 偶数行是对应的最终RS译码后的2^M=16进制符号错误数
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