通信原理实验二QPSK通信系统的-Monte-Carlo-仿真.doc

通信原理实验二

QPSK通信系统的 MonteCarlo 仿真

一、实验目的
1、提高独立学习的能力；
2、培养发现问题、解决问题和分析问题的能力；
3、学习Matlab 的使用；
4、掌握4PSK通信系统的Monte Carlo仿真方法；
5、掌握4PSK通信系统的的组成原理；
6、比较编码信号和为编码信号在随机信道中的传输，加深对纠错编码的理解；

二、系统框图及编程原理
实验原理
PSK是利用载波的不同相位表示相应的数字信息。对于二进制相位调制（M=2)来说，两个载波相位是0和π。对于M相相位调制来说M=2k，这里k是每个传输符号的信息比特数。4PSK是M=4的载波相位调制。
这里，将理论差错概率与仿真的差错概率比较，进一步观察仿真与理论值之间的差别。同时，用不同的判决准则对接受信号进行判决。并比较两种判别方法的差别。
一．QPSK调制原理
1．信号能量分析
一组M载波相位调制信号波形的一般表示式为
m=[0，M-1]
式中是发送滤波器的脉冲形状，它决定了传输信号的频谱特性，A是信号的幅度。注意到，PSK信号对所有m都具有相等的能量，即

代表每个传输符号的能量。
2．噪声分析
传输信号的信道假设被加性噪声n(t)所污损，这样信号在接收端将产生误码。因为n(t)是功率谱为的白高斯过程的一个样本函数，所以噪声分量就是零均值高斯型的，即

3．信号判决分析
最佳检测器将接收信号向量r投射到M个可能的传输信号向量{}之一上去，并选取对应于最大投影的向量。据此，得到相关准则为
m=[0,M-1]
检测器观察到接收信号向量，并计算r在4种可能的信号向量上的投影。根据选取对应于最大投影的信号点作为判决，从而判决出信号。
同时，检测器的判决准则也可采用最小距离法，即利用星座图上符号间的距离进行判决，从而得到判决结果。
二．Monte Carlo 仿真过程
仿真框图如图（一）

图（一） 用于 Monte Carlo仿真的4PSK系统的方框图

如图所示，利用一个随机数发生器，产生（0，1）范围内的随机数。再将这个范围分成四个相等的区间（0，0.25），（0.25，0.5），（0.5，0.75），（0.75，1.0），这些子区间分别对应于00，01，11，10信息比特对，再用这些比特对来选择信号相位向量。加性噪声的同相分量和正交分量 ，在上面讨论过，即为零均值，方差为 的统计独立的高斯随机变量。在检测器观察到的接收信号向量 ，利用上面讨论的两种检测方法，得到判决结果，并与传输符号作比较，最后对符号差错和比特差错计数
三、实验内容及程序分析（以下程序皆以N=1000为例）
%映射比较子函数
%函数分为四步第一步产生随机序列，进行4PSK映射。
%第二步：调用高斯高斯噪声子函数，产生正交两路高斯噪声，与输出符号序列相加
%第三步：检测接受信号
%第四步：计算误码率和误比特率

