clc
clear all
snr=-10:1:35; %snr是信噪比
PsPn=10.^(0.1*snr); %Ps/Pn是信号功率比噪声功率，也是信号和噪声的能量比
Capacity=[]; %信道容量
log2py=[]; %信道容量公式第一项
signal_space=[]; %信号空间
sample=66666; %样本个数
%QPSK 8PSK 16PSK调制的信道容量曲线
for k=2:4
M=2.^k;
if k==2 %QPSK信号星座图
signal_space(1) = complex(1,1);
signal_space(2) = complex(-1,1);
signal_space(3) = complex(-1,-1);
signal_space(4) = complex(1,-1);
end
if k==3 %8PSK信号星座图
signal_space(1) = complex(1,0);
signal_space(2) = complex(cos(pi/4),sin(pi/4));
signal_space(3) = complex(0,1);
signal_space(4) = complex(cos(3*pi/4),sin(3*pi/4));
signal_space(5) = complex(-1,0);
signal_space(6) = complex(cos(5*pi/4),sin(5*pi/4));
signal_space(7) = complex(0,-1);
signal_space(8) = complex(cos(-pi/4),sin(-pi/4));
end
if k==4 %16PSK信号星座图
signal_space(1) = complex(1,0);
signal_space(2) = complex(cos(pi/8),sin(pi/8));
signal_space(3) = complex(cos(pi/4),sin(pi/4));
signal_space(4) = complex(cos(3*pi/8),sin(3*pi/8));
signal_space(5) = complex(0,1);
signal_space(6) = complex(cos(5*pi/8),sin(5*pi/8));
signal_space(7) = complex(cos(3*pi/4),sin(3*pi/4));
signal_space(8) = complex(cos(7*pi/8),sin(7*pi/8));
signal_space(9) = complex(-1,0);
signal_space(10) = complex(cos(7*pi/8),sin(-7*pi/8));
signal_space(11) = complex(cos(3*pi/4),sin(-3*pi/4));
signal_space(12) = complex(cos(5*pi/8),sin(-5*pi/8));
signal_space(13) = complex(0,-1);
signal_space(14) = complex(cos(3*pi/8),sin(-3*pi/8));
signal_space(15) = complex(cos(pi/4),sin(-pi/4));
signal_space(16) = complex(cos(pi/8),sin(-pi/8));
end
Es=sum((abs(signal_space)).^2)/M; %调制信号能量Es
for i=1:1:length(snr)
No=Es/PsPn(i); %能量比等于功率比
for j=1:1:sample
x(j)=signal_space(floor(rand*M)+1); %输入随机信号
end
y=x+sqrt(No./2)*complex(randn(1,sample),randn(1,sample)); %信号叠加噪声
for j=1:1:sample
pysm=0;
for b=1:1:M
pysm=pysm+exp(-(abs(y(j)-signal_space(b))).^2/No);
end
log2py(j)=log2(1/M/(pi*No)*pysm); %信道容量公式第一项的和
end
Capacity(i)=-mean(log2py)-log2(pi*exp(1)*No); %mean函数可以直接得出矩阵的均值
end
EbNo=10*log10(PsPn./Capacity);
if k==2
plot(EbNo,Capacity,'b--'); %QPSK调制信道容量曲线
end
hold on;
if k==3
plot(EbNo,Capacity,'r-.'); %8PSK调制信道容量曲线
end
hold on;
if k==4
plot(EbNo,Capacity,'k:'); %16PSK调制信道容量曲线
end
hold on;
end
%16QAM调制的信道容量曲线
k=4;
M=2.^k;
signal_space(1) = complex(1,1);
signal_space(2) = complex(1,3);
signal_space(3) = complex(3,1);
signal_space(4) = complex(3,3);
signal_space(5) = complex(-1,1);
signal_space(6) = complex(-1,3);
signal_space(7) = complex(-3,1);
signal_space(8) = complex(-3,3);
signal_space(9) = complex(-1,-1);
signal_space(10) = complex(-1,-3);
signal_space(11) = complex(-3,-1);
signal_space(12) = complex(-3,-3);
signal_space(13) = complex(1,-1);
signal_space(14) = complex(1,-3);
signal_space(15) = complex(3,-1);
signal_space(16) = complex(3,-3);
Es=sum((abs(signal_space)).^2)/M; %调制信号能量Es
for i=1:1:length(snr)
No=Es/PsPn(i); %能量比等于功率比
for j=1:1:sample
x(j)=signal_space(floor(rand*M)+1); %输入随机信号
end
y=x+sqrt(No./2)*complex(randn(1,sample),randn(1,sample)); %信号叠加噪声
for j=1:1:sample
pysm=0;
for b=1:1:M
pysm=pysm+exp(-(abs(y(j)-signal_space(b))).^2/No);
end
log2py(j)=log2(1/M/(pi*No)*pysm); %信道容量公式第一项的和
end
Capacity(i)=-mean(log2py)-log2(pi*exp(1)*No); %mean函数可以直接得出矩阵的值
end
EbNo=10*log10(PsPn./Capacity);
plot(EbNo,Capacity,'c-');
xlabel('Eb/No(dB)');
ylabel('Capacity(比特每信道容量)');
axis([-2,20,0,5]);
grid on;
text(18,2.2,'QPSK') ;
text(18,3.2,'8PSK') ;
text(18,3.8,'16PSK' ) ;
text(18,4.2,'16QAM') ;
