通信原理实验报告(二) 含matlab程序

通信原理实验报告

实验二 QPSK 通信系统的Monte Carlo 仿真

一、实验目的

1. 提高独立学习的能力

2. 培养发现问题，解决问题，分析问题的能力

3. 学习Matlab 的使用

4. 掌握4PSK 通信系统的Monte Carlo 仿真方法

5. 掌握4PSK 通信系统的组成原理

6.

比较编码信号与未编码信号在随机信道中的传输，加深对纠错编码原理的理解。

二、实验原理

1． 调制解调原理 一组M 载波相位调制信号波形的一般表示式为：

()()1,...,1,0),22cos(-=-

=M m M

m

t f t Ag t u c T m ππ 式中()t g T 是发送滤波器的脉冲形状，决定了传输信号的频谱特性，A 是信号的幅度，PSK 信

号对所有的m 都具有相同的能量，即 ()s m

m dt t u εε==

?

∞

∞

-2

其中的s ε代表每个传输符号的能量。

当()t g T 是一个矩形脉冲时，定义为

()T t T

t g T ≤≤=

0,2

这时在符号区间T t ≤≤0内传输的信号波形可以表示为（用s A ε=）

()1,...,1,0),22cos(2-=-=

M m M

m

t f T t u c s m ππε 把式中的余弦函数的相角看成两个相角的和，可以将上式表示为

()()()()()

()()

t s t s t f M m t g t f M m t g t u ms mc c T s c T s m 212sin 2sin 2cos 2cos ψ+ψ=??

?

??-??? ??=ππεππε 这里M

m

s M m s s ms s mc πεπε2sin

,2cos

==，而()t 1ψ和()t 2ψ是两个正交基函数，定义为 ()()()()()()t f t g t t f t g t c T c T ππ2sin ,2cos 21-=ψ=ψ

适当地将()t g T 归一化，就可以将这两个基函数的能量归一化到1。这样一个相位调制信号可以看作两个正交载波，起幅度取决于在每个信号区间内的相位，因此，数字相位调制信号在几何上可用mc s 和ms s 的二维向量来表示，即

??? ?

?

=M m M m s s s m πεπε2sin ,2cos

在AWGN 信道中，在一个区间内接受到的带通信号可以表示为

()()()()()()()()t f t n t f t n t u t n t u t r c s c c m m ππ2sin 2cos -+=+=

其中()t n c 和()t n s 是加性噪声的同相分量和正交分量，将这个接收信号与给出的()t 1ψ和()t 2ψ作相关，两个相关器的输出可以表示为：

??

?

?

?++=+=s s c s m n M

m n M

m n s r πεπε2sin ,2cos

这两个正交的噪声分量()t n c 和()t n s 是零均值，互不相关的高斯随机过程，这样

()()()()()

2

000

2

2N n E n E n n E n E n E s c s c s c =

====，方差为。 最佳检测器将接收到的信号向量r 投射到M 个可能的传输信号向量{}m s 之一上去，并选取对应与最大投影的向量，从而得到相关准则为()

m m s r s r C ?=

,。

由于二相相位调制与二进制PAM 是相同的，所以差错概率为

???

?

??

??=0

22N Q P b

psk ε 式中b ε是每比特能量。4PSK 可以看作是两个在正交载波上的二相相位调制系统，所以1个比特的差错概率与二相相位调制是一致的。

2． 信道纠错编码

在随机信道中，错码的出现是随机的，且错码之间是统计独立的。例如，由高斯白噪声引起的错码就具有这种性质。因此，当信道中加性干扰主要是这种噪声时，就称这种信道为随机信道。由于信息码元序列是一种随机序列，接收端是无法预知的，也无法识别其中有无错码。为了解决这个问题，可以由发送端的信道编码器在信息码元序列中增加一些监督码元。这些监督码元和信码之间有一定的关系，使接收端可以利用这种关系由信道译码器来发现或纠正可能存在的错码。在信息码元序列中加入监督码元就称为差错控制编码，有时也称为纠错编码。不同的编码方法有不同的检错或纠错能力。有的编码就只能检错不能纠错。

那么，为了纠正一位错码，在分组码中最少要增加多少监督位才行呢？编码效率能否提高呢？从这种思想出发进行研究，便导致汉明码的诞生。汉明码是一种能够纠正一位错码且编码效率较高的线性分组码。下面我们介绍汉明码的构造原理。

一般说来，若码长为n ，信息位数为k ，则监督位数r ＝n?k 。如果希望用r 个监督位构造出r 个监督关系式来指示一位错码的n 种可能位置，则要求

2r ? 1 ≥ n 或 2r ≥ k + r + 1 （8－1）

下面我们通过一个例子来说明如何具体构造这些监督关系式。

设分组码(n ，k)中k ＝4，为了纠正一位错码，由式（8－1）可知，要求监督位数r≥3。若取r=3，则n= k + r =7。我们用α6α5…α0 表示这7个码元，用S 1 、S 2 、S 3 表示三个监督关系式中的校正子，则S 1 S 2 S 3 的值与错码位置的对应关系可以规定如表18-1所列。

表 8-1

由表中规定可见，仅当一错码位置在α2 、α4 、α5 或α6 时，校正子S 1 为1；否则S 1 为0。这就意味着α2 、α4 、α5 和α6四个码元构成偶数监督关系

S 1 ＝α6⊕α5 ⊕α4 ⊕α2

（8－2） 同理，α1 、α3 、α5 和α6构成偶数监督关系

S 2 ＝α6⊕α5 ⊕α3 ⊕α1 （8－3）

以及α0 、α3 、α4 和α6构成偶数监督关系

S 3 ＝α6⊕α4 ⊕α3 ⊕α0 （8－4）

在发送端编码时，信息位α6 、α5 、α4 和α3 的值决定于输入信号，因此它们是随机的。监督位α2 、α1 和α0 应根据信息位的取值按监督关系来确定，即监督位应使上三式中S 1 、S 2 和S 3 的值为零（表示变成的码组中应无错码）

?

?

?

??

⊕⊕⊕⊕⊕⊕⊕⊕⊕000034613562456＝αααα＝αααα＝αααα

18－5）

由上式经移项运算，解出监督位

?

?

?

??

⊕⊕⊕⊕⊕⊕αα＝αααα＝ααααα＝α340635164562

（8－6）

给定信息位后，可直接按上式算出监督位，其结果如表8-2所列。

表 8-2

接收端收到每个码组后，先按式（8－2）～（8－4）计算出S 1 、S 2 和S 3 ，再按表8-2判断错码情况。按上述方法构造的码称为汉明码。表8-2中所列的（7，4）汉明码的最小码距d 0 ＝3，因此，这种码能纠正一个错码或检测两个错码。

3、两实验的系统框图如下：

（一）未加信道纠错编码的4PSK 调制通信系统

（二）信道纠错编码（7，4）汉明码+4PSK

调制的通信系统

三、实验内容

（一）未加信道纠错编码的4PSK调制通信系统

1、实验要求

1)根据系统组成框图用Matlab编写程序，并在每条语句后进行注释

