QPSK通信系统的 MonteCarlo 仿真.doc

Monte Carlo仿真 山东大学通信原理实验 QPSK通信系统的 MonteCarlo 仿真实验报告 班级：通信一班 姓名: 荣林林 20071001089 李海强 20071201086一、实验目的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信息比特对，再用这些比特对来选择信号相位向量。加性噪声的同相分量和正交分量 ，在上面讨论过，即为零均值，方差为 的统计独立的高斯随机变量。在检测器观察到的接收信号向量 ，利用上面讨论的两种检测方法，得到判决结果，并与传输符号作比较，最后对符号差错和比特差错计数三、信道纠错编码信道解码QPSK解调AWGN信道QPSK调制 信道编码输入数据 在随机信道中，错码的出现是随机的，且错码之间是统计独立的。有高斯白噪声引起的错码就具有这种性质。所以，当信道中加性噪声为这种噪声时，就成这种信道为随机信道。由于发送码元是随机的，接收端是无法预知的，也无法识别其中有无错码，为此，可在发送端的信息码元序列中，增加一些监督码元，这些监督码元与信息码元之间有一定的关系，接收端利用这种关系，来发现和纠正其中的可能存在的错码。在信息码元中加入监督码元就称为差错控制编码，也称纠错编码。汉明编码就是一种能够纠正一位错码，且编码效率较高的线性分组码。具体的汉明编码过程主要为：以（7，4）汉明编码为例，发送端信息序列以4个一组，再加上3个监督位，构成7为发送信息。其中，a0,a1,a2为监督位，a3,a4,a5,a6为信息序列。用S1，S2，S3为三个校正因子，当S1,S2,S3=0 0 1, 错码位置a0, 0 1 0, 错码位置a1, 1 0 0, 错码位置a2, 0 1 1, 错码位置a3, 1 0 1, 错码位置a4, 1 1 0, 错码位置a5, 1 1 1, 错码位置a6, 0 0 0, 无错码。在发送端编码时，信息位a3,a4,a5,a6的值决定于输入信号，因此他们是随机的，监督位的a2,a1,a0根据信息位的取值按照监督关系（即是校验因子的值全为零），可以解得：a2=a6a5a4,a1=a6a5a3,a0=a6a4a3,给出信息位后可以直接根据上式确定出监督位。在接收端，收到每个码组先计算出对应的校验位，然后判断错码情况，并对应校正过来。四、实验过程中需要用到的一些子函数模块1、高斯白噪声子函数functiongsrv1,gsrv2=gnguass(m,sgma)if nargin=0 m=0; sgma=1;elseif nargin=1 sgma=m; m=0;endu=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 (temp0.25), % With probability 1/4, source output is 00. dsource1(i)=0; dsource2(i)=0; elseif (temp0.5), % With probability 1/4, source output is 01. dsource1(i)=0; dsource2(i)=1; elseif (temp0.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;三、实验内容及程序分析（以下程序皆以N=1000为例）(一)、未加新年到纠错编码的4PSK调制系统（1）、根据系统组成框图用MATLAB编写程序，并在每条语句后进行注释（2）、最大投影点准则进行判决 A、计算方差sgma2分别为0、0.1、0.5、1.0时的符号差错概率和比特差错概率； B、画出在在每种sgma2是，在检测器输入端1000个接收到的信号加噪声的样本 C、分别画出数据点为1000、5000、10000时Monte Carlo 仿真误比特率曲线和理论误比特率曲线，比较差别，分析数据点的数量对仿真结果的影响；（3）、将检测器判决准则改为最下距离法，比较与上面结果的区别。（4）、分析整个通信系统个组成部分的原理、分析实验结果、总结概括。相应程序如下： A、计算方差sgma2分别为0、0.1、0.5、1.0时的符号差错概率和比特差错概率；%函数分为四步第一步产生随机序列，进行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 (temp0.25), % With probability 1/4, source output is 00. dsource1(i)=0; dsource2(i)=0; elseif (temp0.5), % With probability 1/4, source output is 01. dsource1(i)=0; dsource2(i)=1; elseif (temp0.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); pbps结果统计如下：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 B、画出在在每种sgma2是，在检测器输入端1000个接收到的信号加噪声的样本函数如下：%画星座图函数sgma=input(方差=);%输入方差sgma=sqrt(sgma);%求sgmaN=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 (temp0.25), % With probability 1/4, source output is 00. dsource1(i)=0; dsource2(i)=0; elseif (temp0.5), % With probability 1/4, source output is 01. dsource1(i)=0; dsource2(i)=1; elseif (temp0.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(sgma2=);end;编程出现问题解决：语句for i=1:N , % a uniform random variable between 0 and 1中红色逗号，在修改程序时，不小心输入为中文输入形式下的：“，“，从而导致在运行时报错：” ? Undefined function or method r00 for input arguments of type double. “Sgma2=0 Sgma2=0.1 Sgma2=0.5 Sgma2=1.0 C、分别画出数据点为1000、5000、10000时Monte Carlo 仿真误比特率曲线和理论误比特率曲线，比较差别，分析数据点的数量对仿真结果的影响；% 对M=4的PSK通信系统进行蒙特卡罗仿真 %echo onSNRindB1=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,*);holdsemilogy (SNRindB1,smld_symbol_err_prb,O);semilogy (SNRindB2,theo_err_prb);gridxlabel(Eb/No in dB)ylabel(error probability)title(4PSK通信系统的蒙特卡洛仿真)gtext(（注：“”理论误码率；“*”误比特率；“o”误符号率）)左:最大投影点准则判决 右:最小距离判决N= 1000 N= 5000N= 10000 N= 1000003、实验分析：通过比较在不同长度符号数情况下分别采取最大投影点准则判决和最小距离判决法的蒙特卡罗仿真所得到的误比特率曲线和理论误比特率曲线，可以得到如下结论： （1）、最大投影点准则判决和最小距离判决法判别方法的选取对误码率的影响相似。 （2）、随着数据点的数量的增加，接收端仿真结果误码率随之降低。（二）信道纠错编码（7、4）汉明码+4PSK调制的通信系统（1）、根据系统组成框图用MATLAB编写程序，并在每条语句后进行注释（2）、比较编码信号与未编码信号在随机信道中的传输性能，分析实验结果，总结概括信道纠错编码（7、4）汉明码+4PSK调制的通信系统程序如下：1、 用到的相关子程序：（1）、source子函数产生原始序列%产生原始序列，N个function a = source(N)n=0:N-1;a=rand(1,N);for i=1:N if a(i)0.5 a(i)=0; else a(i)=1; endEnd（2） 、gnguass子函数产生两路高斯白噪声(如上)（3）、Hanming子程序，将4为信息位，加上监督位，变为7位%汉明编码,将4位变为7位function hanming,s1=hanming(N)%hanming为7/4编码后的序列，s1为原始发送序列s1=source(N);temp=zeros(1,7);hanming=zeros(1,N*7/4);j=1;%hanming编码后序列角标控制for i=1:4:Ntemp(4:7)=s1(i:i+3);% temp(1)=mod(temp(7)+temp(5)+temp(4),2);% temp(2)=mod(temp(7)+temp(6)+temp(4),2);% temp(3)=mod(temp(7)+temp(6)+temp(5),2);temp(1)=xor(xor(temp(7),temp(5),temp(4);temp(2)=xor(xor(temp(7),temp(6),temp(4);temp(3)=xor(xor(temp(7),temp(6),temp(5);hanming(j:j+6)=temp;j=j+7;end（4）、整个信道编码的子程序：functionR,err_ps,err_pb=hanm(N,snr_in_dB)%N为发送序列长度，R为复数，接收序列对应的M=N*7/4;L=M/2;Eb=1;SNR=10(snr_in_dB/10);%dB为单位的信噪比转换成以1为单位的sgma=sqrt(Eb/SNR/2);%有定义式SNR=Eb/(2*sgma2)，变换一下来计算方差x,s1=hanming(N);y=zeros(1,L); %初始化y，为符号序列，长度为L，即N/2，原始发送符号序列%下面对y赋值，每相邻两个x比特组成一个y符号for i=1:2:M %因为每两个比特位一个符号，步长为2 if x(i) x(i+1)=0 0 y(i+1)/2)=1; %00 -y0 elseif x(i) x(i+1)=0 1 y(i+1)/2)=0; %01 -y1 elseif x(i) x(i+1)=1 1 y(i+1)/2)=3; %11 -y3 elseif x(i) x(i+1)=1 0 y(i+1)/2)=2; %10 -y2 endends=zeros(2,L); %定义s为2行（符号长度），L列（符号个数）的矩阵%s是映射后发送的L个符号组成的，每个符号对应为一列for i=1:L if y(i)=0 s(:,i)=1;0; %0-1;0 对第i列赋值 elseif y(i)=1 s(:,i)=0;1; %1-0;1 elseif y(i)=2 s(:,i)=-1;0; %2-1;0 else s(:,i)=0;-1; %3-0;-1 endend%获得两路高斯信号，长度也为L（与符号个数相同）n1=zeros(1,L); %两路噪声的初始化n2=zeros(1,L);%n1,n2为两路噪声信号for i=1:L n1(i),n2(i)=gnguass(0,sgma);endn=n1;n2; %n1作为第一行，n2作为第二行，长度也都为Lr=s+n; %r为加噪声后的信号，即接收信号r1=zeros(1,L);%r1为r的第一行，即所有符号的对应的第一位r2=zeros(1,L);%r2为r的第二行，即所有符号的对应第二位for i=1:L r1(i)=r(1,i);%r1为r的第一行r2(i)=r(2,i);%r2为r的第二行endR=r1+j*r2;%可利用R的幅度和相位画星座图%最小距离法S0=1 0;S1=0 1;S2=-1 0;S3=0 -1;m00=zeros(1,L);m01=zeros(1,L);m10=zeros(1,L);m11=zeros(1,L);for i=1:L m00(i)=sqrt(r1(i)-S0(1)2+(r2(i)-S0(2)2); m01(i)=sqrt(r1(i)-S1(1)2+(r2(i)-S1(2)2); m10(i)=sqrt(r1(i)-S2(1)2+(r2(i)-S2(2)2); m11(i)=sqrt(r1(i)-S3(1)2+(r2(i)-S3(2)2);enddec=zeros(1,L);for i=1:L m_min(i)=min(m00(i) m01(i) m10(i) m11(i); if (m00(i)=m_min(i), dec(i)=0; elseif(m01(i)=m_min(i), dec(i)=1; elseif(m10(i)=m_min(i), dec(i)=2; else dec(i)=3; end;endc=zeros(1,M); %c为接收的比特序列，M个for i=1:2:M if dec(i+1)/2)=0 c(i)=0; c(i+1)=1; elseif dec(i+1)/2)=1 c(i)=0; c(i+1)=0; elseif dec(i+1)/2)=2 c(i)=1; c(i+1)=0; elseif dec(i+1)/2)=3 c(i)=1; c(i+1)=1; enderr_s=0; %符号差错计数器，置零for i=1:L if y(i)=dec(i) %发送符号等于判决后的符号 err_s=err_s; else err_s=err_s+1; endend%hanming纠错d=c; %d为纠错后序列，初始化 s=zeros(3,N/4);%s为校验因子for i=1:7:M s(1,(i+6)/7)=xor(xor(xor(c(i+6),c(i+5),c(i+4),c(i+2); s(2,(i+6)/7)=xor(xor(xor(c(i+6),c(i+5),c(i+3),c(i+1); s(3,(i+6)/7)=xor(xor(xor(c(i+6),c(i+4),c(i+3),c(i); if s(:,(i+6)/7)=0;0;1 d(i)=not(c(i); elseif s(:,(i+6)/7)=0;1;0 d(i+1)=not(c(i+1); elseif s(:,(i+6)/7)=1;0;0 d(i+2)=not(c(i+2); elseif s(:,(i+6)/7)=0;1;1 d(i+3)=not(c(i+3);elseif s(:,(i+6)/7)=1;0;1 d(i+4)=not(c(i+4); elseif s(:,(i+6)/7)=1;1;0 d(i+5)=not(c(i+5); elseif s(:,(i+6)/7)=1;1;1 d(i+6)=not(c(i+6); end ende=zeros(1,N);%纠错后的信息序列k=1;for i=1:7:M e(k)=d(i+3); e(k+1)=d(i+4); e(k+