clc;clear;
f=1*300 %采样频率
g=(sign(rand(1,300)-0.5+eps)+1)/2 %产生二进制序列
sn=randn(1,100*length(g)); %产生高斯白噪声
dt=2*pi/299;
t=0:dt:2*pi;
si=[];co=[]; %si为正交分量,co为同相分量
sit=[];sqt=[]; %sit为同相分量幅度,sqt为正交分量幅度
sb2=[]; %输入二进制序列
for n=1:3:length(g); %一次取3个二进制数
if g(n)==0 && g(n+1)==0 && g(n+2)==0
%b1b2b3=000时正交分量和同相分量的幅值
it=-0.383*ones(1,300);
qt=-0.924*ones(1,300);
b2=[zeros(1,100) zeros(1,100) zeros(1,100)]
elseif g(n)==0 && g(n+1)==0 && g(n+2)==1 %b1b2b3=001时
it=-0.924*ones(1,300);
qt=-0.383*ones(1,300);
b2=[zeros(1,100) zeros(1,100) ones(1,100)]
elseif g(n)==1 && g(n+1)==0 && g(n+2)==0
it=0.383*ones(1,300);
qt=-0.924*ones(1,300);
b2=[ones(1,100) zeros(1,100) zeros(1,100) ]
elseif g(n)==1 && g(n+1)==0 && g(n+2)==1
it=0.924*ones(1,300);
qt=-0.383*ones(1,300);
b2=[ones(1,100) zeros(1,100) ones(1,100) ]
elseif g(n)==0 && g(n+1)==1 && g(n+2)==0
it=-0.383*ones(1,300);
qt=0.924*ones(1,300);
b2=[zeros(1,100) ones(1,100) zeros(1,100) ]
elseif g(n)==0 && g(n+1)==1 && g(n+2)==1
it=-0.924*ones(1,300);
qt=0.383*ones(1,300);
b2=[zeros(1,100) ones(1,100) ones(1,100) ]
elseif g(n)==1 && g(n+1)==1 && g(n+2)==1
it=0.924*ones(1,300);
qt=0.383*ones(1,300);
b2=[ones(1,100) ones(1,100) ones(1,100) ]
elseif g(n)==1 && g(n+1)==1 && g(n+2)==0
it=0.383*ones(1,300);
qt=0.924*ones(1,300);
b2=[ones(1,100) ones(1,100) zeros(1,100) ]
end
sb2=[sb2 b2];
c=cos(f*t); s=sin(f*t);
sit=[sit it]; sqt=[sqt qt];
co=[co c]; si=[si s];
end
psk=sit.*co+sqt.*si; %调制后的8psk信号
subplot(211);plot(sb2,'LineWidth',1.5);
grid on;title('二进制序列(信源)');
axis([0 3000 -1.5 1.5]);
set(gca,'Xtick',[300,600,900,1200]);
xlabel('(a) t/(ts/300)');
subplot(212);plot(psk,'LineWidth',1.5);
grid on;title('8PSK');axis([0 3000 -1.5 1.5]);
xlabel('(b) t/(Ts/300)');
rpsk=psk+sn; %加入加性高斯白噪声
rs=[]; %rs用来存放解调后的二进制序列
for m=1:300:100*length(g)-300;
rpsk1=rpsk(m:m+299); %取一个码元
sit=rpsk1.*cos(f*t);
it=cumtrapz(sit)*dt;
it=it(end); %相关后得的I路电平
if it>0 %对得到的电平进行判决
rs=[rs ones(1,100)];
elseif it<0
rs=[rs zeros(1,100)];
end
sqt=rpsk1.*sin(f*t);
qt=cumtrapz(sqt)*dt;
qt=qt(end); %相关后得的Q路电平
if qt>0 %对得到的电平进行判决
rs=[rs ones(1,100)] ;
elseif qt<0
rs=[rs zeros(1,100)];
end
sb3=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的电平并判决
if b5<2
rs=[rs ones(1,100)] ;
elseif b5>2
rs=[rs zeros(1,100)];
end
end
figure(2);
plot(rs,'LineWidth',1.5);grid on;
title('解调输出');axis([0 3000 -1.5 1.5]);
set(gca,'Xtick',[300,600,900,1200]);
xlabel(' t/(ts/300)');
%绘眼图
% figure(3);
tt=0:2*300-1; %显示2个码元周期内的眼图
figure(2);
for k=3:40
eyepsk=rpsk(k*300+1:(k+2)*300); % rpsk为接收到的8PSK信号
drawnow
plot(tt,eyepsk);hold on;
end
%蒙特卡罗分析
EsNodb=3:0.5:10; %设置信噪比范围
Es=1;
No=10.^(-EsNodb/10);
sigma=sqrt(No/2); %噪声功率,其值随信噪比而变
error=zeros(1,length(EsNodb)); %错误计数
sdata=zeros(1,length(EsNodb)); %进行比较判决抽样值的总的计数
for i=1:length(EsNodb)
error(i)=0;
sdata(i)=0;
while error(i)<1000 %误码数<1000
d=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)));
for m = 1 : 8
rd(m,:) = abs(r-sqrt(Es)*exp(j*2*pi/8*(m-1)));
%rd有m行,每行对应r与m的差值,8*10000的二维数组
end
for m=1:length(s)
dd(m)=find(rd(:,m)==min(rd(:,m)));
%找到rd的m列中最小的值的行序号(与之相对的判决电平值),
%dd(m)即为接收到的m值,find()函数返回的是行号
if dd(m)~=d(m) %与发送的m相比,进行误码计数
error(i)=error(i)+1;
end
end
sdata(i)=sdata(i)+10000;
end
end
pe=error./sdata; %仿真得的误码率
ps=erfc(sqrt(EsNodb)*sin(pi/8)); %理论误码率
figure(3);
semilogy(EsNodb,pe,'b*:');hold on;
semilogy(EsNodb,ps,'r-');
xlabel('Es/No(db)');ylabel('误码率');
legend('仿真结果','理论计算结果');
没有合适的资源?快使用搜索试试~ 我知道了~
资源推荐
资源详情
资源评论
收起资源包目录
8psk的蒙特卡洛仿真8psk的蒙特卡洛仿真.zip (1个子文件)
matlab_8psk.m 5KB
共 1 条
- 1
资源评论
依然风yrlf
- 粉丝: 798
- 资源: 2959
下载权益
C知道特权
VIP文章
课程特权
开通VIP
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功