8PSK通信系统的蒙特卡罗仿真分析-matlab源程序

f=150%抽样频率g=(sign(rand(1,150)-0.5)+1)/2%产生二进制序列sn=randn(1,50*length(g));%产生加性高斯白噪声dt=2*pi/149;t=0:dt:2*pi;si=[];co=[];%si为正交分量，co为同相分量sit=[];sqt=[];%sit为同相分量幅度，sqt为正交分量幅度sb2=[];%输入二进制序列%8PSK调制过程forn=1:3:length(g);%一次取3个二进制数ifg(n)==0&&g(n+1)==0&&g(n+2)==0%b1b2b3=000时正交分量和同相分量的幅值it=-0.383*ones(1,150);qt=-0.924*ones(1,150);b2=[zeros(1,50)zeros(1,50)zeros(1,50)]elseifg(n)==0&&g(n+1)==0&&g(n+2)==1%b1b2b3=001时it=-0.924*ones(1,150);qt=-0.383*ones(1,150);b2=[zeros(1,50)zeros(1,50)ones(1,50)]elseifg(n)==1&&g(n+1)==0&&g(n+2)==0it=0.383*ones(1,150);qt=-0.924*ones(1,150);b2=[ones(1,50)zeros(1,50)zeros(1,50)]elseifg(n)==1&&g(n+1)==0&&g(n+2)==1it=0.924*ones(1,150);qt=-0.383*ones(1,150);b2=[ones(1,50)zeros(1,50)ones(1,50)]elseifg(n)==0&&g(n+1)==1&&g(n+2)==0it=-0.383*ones(1,150);qt=0.924*ones(1,150);b2=[zeros(1,50)ones(1,50)zeros(1,50)]elseifg(n)==0&&g(n+1)==1&&g(n+2)==1it=-0.924*ones(1,150);qt=0.383*ones(1,150);b2=[zeros(1,50)ones(1,50)ones(1,50)]elseifg(n)==1&&g(n+1)==1&&g(n+2)==1it=0.924*ones(1,150);qt=0.383*ones(1,150);b2=[ones(1,50)ones(1,50)ones(1,50)]elseifg(n)==1&&g(n+1)==1&&g(n+2)==0it=0.383*ones(1,150);qt=0.924*ones(1,150);b2=[ones(1,50)ones(1,50)zeros(1,50)]endsb2=[sb2b2];c=cos(f*t);s=sin(f*t);sit=[sitit];sqt=[sqtqt];co=[coc];si=[sis];endpsk=sit.*co+sqt.*si;%调制后的8psk信号%画图figure(1);plot(sb2,'LineWidth',1.5);gridon;title('二进制序列(信源)');axis([02000-1.51.5]);set(gca,'Xtick',[150:150:1800]);xlabel('(a)t/(ts/150)');figure(2);plot(psk,'LineWidth',1.5);gridon;title('8PSK调制');axis([02000-1.51.5]);xlabel('(b)t/(Ts/150)');%8PSK的解调rpsk=psk+sn;%加入加性高斯白噪声rs=[];%rs用来存放解调后的二进制序列form=1:150:50*length(g)-150;rpsk1=rpsk(m:m+149);%取一个码元sit=rpsk1.*cos(f*t);it=cumtrapz(sit)*dt;it=it(end);%相关后得的I路电平ifit0%对得到的电平进行判决rs=[rsones(1,50)];elseifit0rs=[rszeros(1,50)];endsqt=rpsk1.*sin(f*t);qt=cumtrapz(sqt)*dt;qt=qt(end);%相关后得的Q路电平ifqt0%对得到的电平进行判决rs=[rsones(1,50)];elseifqt0rs=[rszeros(1,50)];endsb3=rpsk1.*cos(f*t-pi/4);b3=cumtrapz(sb3)*dt;b3=b3(end);sb4=rpsk1.*sin(f*t-pi/4);b4=cumtrapz(sb4)*dt;b4=b4(end);b5=abs(b3+b4);%得到b3的电平并判决ifb52rs=[rsones(1,50)];elseifb52rs=[rszeros(1,50)];endend%画图figure(3);plot(rs,'LineWidth',1.5);gridon;title('解调输出');axis([01800-1.51.5]);set(gca,'Xtick',[150:150:1800]);xlabel('t/(ts/150)');%画对比图figure(4);subplot(3,1,1);plot(sb2,'LineWidth',1.5);gridon;title('二进制序列(信源)');axis([02000-1.51.5]);set(gca,'Xtick',[150:150:1800]);xlabel('(a)t/(ts/150)');subplot(3,1,2);plot(psk,'LineWidth',1.5);gridon;title('8PSK调制');axis([02000-1.51.5]);xlabel('(b)t/(Ts/150)');subplot(3,1,3);plot(rs,'LineWidth',1.5);gridon;title('解调输出');axis([01800-1.51.5]);set(gca,'Xtick',[150:150:1800]);xlabel('t/(ts/150)');%蒙特卡罗分析EsNodb=2:0.5:12;%设置信噪比范围Es=1;No=10.^(-EsNodb/10);sigma=sqrt(No/2);%噪声功率，其值随信噪比而变error=zeros(1,length(EsNodb));%错误计数sdata=zeros(1,length(EsNodb));%进行比较判决抽样值的总的计数fori=1:length(EsNodb)error(i)=0;sdata(i)=0;whileerror(i)1000%误码数1000d=ceil(rand(1,10000)*8);%产生信源10000个,返回大于或者等于指定表达式的最小整数s=sqrt(Es)*exp(j*2*pi/8*(d-1));%复基带形式r=s+sigma(i)*(randn(1,length(d))+j*randn(1,length(d)));form=1:8rd(m,:)=abs(r-sqrt(Es)*exp(j*2*pi/8*(m-1)));%rd有m行，每行对应r与m的差值，8*10000的二维数组endform=1:length(s)dd(m)=find(rd(:,m)==min(rd(:,m)));%找到rd的m列中最小的值的行序号（与之相对的判决电平值），%dd(m)即为接收到的m值，find()函数返回的是行号ifdd(m)~=d(m)%与发送的m相比,进行误码计数error(i)=error(i)+1;endendsdata(i)=sdata(i)+10000;endendpe=error./sdata;%仿真得的误码率ps=erfc(sqrt(EsNodb)*sin(pi/8));%理论误码率figure(5);semilogy(EsNodb,pe,'b*:');holdon;semilogy(EsNodb,ps,'r-');axis([2,12,10^(-2),10^0]);gridonxlabel('Es/No(db)');ylabel('误码率');legend('仿真结果','理论计算结果');%绘眼图tt=0:150-1;%显示1个码元周期内的眼图figure(6);fork=1:9eyeprsk=rpsk(k*150:(k+1)*150-1);%rpsk为接收到的8PSK信号drawnowplot(tt,eyeprsk);holdon;axis([040-55]);end