MPSK 通信系统的 Monte Carlo 仿真

实验原理

当发送滤波器的脉冲形状为矩形脉冲时,一组M载波相位调制信号波形的一般表示式为:
在这里插入图片描述
Es代表每个传输符号的能量。因此,数字相位调制信号在几何上可用smc和sms的二维向量来表示,即
在这里插入图片描述
在AWGN信道中,一个区间内接受到的带通信号可用向量表示为
在这里插入图片描述
最大投影点准则判决:
将接收到的信号向量r投射到M个可能的传输信号向量 之一上去,并选取对应于最大投影的向量。
最小欧式距离准则判决:
求出接收到的信号向量与M个传输向量的欧式距离,选取对应的最小欧式距离的向量,该向量对应的符号即为判决输出符号。
在这里插入图片描述
Eb代表每个传输比特的能量。QPSK误比特率与BPSK的一样。
MPSK 的误符号率(相干、最佳接收机)
在这里插入图片描述
高信噪比下,MPSK 的误比特率(采用 Gray 编码:二进制-M 进制变换)
在这里插入图片描述
(一)未加信道纠错编码的QPSK调制通信系统
在这里插入图片描述
(二)未加信道纠错编码的8PSK调制通信系统
在这里插入图片描述
实验内容

(一)未加信道纠错编码的QPSK调制通信系统
1) 最大投影点准则进行判决
a, 画出噪声方差σ2分别为 0、0.1、0.5、1.0 时,在检测器输入端1000个接收到的信号加噪声的样本(星座图);
b, 分别画出数据点为 1000、10000、100000 时的Monte Carlo仿真误比特率曲线和理论误比特率曲线,比较差别,分析数据点的数量对仿真结果的影响(横坐标snr=Eb/N0(dB),格雷映射);
2) 将检测器的判决准则改为最小距离法(星座图上符号间的距离)画出数据点为100000时的Monte Carlo仿真误比特率曲线和理论误比特率曲线,(横坐标是snr=Eb/N0,格雷映射);比较与最大投影点准则进行判决结果的区别。
(二)未加信道纠错编码的8PSK调制通信系统
检测器的判决准则选为最小距离法(星座图上符号间的距离),Gray 格雷编码,比较数据点为100000时8PSK与QPSK的Monte Carlo仿真误符号率曲线、理论误符号率曲线,比较差别(横坐标是snr=Eb/N0(dB))。

函数部分

signalSource.m

function[a,b]=signalSource(N)%生成二进制信源
    source=rand(1,N);%产生随机数
    a=zeros(1,N);b=zeros(1,N);
    for i=1:N
        if(source(i)<0.25)%当随机数<0.25时,规定为00
            a(i)=0; b(i)=0;
        elseif(source(i)<0.5)%0.25<随机数<0.5时,规定为01    
            a(i)=0;b(i)=1;
        elseif(source(i)<0.75)%0.5<随机数<0.75时,规定为10
            a(i)=1;b(i)=0;
        else%0.75<随机数<1时,规定为11
            a(i)=1;b(i)=1;
        end
    end
end

signalSource2.m

function[a,b,c]=signalSource2(N)%生成二进制信源
    source=rand(1,N);%产生随机数
    a=zeros(1,N);b=zeros(1,N);c=zeros(1,N);
    for i=1:N
        if(source(i)<0.125)%当随机数<0.125时,规定为000
            a(i)=0; b(i)=0;c(i)=0;
        elseif(source(i)<0.25)%0.125<随机数<0.25时,规定为001 
            a(i)=0;b(i)=0;c(i)=1;
        elseif(source(i)<0.375)%0.25<随机数<0.375时,规定为011 
            a(i)=0;b(i)=1;c(i)=1;
        elseif(source(i)<0.5)%0.375<随机数<0.5时,规定为010 
            a(i)=0;b(i)=1;c(i)=0;
        elseif(source(i)<0.625)%0.5<随机数<0.625时,规定为110 
            a(i)=1;b(i)=1;c(i)=0;
        elseif(source(i)<0.75)%0.625<随机数<0.75时,规定为111 
            a(i)=1;b(i)=1;c(i)=1;
        elseif(source(i)<0.875)%0.75<随机数<0.875时,规定为101 
            a(i)=1;b(i)=0;c(i)=1;
        else%0.875<随机数<1时,规定为100 
            a(i)=1;b(i)=0;c(i)=0;
        end
    end
