第5章差错控制系统的MATLAB/Simulink仿真 差错控制技术是指在发送端加入差错控制码元,即监督码元,利用监督码元和信息码元之间的某种确定关系来进行自动纠错和检错。一般来说,增加的监督码元越多,传输效率就越低,纠错和检错能力反而就越强。所以,差错控制技术是通过降低传输效率来提高信息在通信系统中传输的可靠性。本章介绍几种常用差错控制编码的基本原理,并通过MATLAB/Simulink对其进行仿真分析。 5.1基于线性分组码的差错控制系统仿真 对信源编码器输出的序列进行分组,并对每一组独立变换,称为分组码,记为(n,k)码,其中k表示每分组输入符号数,n为编码输出符号数。编码后的码组具有抗信道干扰的能力。若这种变换是线性变换,则称变换后的码组为线性分组码; 若变换是非线性的,则称变换后的码组为非线性分组码。常用的是线性分组码。 线性分组码具有如下两个性质: (1) 线性(包含全零码字,封闭性); (2) 最小码距等于除零码外的码字的最小码重。 为了具体说明线性分组码的基本原理,以(7,3)线性分组码为例。设(7,3)线性分组码为C=(c6,c5,c4,c3,c2,c1,c0),其中c6,c5,c4为信息位,c3,c2,c1,c0为监督位。将信息流分成每3位为一组,构成原码,即A=(a2,a1,a0),按下列线性方程进行编码: c6=a2 c5=a1 c4=a0 c3=a1a0 c2=a2a1 c1=a2a1a0 c0=a2a0(51) 写成矩阵形式,则可表示为: C=[c6,c5,c4,c3,c2,c1,c0]=[a2,a1,a0]1 0 00 1 00 0 10 1 11 1 01 1 11 0 1=A·G(52) 式中,A为原码; C为生成的线性分组码; G为生成线性分组码的生成矩阵。一般情况下,生成(n,k)线性分组码的生成矩阵大小为k×n。 由式(51)进一步变换可以得到监督位与信息位的关系: H·GT=0T(53) 式中,H为监督矩阵。监督矩阵中的每一个行向量与线性分组码中的任一码元内积为0,并且其中每一个行向量都线性无关。 1. 线性分组码的编码译码仿真 采用MATLAB仿真程序完成对(7,4)线性分组码的编码、解码,其中信息位和检验位的约束关系为c1=a1; c2=a2; c3=a3; c4=a4; c5=a1+a2+a3; c6=a2+a3+a4; c7=a1+a2+a4; 生成矩阵为G,校验矩阵为H,原码为A,生成码字为C,纠错后的码字为Cr。 clear all; G1=eye(4);   %生成4×4单位阵 G2=[1,0,1;1,1,1;1,1,0;0,1,1];   %约束关系 G=[G1,G2];   %生成矩阵G fprintf('生成矩阵为:G=') disp(G); A=[0,0,0,1;0,0,1,0;0,0,1,1;0,1,0,0;0,1,0,1;0,1,1,0;0,1,1,1;1,0,0,0;1,0,0,1;1,0,1,0;1,0,1,1;1,1,0,0;1,1,0,1;1,1,1,0;1,1,1,1;];   %A=[a1,a2,a3,a4]编码的原码 fprintf('原码为:A=') disp(A); C1=A*G; C=mod(C1,2);   %模2运算 fprintf('输出的编码为:C=') disp(C); H=gen2par(G)   %生成校验矩阵 fprintf('校验矩阵为:H=') %%以下输入接收到的码字,译出原码 Rev=input('请输入7位接收码字,用空格隔开:','s'); Rev=str2num(Rev)   %接收到的码字 S1=Rev*(H');   %S为校阵子 S=mod(S1,2); E=[1,1,1,1,1,1,1]; for i=1:7;   %取出H中的每一列,与S相加 Hi=H(:,[i]); Sum=S+Hi'; Sum=mod(Sum,2); if (all(Sum(:)==0)); %如果S与H的第i列之和Sum为0矩阵,则表示 %Rev中第i个码字有误 fprintf('接收码字中错误码位是第:'); disp(i) else E(1,i)=0; end; end; Cr=mod((Rev+E),2); fprintf('正确接收码字:Cr='); disp(Cr); 程序运行结果: 生成矩阵为: G = 1000101 0100111 0010110 0001011 原码为: A = 0001 0010 0011 0100 0101 0110 0111 1000 1001 1010 1011 1100 1101 1110 1111 输出的编码为: C = 0001011 0010110 0011101 0100111 0101100 0110001 0111010 1000101 1001110 1010011 1011000 1100010 1101001 1110100 1111111 校验矩阵为: H = 1110100 0111010 1101001 请输入7位接收码字,用空格隔开:0 1 1 1 1 0 1 接收码字中错误码位是第:2 正确码字为:Cr=0011101 2. 线性分组码信号的传输仿真 采用MATLAB和Simulink交互的方式仿真实现线性分组码编码后的信号在通信系统中的传输。Simulink仿真模型如图51所示。模型中Bernoulli Binary Generator(伯努利二进制序列产生器)模块产生的信源序列经过Binary Linear Encoder(二进制线性编码器)进行线性分组码编码。编码后的序列经过Binary Symmetric Channel(二元对称信道)传输,该信道具有误码概率。在接收端进行译码。译码后的序列和信源序列输入Error Rate Calculation(误码率统计)模块,统计接收端误码率。 图51线性分组码的Simulink仿真模型 模型中各模块的主要参数设置见表51。 表51线性分组码Simulink仿真参数 模 块 名 称 参 数 名 称 参数值 Bernoulli Binary Generator(伯努利二进制序列产生器) Probability of zero 0.5 Initial seed 10000 Sample time 1 Framebased output Checked Samples per frame 4 Binary Linear Encoder(二进制线性编码器) Generator matrix [[1 1 0; 0 1 1; 1 1 1; 1 01]eye(4)] Binary Symmetric Channel(二元对称信道) Error probability errB Initial seed 2137 Error Rate Calculation(误码率统计) Receive delay 0 Computation delay 0 Computation mode Entire frame Output data port To Workspace(工作区间) Variable name s Limit data point to last inf Decimation 1 Sample time -1 Save format Array 本例仿真采用MATLAB和Simulink交互的方式,分析经过线性分组码编码的信号在不同误码率的信道中传输后,接收端收到信号的误码率情况,代码如下: Clear; x=0:0.01:0.05; y=x; hold off; fori=1:length(x) errB=x(i); sim('m');   %实现MATLAB和Simulink模块的交互 y(i)=mean(s); end plot(x,y); hold on grid on; xlabel('信道误码率'); ylabel('接收端误码率Pe '); 图52为线性分组码的MATLAB和Simulink交互仿真结果。通过对图中信道误码率和接收端信号误码率的分析,可以看出使用线性分组码的差错控制仿真编码以后,差错率明显下降,比如信道误码率在0.04时,接收端信号误码率不足0.015。 图52线性分组码的传输仿真结果 5.2基于循环码的差错控制系统仿真 设C是某(n,k)线性分组码的码字集合,如果将任一码字c=(cn-1,cn-2,...,c0)向左移一位,记为c(1)=(cn-2,cn-3,…,c0,cn-1)也属于码集C,则该线性分组码为循环码。循环码的特点是具有循环性,即任何许用码字的循环移位仍然是一个许用码字。 对于一个任意长为n的码字c=(cn-1,cn-2,…,c1,c0),用多项式的形式表示,即c(x)=cn-1xn-1+cn-2xn-2+…+c1x+c0,此多项式即为码多项式,系数不为0的x的最高次数称为多项式c(x)的次数或阶数。在进行码多项式简单运算时,所得系数需进行模2运算。这就是循环码的编码原理。 构造循环码主要是要找出一组线性分组码的生成矩阵,将信息位与生成矩阵G相乘得到码字。由于循环码的码字多项式是生成多项式g(x)的倍式,且根据线性码生成矩阵的特性,(n,k)码的生成矩阵可以由(n,k)码k个不相关的码组构成。根据以上两点,可以挑选出k个线性不相关的循环码的码多项式。 接收端检测时,用码多项式y(x)除以生成多项式g(x)。若不能除尽,则说明接收码不属于循环码,在传输过程中发生错误,即y(x)g(x)=Q(x)+R(x)g(x)。因此,可用其余式R(x)是否为零来判断接收码组中是否有错。但当码字中错误位数超过循环码的检错能力时,接收码多项式仍有可能被g(x)整除,利用以上方法则不能检测出误码。 1. 循环码的编码仿真 下面程序完成(7,4)循环码的编码。代码中cyclpoly(n,k,'all')返回(n,k)循环码的所有生成多项式(1个生成多项式为返回矩阵的1行); cyclgen(n,g)返回循环码的监督矩阵和生成矩阵,其中g为生成多项式向量; rem(msg*G,2)返回循环码的所有需用码组,其中G为生成矩阵,msg为信息矩阵。 clear all; close all; n=7; k=4; p=cyclpoly(n,k,'all');   %产生循环码的生成多项式 [H,G]=cyclgen(n,p(1,:));   %产生循环码的生成矩阵和校验矩阵 Msg=[0 0 0 0;0 0 0 1;0 0 1 0;0 0 1 1;0 1 0 0;0 1 0 1;0 1 1 0;0 1 1 1; 1 0 0 0;1 0 0 1;1 0 1 0;1 0 1 1;1 1 0 0;1 1 0 1;1 1 1 0;1 1 1 1 ]; C = rem(Msg*G,2) 运行结果: C = 0000000 0110001 1100010 1010011 1110100 1000101 0010110 0100111 1011000 1101001 1101001 0111010 0001011 0101100 0011101 1001110 1111111 2. 循环码信号的传输仿真 循环码的Simulink仿真模型如图53所示。Bernoulli Binary Generator模块产生的信源序列经过Binary Cyclic Encode编码后在Binary Symmetric Channel中传输。Error Rate Calculation模块将接收端Binary Cyclic Decode译码后的序列和信源序列输入进行比较,统计误码率。 图53循环码Simulink仿真模型 表52是模型中各模块的主要参数设置。 表52循环码Simulink仿真参数 模 块 名 称 参 数 名 称 参数值 Bernoulli Binary Generator Probability of zero 0.5 Initial seed 10000 Sample time 1 Samples per frame 4 Binary Cyclic Encoder Codeword length N 7 Message length K 4 Binary Symmetric Channel Error probability errB Initial seed 2137 Binary Cyclic Decoder Codeword length N 7 Message length K 4 Error Rate Calculation Receive delay 0 Computation delay 0 Computation mode Entire frame Output data port 续表 模 块 名 称 参 数 名 称 参数值 To Workspace Variable name s Limit data point to last inf Decimation 1 Sample time -1 Save format Array 本例仿真同样采用MATLAB和Simulink交互的方式,利用5.1节线性分组码的MATLAB和Simulink交互代码调用图53循环码Simulink仿真模型,可以完成不同信道误码率下接收端的误码率统计,仿真结果见图54。从图中可以看出,循环码同样降低了接收端的误码率,提高了通信质量。 图54循环码信号传输的仿真结果 5.3基于卷积码的差错控制系统仿真 卷积码是结合本码组和其他码组之间的关系来进行编码的。卷积码通常记作(n,k,N)的形式,其中n为编码后码组的长度,k为输入的信息位的长度,N为编码约束长度。约束长度就是编码过程中互相约束的码组个数。卷积码编码的码字不仅与本码组中k个码字有关,同时和前(N-1)组输入的信息码字有关。相互有关系的码字一共有N×n个。约束长度越长纠错能力就越强,但码率会因此降低。R=k/n是卷积码的码率,表示信息位在所有要传输的码字中占有的比重。卷积码的一般结构如图55所示。 图55卷积码编码器结构 卷积码中k和n的值通常比较小,延时也相应小,适合串行传输。 下面利用MATLAB和Simulink两种方式进行仿真,对输入信号进行卷积编码,经过AWGN信道传输,接收端解码后统计传输误码率。 1. 卷积码的MATLAB仿真 下面采用(2,1,9)卷积码。程序利用randint函数产生信源序列,经过convenc卷积编码后,进行BPSK调制并在AWGN信道中传输。卷积码生成多项式为G0=561(八进制),G1=753(八进制)。接收端BPSK解调后利用软判决滑动窗维特比译码,译码深度为40。 clear all; close all SNRdB=0:0.5:3;   %设置信噪比范围,0~3dB declen=40;   %译码深度 SNRnum=length(SNRdB);   %信噪比数目 iter=10;   %每个信噪比下的迭代次数 for i=1:SNRnum   %循环内计算每个信噪比下的误码率 for j=1:iter %每个信噪比下迭代计算误码率10次,再求平均误码率 trel=poly2trellis(9,[561 753]);   %卷积码(2,1,9)网格图,约束长度为9 siglen=1000000;   %设置信号长度 msg=randi([0,1],siglen,1);   %生成0,1序列,长度同信号长度 encode=convenc(msg,trel,0);   %从0状态开始作卷积编码 I=0.5*ones(siglen*2,1); y=encode-I; bpsk=sign(y);   %BPSK调制 channelout=awgn(bpsk,SNRdB(i));   %添加高斯白噪声,AWGN信道 debpsk=channelout*0.5+0.5;   %解调 parti=0:.15:.9;   %设置量化等级划分 codebk=0:7;   %设输出等级 [x,qcode]=quantiz(debpsk,parti,codebk); %量化,准备维特比软判决 %维特比译码,量化级数=2^3 decode= vitdec(qcode',trel,declen,'cont','soft',3); %计算本次BER [errorbit,errorrate(j)] = biterr(decode(declen+1:end),msg(1:end-declen)); end BER(i)=sum(errorrate)/iter;   %求平均BER end semilogy(SNRdB,BER);   %绘制不同信噪比下的误码率 xlabel('信噪比SNR(dB)'); ylabel('误码率'); grid on; 图56为卷积码MATLAB仿真结果。从图中可以看出,接收端误码率随着信道信噪比的提高迅速降低。 图56卷积码MATLAB仿真结果 2. 卷积码的Simulink仿真 图57为卷积码的Simulink仿真模型。信源序列由Bernoulli Binary Generator模块产生,经过Convolution Encoder编码后进行BPSK调制,输入AWGN信道。接收端将信号BPSK解调后,采用Viterbi Decoder译码。译码后的序列和信源序列一起输入Error Rate Calculation模块统计误码率。 图57卷积码的Simulink仿真模型 表53为卷积码的Simulink仿真参数设置。 表53卷积码Simulink仿真参数 模 块 名 称 参 数 名 称 参数值 Bernoulli Binary Generator Probability of zero 0.5 Initial seed 61 Sample time 0.02/268 Framebased output Checked Samples per frame 268 BPSK Modulator Baseband(BPSK调制) Phase offset(rad) 0 Samples per symbol 1 AWGN Channel(加性高斯白噪声信道) Initial seed 67 Mode Signal to noise ratio(SNR) SNR(dB) SNR Input signal power(watts) 1 Convolution Encoder Trellis structure STRUCTURE Operation mode Truncated(reset every frame) Selector(数据选通器) Number of input dimensions 1 Index mode Onebased Index Option Index Vector(dialog) Index [1] Viterbi Decoder Trellis structure STRUCTURE Decision type Hard Decision Operation mode Truncated Traceback depth 34 本例通过MATLAB和Simulink交互的方式,完成不同码率、不同信噪比下接收端误码率的计算,代码如下: x=-10:2; %x表示信噪比 %卷积方式分别取1/3卷积和1/2卷积 A=[poly2trellis(9, [557 663 711]),poly2trellis(7, [171 133])]; %不同卷积方式、信噪比下重复调用图5-7的Simulink模型 for j=1:2 STRUCTURE=A(j); for i=1:length(x) SNR=x(i); %信道的信噪比依次取x中的元素 sim('book_c'); %实现MATLAB和Simulink模块的交互 y(j,i)=mean(BitErrorRate); %计算BitErrorRate的均值 end end semilogy(x,y(1,:),'r',x,y(2,:),'b') %绘图,采用对数坐标 xlabel('信道信噪比'),ylabel('接收端误码率'); legend('1/3卷积码','1/2卷积码'); grid on 图58为卷积码的MATLAB和Simulink交互仿真结果。图中两条曲线分别表示不同码率的卷积码在不同信噪比下误码率性能,其中上面的曲线码率为1/2,下面的曲线码率为1/3。 图58卷积码的仿真结果 5.4基于循环冗余码的差错控制系统仿真 为了满足实际中对n、k取值的多样性要求,通常在传送码字的后部预留一定的空间,用于差错校验。循环冗余码是最常见的校验码,由信息位和校验位两部分组成。其编码方法如下: (1) 移位: 将k比特原码左移r位,形成k+r=n位。 (2) 相除: 用生成多项式g(x),以模2除的方式去除移位后的式子,得到的余数就是校验码。 接收端收到数据后对其进行CRC校验,方法是将整个数据串当作一个整体去除生成矩阵,判断循环冗余校验器产生的余数。若余数为零,则说明接收正确; 否则接收错误,并不纠错。 循环冗余校验编解码的MATLAB仿真: clear all; close all; %编码 crcmsg=input('请输入信息原码,空格隔开:','s'); crcmsg=str2num(crcmsg); msglen = length(crcmsg);   %原码长度 crcgen =input('请输入生成多项式g(x),空格隔开:','s'); %如g(x)=x^5+x^4+x+1,输入1 1 0 0 1 1 crcgen=str2num(crcgen); critlen=length(crcgen); critbit = zeros(1,critlen);   %添加冗余比特 crcencode = [crcmsg zeros(1,critlen)];   %原码左移,左移位数为校验位长度 crcmsg = [crcmsg critbit]; divddvec = crcmsg;   %被除数向量 for k = 1:msglen   %开始循环计算长除得到最终余数 added = zeros(1,msglen-k+1); % 生成多项式向量左移位数,准备模2除 divvec = [crcgen added];   %除数向量 if divddvec(1) == 0   %被除数第1位为0,不必除 divvec = zeros(1,length(divvec)); end divddvec = bitxor(divvec,divddvec);   %模2除,等同异或 divvec = crcgen;   %恢复除数 divddvec(1) = [];   %去除被除数第1位 end crclen= length(crcencode) divdlen= length(divddvec); divddvec = [zeros(1,crclen-divdlen),divddvec]; %余数序列 crcencode = crcencode + divddvec   %编码后序列 %msglen = length(crcencode);   % 得到冗余编码的长度 %解码 divddvec = crcencode;   % 余数初始化 crcdecode = crcencode;   % 解码初始化 divnum = crclen-length(crcgen)+1;   % 长除的循环次数 for k = 1:divnum   % 计算余数 added = zeros(1,divnum-k); divvec = [crcgen added];   %除数向量 if divddvec(1) == 0   %被除数第1位为0,不必除 divvec = zeros(1,length(divvec)); end divddvec = bitxor(divvec,divddvec);   %模2除,等同异或 divvec = crcgen;   % 恢复除数 divddvec(1) = [];   %去除被除数第1位 end if sum(divddvec) == 0   % 校验正确 crcdecode = crcencode(1:divnum-1) else fprintf('校验错误'); end 运行结果: 请输入信息原码,空格隔开:1 1 1 1 0 0 0 1 0 请输入生成多项式g(x),空格隔开:1 1 0 1 1 crclen =14 crcencode =1 1 1 1 0 0 0 1 0 0 0 1 0 0 crcdecode = 1 1 1 1 0 0 0 1 0 5.5基于LDPC码的差错控制系统仿真 低密度校验码也称LDPC码,属于信道编码的一种,分为规则码和非规则码。因为低密度校验码编码方式不止一种,所以低密度校验码码的结构也会有所分类,码结构的不同或者选用低密度校验码译码算法不同,这些因素都会影响LDPC码的性能表现。在传统的信道中,如果信道中的噪声水平低于一个特定值,LDPC码的传输错误概率会随着码长的增加而减小。理论上,如果码长变得无穷大,传输的错误概率会无限趋近于零,只是在物理上这个理念并不能被实现。受益于LDPC码本身的编码方式,其自带抗突发性错误这个特性,所以当其运用在通信系统中时是不需要交织器的,结果就是减少了通信系统中的时延,从而提高通信系统的性能。 低密度校验(LDPC)码具有以下特点: (1) LDPC码译码具有很低的复杂度,所以不会因为码长的增加而急剧增加运算量,译码算法不仅在理论上,在物理上也是可以实现的,并且实际译码的仿真性接近于理论分析。 (2) LDPC码译码方法采用迭代译码算法,在物理上可以实现并行操作,并且译码速度高于其他编码。 (3) LDPC码拥有大的吞吐量,可以提高传输系统的传输效率,尤其在硬件模块更能体现这个优点。 (4) LDPC码的奇偶校验矩阵具有稀疏性,译码复杂度随着码长的改变而线性改变,所以不会出现码长过大而导致计算复杂度指数型增长这种情况。因为编码的码结构特性本身就自带抗突发差错这个能力,不像Turbo码那样需要交织器引入就拥有随机性,因此没有因交织器的存在而带来延时。 LDPC码也是一种分组码,其校验矩阵只含有很少量非零元素。正是校验矩阵H的这种稀疏性,保证了译码复杂度和最小码距都只随码长呈现线性增加。校验矩阵H的每一行对应一个校验方程,每一列对应码字中的一比特。因此,对于一个二进制码,如果它有m个奇偶校验约束关系,码字的长度为n,则校验矩阵是m×n的二进制矩阵。对于m×n维校验矩阵为H,当且仅当向量c=[c(1) c(2) … c(m)]满足 H·cT=0(54) 时,它才是该码的一个有效码字。 LDPC码可通过伪随机的方法构造,需要在给定一些设计规范后得到所有LDPC码的集合,而不仅限于如何选择一些特殊的校验矩阵使之满足设计规范。LDPC码常常通过Tanner图来表示,而Tanner图所表示的其实是LDPC码的校验矩阵。Tanner图包含两类顶点: n个码字比特顶点(称为比特节点),分别与校验矩阵的各列相对应; m个校验方程顶点(称为校验节点),分别与校验矩阵的各行相对应。校验矩阵的每行代表一个校验方程,每列代表一个码字比特。所以,如果一个码字比特包含在相应的校验方程中,那么就用一条连线将所涉及的比特节点和校验节点连起来,所以Tanner图中的连线数与校验矩阵中1的个数相同。图59是式(55)校验 图59LDPC码的Tanner图 矩阵的Tanner图,其中上面的圆形节点表示校验节点,下面的圆形节点表示比特节点,黑线表示的是一个6循环。 H=110100 011010 100011 001101(55) 除了校验矩阵H是稀疏矩阵外,LDPC码本身与任何其他的分组码并无二致。如果现有的分组码可以被稀疏矩阵所表达,那么用于LDPC码的迭代译码算法也可以成功地移植到它身上。不同的是,LDPC码的设计是以构造一个校验矩阵开始的,然后才通过它确定一个生成矩阵进行后续编码。 LDPC编码方法采用高斯消元算法,算法流程图如图510所示。 图510高斯消元算法流程图 LDPC码的译码方法采用概率BP译码算法,流程图如图511所示。 图511概率BP译码算法过程 LDPC码的Simulink仿真模型如图512所示。Bernoulli Binary Generator模块采样频率是32400。信道模块和程序设计选用的信道模型一样,都是选用的加性高斯白噪声信道(AWGN),调制/解调模块选用的是二进制相移键控(BPSK),误码率的计算就用Error Rate Calculation模块,这个模块的作用原理就是把信源的原始数据接入Tx端口,经过信道传输后的输出信号接入Rx端口,将输出信号与原始信号进行对比来计算误码率。对比后的数据就用输出端口接入Display模块来显示误码率的大小。 图512LDPC码Simulink仿真模型 表54是模型中各模块的主要参数设置。 表54LDPC码Simulink仿真参数 模 块 名 称 参 数 名 称 参数值 Bernoulli Binary Generator Probability of zero 0.7 Source of initial seed 61 Sample time 1 Samples per frame 32400 LDPC Encoder Paritycheck matrix dvbs2ldpc(1/2) BPSK Modulator Phase offset (rad) 0 Output type double AWGN Channel Initial seed 67 Mode Signal to noise ratio(Eb/No) Eb/No(dB) 10 Number of bits per symbol 1 Input signal power(watts) 1 Symbol period (s) 1 LDPC Decoder Output format Information part Decision type Soft decision Number of interation(迭代次数) 50 BPSK Demodulator Decision type Soft decision Phase offset (rad) 0 Error Rate Calculation Receive delay 0 Computation delay 0 Computation mode Entire frame Output data port 图513是信噪比为10dB时LDPC码Simulink仿真系统的输出结果图,从图中可以看出LDPC码经过BPSK调制/解调后输出消息的误码率较低。 图513LDPC码信号传输的仿真结果