N=1000; %符号长度
E=1; %计算噪声方差
sgma=input('方差=');
sgma=sqrt(sgma);
% 4PSK比特映射
s00=[1 0];
s01=[0 1];
s11=[-1 0];
s10=[0 -1];
%第一步产生随机序列，进行4PSK映射。生成随机信源
for i=1:N, %生成随机信源
temp=rand;
if (temp<0.25), % With probability 1/4, source output is "00."
dsource1(i)=0;
dsource2(i)=0;
elseif (temp<0.5), % With probability 1/4, source output is "01."
dsource1(i)=0;
dsource2(i)=1;
elseif (temp<0.75), % With probability 1/4, source output is "10."
dsource1(i)=1;
dsource2(i)=0;
else % With probability 1/4, source output is "11."
dsource1(i)=1;
dsource2(i)=1;
end;
end;
numofsymbolerror=0; %检测错误并计算错误率
numofbiterror=0;
%第二步：调用高斯高斯噪声子函数，产生正交两路高斯噪声，与输出符号序列相加
for i=1:N,
[gsrv1,gsrv2]=gnguass(0,sgma); %调用高斯噪声子函数
n(1)=gsrv1;
n(2)=gsrv2;
if ((dsource1(i)==0) & (dsource2(i)==0)),
r=s00+n;
elseif ((dsource1(i)==0) & (dsource2(i)==1)),
r=s01+n;
elseif ((dsource1(i)==1) & (dsource2(i)==0)),
r=s10+n;
else
r=s11+n;
end;
%第三步：检测接受信号
c00=sqrt((r(1)-s00(1))^2+(r(2)-s00(2))^2); %最小距离判决
c01=sqrt((r(1)-s01(1))^2+(r(2)-s01(2))^2);
c10=sqrt((r(1)-s10(1))^2+(r(2)-s10(2))^2);
c11=sqrt((r(1)-s11(1))^2+(r(2)-s11(2))^2);
c_min=min([c00 c01 c10 c11]);
if (c00==c_min),
decis1=0; decis2=0;
elseif (c01==c_min),
decis1=0; decis2=1;
elseif (c10==c_min),
decis1=1; decis2=0;
else
decis1=1; decis2=1;
end;

%第四步：计算误码率和误比特率
symbolerror=0; %设置符号错误标志，以统计错误个数
if (decis1~=dsource1(i)),
numofbiterror=numofbiterror+1;
symbolerror=1;
end;
if (decis2~=dsource2(i)),
numofbiterror=numofbiterror+1;
symbolerror=1;
end;
if (symbolerror==1),
numofsymbolerror = numofsymbolerror+1;
end;
end;
ps=numofsymbolerror/N; % 错误率计算
pb=numofbiterror/(2*N);
pb
ps
结果统计如下：
sgma=0
最大投影准则: pb =0.0 ps=0.0
最小距离判决: pb =0.0 ps=0.0
sgma=0.1
最大投影准则: pb =5.0000e-004 ps=1.0000e-003
最小距离判决: pb =1.0000e-003 ps=0.0020
sgma=0.5
最大投影准则: pb =0.0690 ps=0.1320
最小距离判决: pb =0.0690 ps=0.1360

sgma=1.0
最大投影准则: pb =0.1345 ps=0.2450
最小距离判决: pb =0.1455 ps=0.2680
>> %j化星座图函数
sgma=input('方差=');%输入方差
sgma=sqrt(sgma);%求sgma
N=1000;%输入数据点数
E=1;% 单个符号能量

%4PSK映射规则%
s00=[1 0];
s01=[0 1];
s11=[-1 0];
s10=[0 -1];
% 化星座图
%产生信源序列%
for i=1:N,% a uniform random variable between 0 and 1
temp=rand;%随机序列
if (temp<0.25), % With probability 1/4, source output is "00."
dsource1(i)=0;
dsource2(i)=0;
elseif (temp<0.5), % With probability 1/4, source output is "01."
dsource1(i)=0;
dsource2(i)=1;
elseif (temp<0.75), % With probability 1/4, source output is "10."
dsource1(i)=1;
dsource2(i)=0;
else % With probability 1/4, source output is "11."
dsource1(i)=1;
dsource2(i)=1;
end;
end;
for i=1:N,
% The received signal at the detector, for the ith symbol, is:
[gsrv1,gsrv2]=gnguass(0,sgma);
n(1)=gsrv1;
n(2)=gsrv2;
if ((dsource1(i)==0) & (dsource2(i)==0)),
r00=s00+n;
elseif ((dsource1(i)==0) & (dsource2(i)==1)),
r01=s01+n;
elseif ((dsource1(i)==1) & (dsource2(i)==0)),
r11=s10+n;
else
r10=s11+n;
end;
hold on;
plot(r00(1),r00(2),'*g',r01(1),r01(2),'*r',r11(1),r11(2),'*y',r10(1),r10(2),'*b');
title('sgma^2=');
end;