end

orthographicProjection.m

function[sm]=orthographicProjection(a,b,N)%二进制信源转换成4PSK两路正交信号
    sm=zeros(2,N);
    for i=1:N
        if((a(i)==0)&&(b(i)==0))%根据格雷码对照进行映射
            sm(:,i)=[1;0];%00时,sm=(1 0)[1;0]中的;表示分为两行    
        elseif((a(i)==0)&&(b(i)==1))
            sm(:,i)=[0;1];%01时,sm=(0 1)    
        elseif((a(i)==1)&&(b(i)==0))
            sm(:,i)=[-1;0];%10时,sm=(-1 0)    
        else
            sm(:,i)=[0;-1];%11时,sm=(0 -1)    
        end
    end
end

orthographicProjection2.m

function[sm]=orthographicProjection2(a,b,c,N)%二进制信源转换成8PSK两路正交信号
    sm=zeros(2,N);
    for i=1:N
        if((a(i)==0)&&(b(i)==0)&&(c(i)==0))%根据格雷码对照进行映射
            sm(:,i)=[sqrt(3)/2;1/2];%000时,sm=(sqrt(3)/2 1/2)[sqrt(3)/2;1/2]中的;表示分为两行    
        elseif((a(i)==0)&&(b(i)==0)&&(c(i)==1))
            sm(:,i)=[1/2;sqrt(3)/2];%001时,sm=(1/2 sqrt(3)/2)    
        elseif((a(i)==0)&&(b(i)==1)&&(c(i)==1))
            sm(:,i)=[-1/2;sqrt(3)/2];%011时,sm=(-1/2 sqrt(3)/2) 
        elseif((a(i)==0)&&(b(i)==1)&&(c(i)==0))
            sm(:,i)=[-sqrt(3)/2;1/2];%010时,sm=(-sqrt(3)/2 1/2)
        elseif((a(i)==1)&&(b(i)==1)&&(c(i)==0))
            sm(:,i)=[-sqrt(3)/2;-1/2];%110时,sm=(-sqrt(3)/2 -1/2)
        elseif((a(i)==1)&&(b(i)==1)&&(c(i)==1))
            sm(:,i)=[-1/2;-sqrt(3)/2];%111时,sm=(-1/2 -sqrt(3)/2)
        elseif((a(i)==1)&&(b(i)==0)&&(c(i)==1))
            sm(:,i)=[1/2;-sqrt(3)/2];%101时,sm=(1/2 -sqrt(3)/2)
        else
            sm(:,i)=[sqrt(3)/2;-1/2];%100时,sm=(sqrt(3)/2 -1/2)    
        end
    end
end

guass.m

function[n]=guass(N,sigma)%生成两路正交高斯噪声
    nc=zeros(1,N);    
    ns=zeros(1,N); 
    for i=1:N
        u=rand;
        z=sigma*sqrt(2*log10(1/(1-u)));         
        u=rand;        
        gsrv1=z*cos(2*pi*u);         
        gsrv2=z*sin(2*pi*u);         
        nc(i)=gsrv1;        
        ns(i)=gsrv2; 
    end
    n=zeros(2,N); %写成矩阵模式便于和之前求得的sm相加 
    n(1,:)=nc; %(1,:)表示第一行的元素
    n(2,:)=ns; 
end

maximalProjection.m

function[c]=maximalProjection(r,N)%最大投影准则
    d=zeros(1,4);
    c=zeros(1,N);
    for i=1:N
        d(1)=1*r(1,i)+0*r(2,i);
        d(2)=0*r(1,i)+1*r(2,i);
        d(3)=(-1)*r(1,i)+0*r(2,i);
        d(4)=0*r(1,i)+(-1)*r(2,i);%将信道输出的r和原信号s逐个作向量积
        dm=d(1);
        for k=2:4 %通过for循环逐个比较求出最大向量积
            if dm<d(k)  
                dm=d(k);
            end
        end
        if dm==d(1) %根据映射的逆过程进行判决,恢复原信号
            c(i)=0;
        elseif dm==d(2)
            c(i)=1;
        elseif dm==d(3)
            c(i)=2;
        elseif dm==d(4)
            c(i)=3;
        end
    end
end

minimumDistance.m