2)最大投影点准则进行判决

a, 计算噪声方差2σ分别为0、0.1、0.5、1.0时的符号差错概率和比特差错概率；

b, 画出在每种2σ时，在检测器输入端1000个接收到的信号加噪声的样本（星座图）；

c, 分别画出数据点为1000、5000、10000、100000时的Monte Carlo仿真误比特率曲线和理论误比特率曲线，比较差别，分析数据点的数量对仿真结果的影响；

3)将检测器的判决准则改为最小距离法（星座图上符号间的距离），比较与上面结果的区

别。

4)分析整个通信系统各组成部分的原理、分析实验结果，总结概括。

2、实验程序

（1）主程序1

N=1000;%数据点数为1000

xulie=subserial(N);%生成序列

m=psktiaozhi(xulie);%4psk调制第一步，算出码元所对应的m

for j=1:N/2 %两路信号

a(j)=cos(pi*m(j)/2);

b(j)=sin(pi*m(j)/2);

end;

t=1;

G=char('*r','*g','*b','*y');%！！！！！！！！！！！！

for R=[0,0.1,0.5,1.0]

for i=1:N/2%生成噪声

[nc(i),ns(i)]=gnguass(0,sqrt(R));%生成噪声序列

end;

for j=1:N/2 %两路信号加噪声