主函数
% 对M=4的PSK通信系统进行蒙特卡罗仿真 %
echo on
SNRindB1=0:2:10;
SNRindB2=0:0.1:10;
for i=1:length (SNRindB1),%计算信噪比区间大小
[pb , ps]=cm_sm32(SNRindB1 (i)); %仿真比特和符号误码率
smld_bit_err_prb (i)=pb;%将误比特率赋值
smld_symbol_err_prb (i)=ps;%将误码率赋值
end;
for i=1:length (SNRindB2),%计算信噪比区间大小
SNR=exp(SNRindB2(i)*log(10)/10); %信噪比
theo_err_prb(i)=erfc(sqrt(2*SNR)); %理论比特误码率
end;
%随后绘图曲线
semilogy (SNRindB1,smld_bit_err_prb,'*');
hold
semilogy (SNRindB1,smld_symbol_err_prb,'O');
semilogy (SNRindB2,theo_err_prb);
grid
xlabel('Eb/No in dB')
ylabel('error probability')
title('4PSK通信系统的蒙特卡洛仿真')
gtext('（注：“—”理论误码率；“*”误比特率；“o”误符号率）')
左:最大投影点准则判决 右:最小距离判决
N= 1000

N= 5000

N= 10000

%函数分为七步
%第一步：产生随机序列
%第二步：实现4码到7码的变换
%第三步：进行4PSK映射。
%第四步：调用高斯高斯噪声子函数，产生正交两路高斯噪声，与输出符号序列相加
%第五步：7码到4码的解码
%第六步：进行码元检测%
%第七步：计算误码率和误比特率
sgma=input('噪声方差=');%输入噪声方差
sgma=sqrt(sgma);%求sgma
Eb=1; %每比特能量
N=input('数据点数=');%输入数据点数
%第一步产生随机序列
for i=1:N
temp=rand;%随机序列
if (temp>0.5)
source(i)=1;
else
source(i)=0;
end
end
%(7,4)汉明编码
%第二步：实现4码到7码的变换
for k=1:N/4
out(7*k-6)=source(4*k-3);
out(7*k-5)=source(4*k-2);
out(7*k-4)=source(4*k-1);
out(7*k-3)=source(4*k);
out(7*k-2)=xor(xor(out(7*k-6),out(7*k-5)),out(7*k-4));
out(7*k-1)=xor(xor(out(7*k-6),out(7*k-5)),out(7*k-3));
out(7*k)=xor(xor(out(7*k-6),out(7*k-4)),out(7*k-3));
end
%第三步：进行4PSK映射。
s00=[1 0];%s00=[1 0]
s01=[0 1];%s01=[0 1]
s11=[-1 0];%s11=[-1 0]
s10=[0 -1];%s10=[0 -1]
%进行4PSK映射%
for i=1:N/2
if(source(2*i-1)==0&&source(2*i)==0)
s=s00;
elseif(source(2*i-1)==0&&source(2*i)==1)
s=s01;
elseif(source(2*i-1)==1&&source(2*i)==0)
s=s10;
elseif(source(2*i-1)==1&&source(2*i)==1)
s=s11;
end
%第四步：调用高斯高斯噪声子函数，产生正交两路高斯噪声，与输出符号序列相加[n(1),n(2)]=gnguass(0,sgma);%调用gnguass函数
r=s+n;
%(7,4)汉明解码%
%第五步：7码到4码的解码
for k=1:N/4
jiema(4*k-3)=out(7*k-6);
jiema(4*k-2)=out(7*k-5);
jiema(4*k-1)=out(7*k-4);
jiema(4*k)=out(7*k-3);
end
%第六步：进行码元检测%
%最大投影点准则%
c00=dot(r,s00);%取r在s00上的投影
c01=dot(r,s01);%取r在s01上的投影
c10=dot(r,s10);%取r在s10上的投影
c11=dot(r,s11);%取r在s11上的投影
c_max=max([c00 c01 c10 c11]);%取c00，c01，c10，c11中的最大值
if (c00==c_max),
decis1=0;decis2=0;
elseif(c01==c_max),
decis1=0;decis2=1;
elseif(c10==c_max),
decis1=1;decis2=0;
else
decis1=1;decis2=1;
end
out(2*i-1)=decis1;
out(2*i)=decis2;
end
%第七步：计算误码率和误比特率
%计算符号差错概率%
symbolerror=0;
for i=1:N/2
if(out(2*i-1)~=source(2*i-1)||out(2*i)~=source(2*i))
symbolerror=symbolerror+1;
end
end
ps=2*symbolerror/N;
%计算比特差错概率
biterror=0;
for i=1:N
if(out(i)~=source(i))
biterror=biterror+1;
end
end
pb=biterror/N;