function[c]=minimumDistance(r,N)%最小欧氏距离准则
    d=zeros(1,4);
    c=zeros(1,N);
    for i=1:N
        d(1)=(r(1,i)-1)^2+(r(2,i)-0)^2;
        d(2)=(r(1,i)-0)^2+(r(2,i)-1)^2;
        d(3)=(r(1,i)-(-1))^2+(r(2,i)-0)^2;
        d(4)=(r(1,i)-0)^2+(r(2,i)-(-1))^2;
        dm=d(1);
        for k=2:4 %通过for循环逐个比较求出最小距离
            if dm>d(k)
                dm=d(k);
            end
        end
        if dm==d(1) %根据映射的逆过程进行判决,恢复原信号
            c(i)=0;
        elseif dm==d(2)
            c(i)=1;
        elseif dm==d(3)
            c(i)=2;
        elseif dm==d(4)
            c(i)=3;
        end
    end
end

minimumDistance2.m

function[d]=minimumDistance2(r,N)%最小欧氏距离准则
    e=zeros(1,8);
    d=zeros(1,N);
    for i=1:N
        e(1)=(r(1,i)-sqrt(3)/2)^2+(r(2,i)-1/2)^2;
        e(2)=(r(1,i)-1/2)^2+(r(2,i)-sqrt(3)/2)^2;
        e(3)=(r(1,i)-(-1/2))^2+(r(2,i)-sqrt(3)/2)^2;
        e(4)=(r(1,i)-(-sqrt(3)/2))^2+(r(2,i)-1/2)^2;
        e(5)=(r(1,i)-(-sqrt(3)/2))^2+(r(2,i)-(-1/2))^2;
        e(6)=(r(1,i)-(-1/2))^2+(r(2,i)-(-sqrt(3)/2))^2;
        e(7)=(r(1,i)-1/2)^2+(r(2,i)-(-sqrt(3)/2))^2;
        e(8)=(r(1,i)-sqrt(3)/2)^2+(r(2,i)-(-1/2))^2;
        dm=e(1);
        for k=2:8 %通过for循环逐个比较求出最小距离
            if dm>e(k)
                dm=e(k);
            end
        end
        if dm==e(1) %根据映射的逆过程进行判决,恢复原信号
            d(i)=0;
        elseif dm==e(2)
            d(i)=1;
        elseif dm==e(3)
            d(i)=2;
        elseif dm==e(4)
            d(i)=3;
        elseif dm==e(5)
            d(i)=4;
        elseif dm==e(6)
            d(i)=5;
        elseif dm==e(7)
            d(i)=6;
        elseif dm==e(8)
            d(i)=7;
        end
    end
end

rebuild.m

function[y]=rebuild(c,N)%QPSK重新建立原信号
    y=zeros(1,2*N);
    for i=1:N
        if c(i)==0
            y(2*i-1)=0;y(2*i)=0;  
        elseif c(i)==1
            y(2*i-1)=0;y(2*i)=1;   
        elseif c(i)==2
            y(2*i-1)=1;y(2*i)=0;   
        elseif c(i)==3
            y(2*i-1)=1;y(2*i)=1;
        end
    end
end

rebuild2.m

function[y]=rebuild2(d,N)%8PSK重新建立原信号
    y=zeros(1,3*N);
    for i=1:N
        if d(i)==0
            y(3*i-2)=0;y(3*i-1)=0;y(3*i)=0;  
        elseif d(i)==1
            y(3*i-2)=0;y(3*i-1)=0;y(3*i)=1;   
        elseif d(i)==2
            y(3*i-2)=0;y(3*i-1)=1;y(3*i)=1;   
        elseif d(i)==3
            y(3*i-2)=0;y(3*i-1)=1;y(3*i)=0;
        elseif d(i)==4
            y(3*i-2)=1;y(3*i-1)=1;y(3*i)=0;   
        elseif d(i)==5
            y(3*i-2)=1;y(3*i-1)=1;y(3*i)=1;   
        elseif d(i)==6
            y(3*i-2)=1;y(3*i-1)=0;y(3*i)=1;
        elseif d(i)==7
            y(3*i-2)=1;y(3*i-1)=0;y(3*i)=0;      
        end
    end
end

error1.m

function[pb,ps]=error1(y,a,b,N)%QPSK计算误比特率和误码率
    bitNum=0;symbolNum=0;
    for i=1:N  
        symbol=0;
        if(y(2*i-1)~=a(i))     
            bitNum=bitNum+1;    
            symbol=1;
        end
        if(y(2*i)~=b(i))
            bitNum=bitNum+1;
            symbol=1;
        end
        if(symbol==1)     
            symbolNum=symbolNum+1;
        end
    end 
    pb=bitNum/(2*N);ps=symbolNum/N;
