clc
clear all
%首先我们将QPSK在星座图上的位置分别用s00,s01,s10,s11表示
s00=1/sqrt(2)*(-1-j);
s01=1/sqrt(2)*(-1+j);
s10=1/sqrt(2)*(1-j);
s11=1/sqrt(2)*(1+j);
ber=zeros(10,15);%对误码率进行初始化
for t=1:10;%蒙特卡洛循环次数为10
%产生信源s
num=10000;%信号点数为1000000
s=randint(1,num);%用randint函数产生01比特流
x=[];
k=1;
%将每两个比特一组对01比特流进行QPSK编码，信息存储在x中
for i=1:2:length(s)-1
if(s(i)==0&&s(i+1)==0)
x(k)=s00;
k=k+1;
elseif(s(i)==0&&s(i+1)==1)
x(k)=s01;
k=k+1;
elseif(s(i)==1&&s(i+1)==0)
x(k)=s10;
k=k+1;
elseif (s(i)==1&&s(i+1)==1)
x(k)=s11;
k=k+1;
end
end
%将信号通过高斯信道，信噪比为0：14dB
snr=0:14;
for i=1:length(snr)
N(i)=power(10,-snr(i)/10);%计算噪声功率
n=sqrt(N(i)/2).*randn(1,length(x))+j.*sqrt(N(i)/2).*randn(1,length(x));
r=x+n;
%r=awgn(x,snr(i),'measured');
%采用最大似然检测算法对接收信号r进行判决
for k=1:length(x)
d00=abs(r(k)-s00);
d01=abs(r(k)-s01);
d10=abs(r(k)-s10);
d11=abs(r(k)-s11);
if((d00<d01)&(d00<d10)&(d00<d11))
y(2*k-1)=0;
y(2*k)=0;
elseif((d01<d00)&(d01<d10)&(d01<d11))
y(2*k-1)=0;
y(2*k)=1;
elseif((d10<d00)&(d10<d01)&(d10<d11))
y(2*k-1)=1;
y(2*k)=0;
else
y(2*k-1)=1;
y(2*k)=1;
end
end
%计算仿真误码率
for k=1:length(s)
if(s(k)~=y(k)) ;
ber(t,i)=ber(t,i)+1;
end
end
ber(t,i)=ber(t,i)/num;
end
end
bern=sum(ber);
bern=bern/10;
snrn=power(10,snr/10);
q=qfunc(sqrt(snrn));%计算误码率理论值
%画图
semilogy(snr,bern,'-',snr,q,'+')
grid on
xlabel('Es/No(dB)');
ylabel('误码率');
title('QPSK调制信号在AWGN下的蒙特卡罗仿真');
legend('仿真误码率','理论误码率');
grid on