1、高斯白噪声子函数
function[gsrv1,gsrv2]=gnguass(m,sgma)
if nargin==0
m=0;
sgma=1;
elseif nargin==1
sgma=m;
m=0;
end
u=rand;
z=sgma*sqrt(2*log10(1/(1-u)));
u=rand;
gsrv1=m+z*cos(2*pi*u);
gsrv2=m+z*sin(2*pi*u);
函数评注：编程产生正交两路高斯白噪声
2、映射检测比较子函数
（1）、采用最大投影准则
function [pb,ps]=cm_sm32(snr_in_dB)
% [pb,ps]=cm_sm32(snr_in_dB)
% pb误符号率
% ps误比特率
N=1000; %符号长度
E=1; %计算噪声方差
snr=10^(snr_in_dB/10);
sgma=sqrt(E/snr/2);
s00=[1 0]; % 比特映射
s01=[0 1];
s11=[-1 0];
s10=[0 -1];
for i=1:N, %生成随机信源
temp=rand;
if (temp<0.25), % With probability 1/4, source output is "00."
dsource1(i)=0;
dsource2(i)=0;
elseif (temp<0.5), % With probability 1/4, source output is "01."
dsource1(i)=0;
dsource2(i)=1;
elseif (temp<0.75), % With probability 1/4, source output is "10."
dsource1(i)=1;
dsource2(i)=0;
else % With probability 1/4, source output is "11."
dsource1(i)=1;
dsource2(i)=1;
end;
end;
numofsymbolerror=0; %检测错误并计算错误率
numofbiterror=0;
for i=1:N,
[gsrv1,gsrv2]=gnguass(0,sgma); %调用高斯噪声子函数
n(1)=gsrv1;
n(2)=gsrv2;
if ((dsource1(i)==0) & (dsource2(i)==0)),
r=s00+n;
elseif ((dsource1(i)==0) & (dsource2(i)==1)),
r=s01+n;
elseif ((dsource1(i)==1) & (dsource2(i)==0)),
r=s10+n;
else
r=s11+n;
end;
c00=dot(r,s00); %最大投影点准则，即向量点乘积
c01=dot(r,s01);
c10=dot(r,s10);
c11=dot(r,s11);
c_max=max([c00 c01 c10 c11]);
if (c00==c_max),
decis1=0; decis2=0;
elseif (c01==c_max),
decis1=0; decis2=1;
elseif (c10==c_max),
decis1=1; decis2=0;
else
decis1=1; decis2=1;
end;
symbolerror=0; %设置符号错误标志，以统计错误个数
if (decis1~=dsource1(i)),
numofbiterror=numofbiterror+1;
symbolerror=1;
end;
if (decis2~=dsource2(i)),
numofbiterror=numofbiterror+1;
symbolerror=1;
end;
if (symbolerror==1),
numofsymbolerror = numofsymbolerror+1;
end;
end;
ps=numofsymbolerror/N; % 错误率计算
pb=numofbiterror/(2*N);
（2）最小距离判别法（只需将上面红色部分改为下面程序即可）
c00=sqrt((r(1)-s00(1))^2+(r(2)-s00(2))^2); %最小距离判决
c01=sqrt((r(1)-s01(1))^2+(r(2)-s01(2))^2);
c10=sqrt((r(1)-s10(1))^2+(r(2)-s10(2))^2);
c11=sqrt((r(1)-s11(1))^2+(r(2)-s11(2))^2);
c_min=min([c00 c01 c10 c11]);
if (c00==c_min),
decis1=0; decis2=0;
elseif (c01==c_min),
decis1=0; decis2=1;
elseif (c10==c_min),
decis1=1; decis2=0;
else
decis1=1; decis2=1;
end;

24