Rc(j)=a(j)+nc(j);

Rs(j)=b(j)+ns(j);

end;

subplot(2,2,t)

hold on;

for i=1:N/2

plot(Rs(i),Rc(i),G(m(i)+1,:));axis([-4 4 -4 4]);%根据不同的m（i）值可以判断某一码元的来源，

从而赋予其相应颜色

end

hold off;

title('星座图');

s1=zuixiaojuli(Rc,Rs)%最小距离法

shuchuxulie1=pskjietiao(s1);%4psk解调

errorate1juli=byterror(xulie,shuchuxulie1,t)%比特错误

errorate2juli=fuhaoerror(xulie,shuchuxulie1,t)%符号错误

s2=zuidatouying(Rc,Rs);%最大投影法

shuchuxulie2=pskjietiao(s2);%4psk解调

errorate1touying=byterror(xulie,shuchuxulie2,t)%比特错误

errorate2touying=fuhaoerror(xulie,shuchuxulie2,t)%符号错误

t=t+1;

end;

（2）主程序2

N0=0:0.1:2;

Eb=1;%平均码元能量为1

Pb=0.5*erfc(sqrt(Eb./N0));%理论差错概率曲线公式

semilogy(N0,Pb);%画出N0与Pb（差错概率）的关系曲线，理论为蓝色，1000点红色，5000点黑色，10000点绿色，100000点黄色

hold on;

k=1;%控制大循环，每次循环改变点数

G=char('rd-','kd-','gd-','yd-');

for N=[1000,5000]%,10000,100000]

xulie=subserial(N);%生成序列

m=psktiaozhi(xulie);%4psk调制第一步，算出码元所对应的m

for j=1:N/2 %两路信号

a(j)=cos(pi*m(j)/2);

b(j)=sin(pi*m(j)/2);

end;

t=1;%控制循环，每次循环改变N值

for RR=0:0.1:2 %不同N值

R=sqrt(RR/2);%不同N值所对应的不同方差值

for i=1:N/2%生成噪声

[nc(i),ns(i)]=gnguass(0,R);

end;

for j=1:N/2 %两路信号加噪声

Rc(j)=a(j)+nc(j);

Rs(j)=b(j)+ns(j);

end;

s=zuidatouying(Rc,Rs);%最大投影法

shuchuxulie=pskjietiao(s);%4psk解调

byteerro(t)=0;%计算比特错误，计数变量清零

for i=1:N

if (xulie(i)==shuchuxulie(i))%比较输出序列与原序列，若相同无操作，若不同计数变量增一

byteerro(t)=byteerro(t);

else byteerro(t)=byteerro(t)+1;

end;

end;

errorate(t)=byteerro(t)/N;%算出误比特率

t=t+1;

end;

RR=0:0.1:2;

semilogy(RR,errorate,G(k));%画图显示误比特率与N值的函数关系

k=k+1;

end;

（3）高斯噪声产生子函数

function [gsrv1,gsrv2]=gnguass(m,sgma)

%本函数生成以高斯分布生成两相互正交的随机数

%m为期望，sgma为方差

if nargin==0

m=0;

sgma=1;

elseif nargin==1

sgma=m;

m=0;

end%以上几句的含义是，若未输入参量，则期望为0，方差为1，即正态分布；如输入一个参量，则该参量为sgma，期望0；

u=rand;%u为随机产生的一个数，在0~1中服从均匀分布

z=sgma*sqrt(2*log10(1/(1-u)));

u=rand;

gsrv1=m+z*cos(2*pi*u);

gsrv2=m+z*sin(2*pi*u);%通过上述四行公式将均匀分布映射成高斯分布