end

error2.m

function[pb,ps]=error2(y,a,b,c,N)% 8PSK计算误比特率和误码率
    bitNum=0;symbolNum=0;
    for i=1:N  
        symbol=0;
        if(y(3*i-2)~=a(i))     
            bitNum=bitNum+1;    
            symbol=1;
        end
        if(y(3*i-1)~=b(i))
            bitNum=bitNum+1;
            symbol=1;
        end
        if(y(3*i)~=c(i))
            bitNum=bitNum+1;
            symbol=1;
        end
        if(symbol==1)     
            symbolNum=symbolNum+1;
        end
    end 
    pb=bitNum/(3*N);ps=symbolNum/N;
end

errorRateCurve.m

function[pb,ps]=errorRateCurve(N,a,b,sigma)%QPSK误比特率/误码率曲线    
    sm=orthographicProjection(a,b,N);
    n=guass(N,sigma); 
    r=sm+n; 
    c=maximalProjection(r,N);   %c=minimumDistance(r,N);  
    y=rebuild(c,N);  
    [pb,ps]=error1(y,a,b,N);
end

errorRateCurve2.m

function[pb,ps]=errorRateCurve2(N,a,b,c,sigma)%8PSK误比特率/误码率曲线    
    sm=orthographicProjection2(a,b,c,N);
    n=guass(N,sigma); 
    r=sm+n;    
    d=minimumDistance2(r,N);  
    y=rebuild2(d,N);  
    [pb,ps]=error2(y,a,b,c,N);
end

主程序

(一)
1) a,
planisphere.m

N=input('N=');%输入信源长度
s=input('方差=');%输入噪声方差
sigma=sqrt(s);%计算标准差
[a,b]=signalSource(N);%信源信号
sm=orthographicProjection(a,b,N);%将信源信号映射成4PSK两路正交信号 
n=guass(N,sigma);%产生噪声
r=sm+n;%加入噪声
c=maximalProjection(r,N);%利用最大投影准则判决
for i=1:N
    if(c(i)==0)
        plot(r(1,i),r(2,i),'B*');         
        hold on;    
    elseif(c(i)==1)
        plot(r(1,i),r(2,i),'R*');         
        hold on;            
    elseif(c(i)==2)
        plot(r(1,i),r(2,i),'Y*');          
        hold on;              
    elseif(c(i)==3)
        plot(r(1,i),r(2,i),'G*');           
        hold on;
    end
end
axis([-2 2 -2 2]);
line([2,-2],[0,0],'linewidth',1,'color','black')
line([0,0],[2,-2],'linewidth',1,'color','black')
title('星座图');
hold off

b,
bitErrorRateCurve.m

N=input('N=');%输入信源长度
snr=0:0.5:10;%信噪比范围 
Eb=1;  
sigma=zeros(1,length(snr));  %产生空序列用来定义仿真图的sigma 
pb1=zeros(1,length(snr));    %产生空序列用来定义仿真图的未加纠错码的仿真误比特率
pb2=zeros(1,length(snr));    %产生空序列用来定义仿真图的未加纠错码的理论误比特率
for i=1:length(snr)    %求不同i代表的信噪比对应的误比特率         
    [a,b]=signalSource(N);   %生成要发送序列             
    sigma(i)=sqrt((Eb/(10^(snr(i)/10)))/2);  %由信噪比求噪声标准差       
    [pb,ps]=errorRateCurve(N,a,b,sigma(i));  %求未加纠错码误比特率 
    pb1(i)=pb;
end  
semilogy(snr,pb1,'r'); %画出未加纠错码的仿真误比特率曲线 ,仿真为红色
hold on; 
for i=1:length(snr)%计算信噪比区间大小
    EN=(10^(snr(i)/10));%即Eb/N0
    pb2(i)= qfunc(sqrt(2*EN));%Q函数
end
semilogy(snr,pb2,'b'); %画出未加纠错码的理论误比特率曲线 ,仿真为蓝色
grid on
xlabel('信噪比SNR/dB')
ylabel('误比特率') 
title('QPSK通信系统的Monte Carlo仿真')
legend('仿真误比特率','理论误比特率')
hold off;

