multi_rayleigh_vita_84encoded_all.asv
来自「短波信道抗多音干扰的性能分析及其仿真」· ASV 代码 · 共 284 行
ASV
284 行
%function pb=multi_rayleigh_vita_84encoded_all(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
% 这个程序是正确的,带(8,4)纠错码,编码后在后面只补上L个0,即符号流总长度是2N+L而不是2(N+L),
% 而multi_awgn_vita_encode.m是带(8,4)纠错码但符号流总长度是2(N+L),这两个程序都对
%
% 这个程序是在带(8,4)纠错码的multi_awgn_vita_encodev.m基础上加了交织,并且是基于比特的分组(block)交织
%
%********系统参数********%
% Eb_to_Nj_in_dB=11;
Eb_to_Nj_in_dB=[0,3,6];
Eb_to_No_in_dB=13.35;
BPH=2;
number_of_states=16;
Q=4; % 干扰音的个数
N=1000; % 每次符号流长度
times=10; % 重复做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级是其状态位,与网格图中的状态一一对应(注意并不是与跳频频点一一对应)
% pb=zeros(1,times);
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); %%%%%%%%%%%%其实这里是把 2^BPH 进制的符号source(i)转换成BPH个二进制比特,前面产生的符号流source=randint(1,N,fanout)
%产生的符号流原来本身也就是fanout=2^BPH进制的,这正好相吻合,所以这里并不是把十进制的符号source(i)转化为
%BPH位二进制比特,而是把2^BPH进制的符号source(i)(显然source=randint(1,N,fanout)表明source(i)
%的取值范围被限定在0到fanout-1之间)转化为BPH位二进制比特,
end
else
dsource=source(1:N);%*****
end
%******************************************************%
% *******************差错控制编码部分(8,4)码**********************%
% G=[1 0 0 0 1 0 1 1;
% 0 1 0 0 1 1 1 0;
% 0 0 1 0 1 1 0 1;
% 0 0 0 1 0 1 1 1];
G=[1 0 0 0 1 0 1 1;0 1 0 0 1 1 0 1;0 0 1 0 0 1 1 1;0 0 0 1 1 1 1 0]; % (8,4,4)扩展汉明码的生成矩阵G
% H=[1 1 1 1 1 1 1 1;1 1 0 1 1 0 0 0;0 1 1 1 0 1 0 0;1 0 1 1 0 0 1 0]; % (8,4,4)扩展汉明码的一致校验矩阵H
% trt=syndtable(H); % (8,4,4)扩展汉明码的伴随式译码表
dsource_coded=encode(dsource,8,4,'linear',G)';
% source_coded0=zeros(1,2*N);
% *****************************************************************%
% *********************交织部分(块交织)(基于比特的交织)****************%
% 仅适用于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); % 其实是把BPH个二进制比特转化为一个2^BPH进制的符号,这里BPH=2,就是把2个二进制比特转化为1个四进制符号(是对的)
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符号,长度为2N+L个fanout进制符号
% ******************* 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);%*******************% depth_of_trellis长度为2N+L
Eb=1;
Es=Eb*BPH*(1/2);
% 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]
thyta=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(thyta);%改了,原来是rc(j+1)+sqrt(Ej0)*cos(thyta)*2222222222222222222222222222222222222222222222222222222222*
rs(j+1)=rs(j+1)+jam_rayleigh*sin(thyta);%改了,原来是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
% ******************************************************************** %
% ***********************维特比译码模块******************************* %
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
for i=1:depth_of_trellis-L % i=1:2N
dec_output_deci=input(state_sequence(1,i)+1,state_sequence(1,i+1)+1);% 输出的十进制符号寄存器dec_output_deci(其实是fanout进制的符号)
% 因为下面一句dec_output_bin=deci2change(dec_output_deci,BPH,2)显然是
% 把dec_output_deci转化为 BPH 位二进制的,所以dec_output_deci
% 应该是2^BPH=fanout进制而不是十进制的
vitdec_output(i)=dec_output_deci; % 存放维特比译码输出的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
% % ***********************随机解交织*********************%
% 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(80,50);
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
% ************************************************************************%
% ***********(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)个比特,行矢量)
% ***********************************%
for i=1:2*N % % 计算维特比译码后,(8,4)译码之前的fanout=4进制符号错误数,共有2N个符号
if(vitdec_output(i)~=source_coded(i))
vita_symbol_err(pp,rep)=vita_symbol_err(pp,rep)+1;% 这时vita_symbol_err是维特比译码后,(8,4)译码之前的4进制符号错误数
end
end
for i=1:N*BPH
if(dsource(i)~=deco_output(i))
num_of_err(pp,rep)=num_of_err(pp,rep)+1;
end
end
end % 与最外层"多少遍rep" for循环对应的end
end % 对应于最最外面pp信干比循环的end
PB_all=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矩阵中奇数行是某一信干比下维特比译码后,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⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?