（4）序列产生子函数

function [source]=subserial(N)

%本函数的功能是产生一个长度为N的随机0、1序列

%N为序列长度，source为生成序列

source=zeros(1,N);

for i=1:N;

temp=rand;%产生随机数，在0~1上服从均匀分布

if (temp<0.5)

source(i)=0;%小于0.5时判决为0，否则为1

else

source(i)=1;

end

end

（5）4psk调制子函数

function [m]=psktiaozhi(xulie)

%本函数实现4psk调制第一步，即算出m序列

%xulie为原序列

N=length(xulie);

j=1;

for i=1:2:N-1

if (xulie(i)==0&xulie(i+1)==0)%通过码元的值判断得出m，共四种情况00、01、10、11，分别对应m=0、1、2、3

m(j)=0;

elseif (xulie(i)==0&xulie(i+1)==1)

m(j)=1;

elseif (xulie(i)==1&xulie(i+1)==0)

m(j)=2;

else (xulie(i)==1&xulie(i+1)==1)

m(j)=3;

end;

j=j+1;

end;

（6）4psk解调子函数

function [shuchuxulie]=pskjietiao(s)

%该函数实现4psk解调

%s为调制后序列，shuchuxulie为解调后序列

N=2*length(s);

j=1;

for n=1:2:N-1

if s(j)==0

shuchuxulie(n)=0;%通过s的值判断得出原码元的值，共四种情况s=0、1、2、3，分别对

应00、01、10、11

shuchuxulie(n+1)=0;

elseif s(j)==1

shuchuxulie(n)=0;

shuchuxulie(n+1)=1;

elseif s(j)==2

shuchuxulie(n)=1;

shuchuxulie(n+1)=0;

else shuchuxulie(n)=1;

shuchuxulie(n+1)=1;

end;

j=j+1;

end;

（7）比特错误统计及误比特率计算子函数

function [errorate1]=byterror(xulie,shuchuxulie,t)

%本函数实现序列误码数与误码率的计算

%errorate1为误码率，xulie、shuchuxulie为两序列，t为结果数组的下标，若误码率仅为一个数，可使t为一定值

N=length(xulie);

byteerro(t)=0;%比特错误

for i=1:N

if (xulie(i)==shuchuxulie(i))

byteerro(t)=byteerro(t);

else byteerro(t)=byteerro(t)+1;%比较输出序列与原序列，若相同无操作，若不同计数变量增一

end;

end;

errorate1(t)=byteerro(t)/N;%算出误比特率

（8）符号错误统计及误符号率计算子函数

function [errorate2]=fuhaoerror(xulie,shuchuxulie,t)

%本函数实现序列误符号数与误符号率的计算

%errorate1为误码率，xulie、shuchuxulie为两序列，t为结果数组的下标，若误码率仅为一个数，可使t为一定值

N=length(xulie);

fuhaoerro(t)=0;%符号错误

for i=1:2:N-1

if ((xulie(i)==shuchuxulie(i))&(xulie(i+1)==shuchuxulie(i+1)))

fuhaoerro(t)=fuhaoerro(t);

else fuhaoerro(t)=fuhaoerro(t)+1;%比较输出序列与原序列的码元，码元中的两符号若都相同无操作，否则计数变量增一

end;

end;

errorate2(t)=fuhaoerro(t)/(N/2);%误符号率

3、程序运行结果

（1）主程序1运行结果如下：

σ经4psk调制及在信道中叠加噪声后的星座图，图中不同颜色图1为序列在不同噪声方差2

决定于该数据点在调制之后加噪之前的调制值，红色点为（0，1），黄色点为（-1，0），蓝

σ值分别为0、色点为（0，-1），绿色点为（1，0）。图中，左上、右上、左下、右下图中2

0.1、0.5、1.0，由图可知，噪声方差2σ越大，星座图中数据点的离散程度就越大。

图二为两种判决方法在各种情况下的误比特率与误码率的统计结果，其中errorate1juli 代表最小距离法的误比特率，errorate2juli代表最小距离法的误码率，errorate1touying代表最大投影法的误比特率，errorate2touying代表最大投影法的误码率。每项由四列组成，从