2)
代码与1)b类似,只是将errorRateCurve函数里的maximalProjection()改为minimumDistance()

(二)
symbolErrorRateCurve1.m

N=input('N=');%输入信源长度
snr=0:0.5:10;%信噪比范围 
Eb=1;
M=4;
sigma=zeros(1,length(snr));  %产生空序列用来定义仿真图的sigma 
ps1=zeros(1,length(snr));    %产生空序列用来定义仿真图的未加纠错码的仿真误符号率
ps2=zeros(1,length(snr));    %产生空序列用来定义仿真图的未加纠错码的理论误符号率
for i=1:length(snr)    %求不同i代表的信噪比对应的误符号率         
    [a,b]=signalSource(N);   %生成要发送序列             
    sigma(i)=sqrt((Eb/(10^(snr(i)/10)))/2);  %由信噪比求噪声标准差       
    [pb,ps]=errorRateCurve(N,a,b,sigma(i));  %求未加纠错码误符号率 
    ps1(i)=ps;
end  
semilogy(snr,ps1,'r'); %画出未加纠错码的仿真误符号率曲线 ,仿真为红色
hold on; 
for i=1:length(snr)%计算信噪比区间大小
    EN=(10^(snr(i)/10));%即Eb/N0
    ps2(i)=2*qfunc(sqrt(2*log2(M)*EN)*sin(pi/M));%Q函数
end
semilogy(snr,ps2,'b'); %画出未加纠错码的理论误符号率曲线 ,仿真为蓝色
grid on
xlabel('信噪比SNR/dB')
ylabel('误符号率') 
title('QPSK通信系统的Monte Carlo仿真')
legend('仿真误符号率','理论误符号率')
hold off;

symbolErrorRateCurve2.m

N=input('N=');%输入信源长度
snr=0:0.5:10;%信噪比范围 
Eb=1;
M=8;
sigma=zeros(1,length(snr));  %产生空序列用来定义仿真图的sigma 
ps1=zeros(1,length(snr));    %产生空序列用来定义仿真图的未加纠错码的仿真误符号率
ps2=zeros(1,length(snr));    %产生空序列用来定义仿真图的未加纠错码的理论误符号率
for i=1:length(snr)    %求不同i代表的信噪比对应的误符号率         
    [a,b,c]=signalSource2(N);   %生成要发送序列             
    sigma(i)=sqrt((Eb/(10^(snr(i)/10)))/2);  %由信噪比求噪声标准差       
    [pb,ps]=errorRateCurve2(N,a,b,c,sigma(i));  %求未加纠错码误符号率 
    ps1(i)=ps;
end  
semilogy(snr,ps1,'r'); %画出未加纠错码的仿真误符号率曲线 ,仿真为红色
hold on; 
for i=1:length(snr)%计算信噪比区间大小
    EN=(10^(snr(i)/10));%即Eb/N0
    ps2(i)=2*qfunc(sqrt(2*log2(M)*EN)*sin(pi/M));%Q函数
end
semilogy(snr,ps2,'b'); %画出未加纠错码的理论误符号率曲线 ,仿真为蓝色
grid on
xlabel('信噪比SNR/dB')
ylabel('误符号率') 
title('8PSK通信系统的Monte Carlo仿真')
legend('仿真误符号率','理论误符号率')
hold off;

实验结果

(一)
1)
a,σ^2为0:
在这里插入图片描述
σ^2为0.1:
在这里插入图片描述
σ^2为0.5:
在这里插入图片描述
σ^2为1.0:
在这里插入图片描述
随着噪声方差的增大,所得信号的离散程度增加,说明噪声对原信号的影响增加

b,snr以dB为单位
在这里插入图片描述
数据点为 1000:
在这里插入图片描述
数据点为 10000:
在这里插入图片描述
数据点为 100000:
在这里插入图片描述
随着数据点增加,仿真曲线和理论曲线越来越吻合。

  1. 数据点为 100000:
    在这里插入图片描述
    两种判决准则结果区别不大

(二)
数据点为100000:
在这里插入图片描述
在这里插入图片描述
二者在低SNR时区别不大,但是高SNR时,QPSK可以达到非常低的误符号率,而8PSK则不行,所以在高SNR时QPSK优于8PSK。

  • 2
    点赞
  • 22
    收藏
    觉得还不错? 一键收藏
  • 5
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 5
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值