σ值分别为0、0.1、0.5、1.0。由数据知，误码率总大于误比特率，且两种判决方法左至右2

所得最终两种错误率完全相同。

在我们看来，两种判决方法在本质上是统一的，对于最大投影点法实为向量点乘，在两向量模值相同时，大小即由夹角θ决定，夹角越小值越大。对应的，在最小距离法中，根据余弦定理，同样可得大小由夹角θ决定，夹角越小值越大，所以，两种方法的判决结果必然完全相同。

（2）主程序2运行结果

图三为序列点数为1000、5000点时N-Pe（噪声功率-误码率）曲线，其中蓝色为理论值，红色为1000点，黑色为5000点，很明显，黑色更接近理论曲线，即点数越大，实验曲线越接近理论曲线。

图四为序列点数为1000、5000、10000、100000点时的N-Pe（噪声功率-误码率）曲线，其中蓝色为理论值，红色为1000点，黑色为5000点、绿色为10000点、黄色为100000点。结果分析如图三。

（二）信道纠错编码（7，4）汉明码+4PSK调制的通信系统

1、实验要求

1）根据系统组成框图用Matlab编写程序，并在每条语句后进行注释

2）比较编码信号与未编码信号在随机信道中的传输性能，分析实验结果，总结概括

2、实验程序

（1）主程序

%本程序实现编码信号与未编码信号的性能比较

Ts=1;%码元宽度为1

N=1000;%数据点数为1000

xulie=subserial(N);%生成1000点序列

m=psktiaozhi(xulie);%4psk调制第一步，算出m

for j=1:N/2 %两路信号，4psk映射

a(j)=cos(pi*m(j)/2);

b(j)=sin(pi*m(j)/2);

end;

t=1;%控制大循环，每次循环改变R（sigma）

G=char('*r','*g','*b','*y');%生成字符串数组，在循环中方便显示参数的添加

for R=[0,0.5,1,2.0]

for i=1:N/2%生成噪声

[nc(i),ns(i)]=gnguass(0,sqrt(R));%生成噪声序列

end;

for j=1:N/2 %两路信号加噪声

Rc(j)=a(j)+nc(j);

Rs(j)=b(j)+ns(j);

end;

subplot(2,2,t)

hold on;

for i=1:N/2

plot(Rs(i),Rc(i),G(m(i)+1,:));%画星座图，根据不同的m（i）值可以判断某一码元的来源，从而赋予其相应颜色

axis([-4 4 -4 4]);

end

hold off;

title('星座图');

s=zuidatouying(Rc,Rs);%最大投影法

shuchuxulie=pskjietiao(s);%4psk解调

errorate1touying=byterror(xulie,shuchuxulie,t)%比特错误

t=t+1;

end;

figure;%新图，用来显示hamming编码后的星座图

hammingout=hamm(xulie);%调用函数生成汉明序列

M=length(hammingout);%M标志汉明序列长度

m_hamm=psktiaozhi(hammingout);%4psk调制第一步

for j=1:M/2 %两路信号

a(j)=cos(pi*m_hamm(j)/2);

b(j)=sin(pi*m_hamm(j)/2);

end;

t=1;%控制大循环，每次循环改变R（sigma）

for R=[0,0.1,0.5,1.0]

for i=1:M/2%生成噪声

[nc(i),ns(i)]=gnguass(0,sqrt(R));%生成噪声序列

end;

for j=1:M/2 %两路信号加噪声

Rc(j)=a(j)+nc(j);

Rs(j)=b(j)+ns(j);

end;

subplot(2,2,t)

hold on;

for i=1:M/2

plot(Rs(i),Rc(i),G(m_hamm(i)+1,:));axis([-4 4 -4 4]);%根据不同的m（i）值可以判断某一码元的来源，从而赋予其相应颜色

end

hold off;

title('星座图');

s_hamm=zuidatouying(Rc,Rs);%最大投影法

shuchuxulie_hamm=pskjietiao(s_hamm);%4psk解调

shuchuxulie2=decoder(shuchuxulie_hamm);%调用函数，汉明码解码

errorate1_hamm=byterror(xulie,shuchuxulie2,t)%比特错误

t=t+1;

end;

（2）汉明编码子函数

function [hammingout]=hamm(sourc)

%此函数实现汉明码的编码过程

%sourc为原序列，hamming为生成的（7，4）汉明码序列

%模二加算法也可先进行加运算再模二

N=length(sourc);

subsourc1=reshape(sourc,4,N/4);%四行N/4列

subsourc=subsourc1'%变为四列

G=[ 1 0 0 0 1 1 1;

0 1 0 0 1 1 0;

0 0 1 0 1 0 1;

0 0 0 1 0 1 1;];%生成矩阵

subout=subsourc*G;%生成新序列，七位一行

subout1=subout';%变为七位一列

subhammingout=reshape(subout1,1,7*N/4);%将数组变为一行

hammingout=rem(subhammingout,2);%模2，生成最终01序列

（3）汉明解码子函数

function [sourout]=decoder(hamm)

%此函数完成汉明码的译码，具体分三步：1、计算S1S2S3，2、纠错，3、去冗余，生成序

列

%hamm为七位一组的汉明序列，sourout为译码输出的最终序列

N=length(hamm);

H=[ 1 1 1 0 1 0 0;

1 1 0 1 0 1 0;

1 0 1 1 0 0 1;];%监督矩阵

Scase=[0 0 0;0 0 1; 0 1 0; 0 1 1; 1 0 0; 1 0 1;1 1 0;1 1 1;];

lala1=reshape(hamm,7,N/7);%汉明序列变为7列一行

lala=lala1';

subSSS=lala*H';%计算S1S2S3

SSS=rem(subSSS,2);%得S1S2S3序列

for i=1:N/7

for n=1:8

if SSS(i,:)==Scase(1,:)%等于000时，无错

elseif SSS(i,:)==Scase(n,:)

lala(i,9-n)=rem(lala(i,9-n)+1,2);%n=2时a0错，a0为第七个数；n=3时a1错，a1为第六个数；所以，n=m时为第（9-m）个数

end

end

end

subout=lala(:,1:4);%取前四列信息位，后三列监督位舍去

subout1=subout';

sourout=reshape(subout1,1,N/7*4);%变为一行，即译码后最终序列

3、程序运行结果

图一中，左图为未进行汉明编码序列在判决前的星座图，右图为进行汉明编码的序列在判决前的星座图。很明显，在加入相同强度噪声时，进行汉明编码的序列较稳定，受噪声影响相

值分别为0、0.1、1.0、2.0。

对较小。每图中的四幅图2

σ值很小时，两序列误比特率都很对以上两方法产生的序列进行误比特率比较，可见，在2

σ值较大时，汉明编码对减小序列误比特率方面效果明显，能够较大程度提高系统小；当2

传输性能。

四、经验与收获

（1）在循环中我们加入了显示语句，但发现显示出来的不同曲线的颜色、线型都完全相同，最后实验中我们发现，显示语句中的修饰部分参数字符可以用前面定义的字符串代替，从而，我们定义了一个多维字符串数组，使每次循环时显示的图形颜色有了区别，方便进行比较。

（2）在试验一主程序2中，编程过程初期我们的理论曲线与实验曲线差别很大，总体趋势相同但相差大约一倍。后来我们找到错误——在显示时，我们对理论曲线是取N0=0:0.1:2;，而

σ=0:0.1:2;最终画图必然是不相符的。

对实验曲线却是取2

σ再进行相应转换，解决了原先的问题。

通过上述语句，将RR转化到2

（3）在汉明编码解码程序中，起初汉明编码后信噪比竟然大大超过编码前，后来我们发现，reshape函数在改变矩阵大小时，是一列一列从上到下扫描的，如下面例子：

要想解决上述问题，可以采用对原序列先进行转置在reshape的方法。

即可实现矩阵->序列的转换。