前言
记录一下我做通信原理软件实验的过程.
常来我的博客看作业的同学可能会发现,我在别的课程都是直接粘贴上我的实验内容和结果.从来不写做实验的过程.
因为我觉得没有必要,但是这次的通信软件实验特殊,原因就是:我通信原理学的不好,说实话在实验开始后,到我码下这些字为止,我都没有具体的思路框架.好吧话不多说,开始吧.
准备:
四相相移键控信号简称“QPSK”。它分为绝对相移和相对相移两种。由于绝对相移方式存在相位模糊问题，所以在实际中主要采用相对移相方式QDPSK.
QPSK四相绝对相移键控，用载波信号的4个初始相位对应4进制码元。因此，对于输入的二进制数字序列每2比特分为一组，称为双比特码元AB，然后用4种不同的载波相位分别表征这4种数字码元。(比如说我们产生的二进制序列,00 01 10 11 一共四个双比特码元,用四种不同的载波相位分别对应这四个码元.)
四相相移调制是利用载波的四种不同相位差来表征输入的数字信息，是四进制移相键控QPSK在M=4时的调相技术，它规定了四种载波相位，分别为45°，135°，225°，315°，调制器输入的数据是二进制数字序列，为了能和四进制的载波相位配合起来，则需要把二进制数据变换为四进制数据，这就是说需要把二进制数字序列中每两个比特分成一组，共有四种组合，即00，01，10，11，其中每一组称为双比特码元。每一个双比特码元是由两位二进制信息比特组成，它们分别代表四进制四个符号中的一个符号。QPSK中每次调制可传输2个信息比特，这些信息比特是通过载波的四种相位来传递的。解调器根据星座图及接收到的载波信号的相位来判断发送端发送的信息比特。
QPSK相对于PSK技术的通信速率提升两倍，且可以有效地避免频率干扰
数字调制用“星座图”来描述，星座图中定义了一种调制技术的两个基本参数：（1）信号分布；（2）与调制数字比特之间的映射关系。星座图中规定了星座点与传输比特间的对应关系，这种关系称为“映射”，一种调制技术的特性可由信号分布和映射完全定义，即可由星座图来完全定义。
参考链接:
郭文博.基于Matlab的QPSK通信系统设计与仿真[J].信息通信,2019(11):211-212.
蒙特卡洛模拟(Monte Carlo Simulation)浅析
matlab生成随机数的rand、randi和randn三种形式
Matlab Code for Performance Analysis (BER vs Eb/N0) of BPSK, QAM, M-PSK, M-QAM, D-PSK, D-QAM etc..
QPSK Quadrature Phase Shift Keying (Basics, Modulator, Waveforms, Demodulator & Applications)
QPSK通信系统的Monte_Carlo仿真(这个判决的时候写的很好,很简练,我做的时候就没想不出来这么美的语句)
实验目的:
掌握QPSK通信系统的Monte Carlo仿真方法
掌握QPSK通信系统的组成原理
实验内容
未加信道纠错编码的QPSK调制通信系统
1)最大投影点准则进行判决
a.画出噪声方差\(\sigma^{2}\)分别为0、0.1、0.5、1.0时，在检测器输入端1000个接收到的信号加噪声的样本（星座图）；
b, 在AWGN信道下，分别画出数据点为1000、10000、100000时QPSK的Monte Carlo 仿真误比特率曲线和理论误比特率曲线，比较差别，分析数据点的数量对仿真结果的影响（横坐标是snr=Eb/No dB，格雷映射）
2）将检测器的判决准则改为最小欧式距离法（星座图上符号间的距离）
a,画出数据点为100000时QPSK的Monte Carlo仿真误比特率曲线和理论误比特率曲线（横坐标是snr=Eb/No dB，格雷映射，在AWGN信道下）；
b.比较与最大投影点准则进行判决结果的区别。
c.并理论分析证明两种判决方法是等价的。
实验框图及步骤
主函数及子函数编程思路
一步一步按照上面给的步骤来
产生二进制信号序列
例:X = rand(sz1,…,szN) 返回由随机数组成的 sz1×…×szN 数组，其中 sz1,…,szN 指示每个维度的大小。例如：rand(3,4) 返回一个 3×4 的矩阵,每个值位于0~1之间.
Binary_signal_sequence.m
%产生二进制信号序列,长度为N
function [a,b]=Binary_signal_sequence(length)
rand_num=randi([0 1],1,length);%产生跟序列长度个数一样的0到1的随机数
%ab表示一个四进制码元
%创建一个1×length的全零矩阵
a=zeros(1,length);
b=zeros(1,length);
for i=1:length
if(rand_num(i)<=0.25)
a(i)=0;b(i)=0;%随机数<0.25，规定码元值为为00
elseif(rand_num(i)<=0.5)
a(i)=0;b(i)=1;%随机数<0.5，规定码元值为为01
elseif(rand_num(i)<0.75)
a(i)=1;b(i)=0;%随机数<0.75，规定码元值为为10
else
a(i)=1;b(i)=1;%随机数<1，规定码元值为为11
end
end
end
第一次修改:
错误原因:这里我用的是使用 randi 函数（而不是 rand）生成在 0 和 1 之间均匀分布的随机整数。所以只能产生0或1.之所以当时用了randi的原因是在之前做Java的时候,要做随机,然后用的就是整数,之后再做通信软件实验的时候下意识给用成整数了.
修改后代码如下:
%产生二进制信号序列,长度为N
function [a,b]=Binary_signal_sequence(length)
rand_num=rand(1,length);%产生跟序列长度个数一样的0到1的随机数
%ab表示一个四进制码元
%创建一个1×length的全零矩阵
a=zeros(1,length);
b=zeros(1,length);
for i=1:length
if(rand_num(i)<=0.25)
a(i)=0;b(i)=0;%随机数<0.25，规定码元值为为00
elseif(rand_num(i)<=0.5)
a(i)=0;b(i)=1;%随机数<0.5，规定码元值为为01
elseif(rand_num(i)<0.75)
a(i)=1;b(i)=0;%随机数<0.75，规定码元值为为10
else
a(i)=1;b(i)=1;%随机数<1，规定码元值为为11
end
end
end
通过QPSK映射为四个信号点,采用gray码映射
参考:
（假设以二进制为0的值做为格雷码的0）
G：格雷码 B：二进制码 n：正在计算的位
根据格雷码的定义可得：
G(n) = B(n+1) XOR B(n)
即
G(n) = B(n+1) + B(n)
自低位至高位运算即可，无需考虑进位
gray_QPSK_mapping.m
function[gQm]=gray_QPSK_mapping(a,b,length)
gQm=zeros(2,length);%生成2×length矩阵
%通过格雷码映射
for i=1:length
if(a(i)==0&&b(i)==0)
gQm(1,i)=1;gQm(2,i)=0;
elseif(a(i)==0&&b(i)==1)
gQm(1,i)=0;gQm(2,i)=1;
elseif(a(i)==1&&b(i)==1)
gQm(1,i)=-1;gQm(2,i)=0;
elseif(a(i)==1&&b(i)==0)
gQm(1,i)=0;gQm(2,i)=-1;
end
end
end
产生两路正交的零均值高斯白噪声序列
这个没查资料,直接用的老师给的,在上面改了改
gaussian_sigma.m
%产生两路正交的零均值高斯白噪声序列，sigma为标准方差
function [noise]=gaussian_sigma(length,sigma)
noise=zeros(2,length); %写成矩阵模式便于和之前求得的gQm相加
for k=1:length
u=rand;
z=sigma*sqrt(2*log(1/(1-u)));
u=rand;
nc(k)=z*cos(2*pi*u);
ns(k)=z*sin(2*pi*u);
end
noise(1,:)=nc; %(1,:)表示第一行的元素
noise(2,:)=ns;
end
主函数:星座图(在检测器输入端1000个接收到的信号加噪声的样本)
这个前面弄完了就很简单了.调用子函数把图打出来就可以了.
exam_1_a.m
length=input('length='); %输入信源长度
s=input('方差='); %输入噪声方差
sigma=sqrt(s); %计算标准差
[a,b]=Binary_signal_sequence(length); %信源信号
gQm=gray_QPSK_mapping(a,b,length); %采用格雷码qpsk映射
noise=gaussian_sigma(length,sigma);%产生加性高斯白噪声
r=gQm+noise; %信号加噪
figure;
for i=1:length
if (gQm(1,i)==1&&gQm(2,i)==0)
plot(r(1,i),r(2,i),'*','color','red'); hold on;
elseif (gQm(1,i)==0&&gQm(2,i)==1)
plot(r(1,i),r(2,i),'*','color','yellow'); hold on;
elseif (gQm(1,i)==-1&&gQm(2,i)==0)
plot(r(1,i),r(2,i),'*','color','blue'); hold on;
else (gQm(1,i)==0&&gQm(2,i)==-1)
plot(r(1,i),r(2,i),'*','color','green'); hold on;
end
end
axis([-4 4 -4 4]);
line([4,-4],[0,0],'linewidth',2,'color','red')
line([0,0],[4,-4],'linewidth',2,'color','red')
title('星座图');
hold off最大投影准则判决
由相关度量准则\(C\left(r,\vec{s_{m}}\right)\)来看,最大投影就是加噪后的信号r和可能的gQm点乘最大.
可能的gQm有10、 01、 -10、 0-1四种,点乘后记录最大时所对应的值
参考:
之后的判决就是gary反编码,分别放在mj_a和mj_b里,再看他和原来信号ab的差异来计算误比特率即可
max_projection.m
%最大投影准则，将信道输出的r和原信号gQm逐个做向量积，for循环逐个比较求出最大向量积后判决
function [mj_a,mj_b]=max_projection(r,length)
mj_a=zeros(1,length);
mj_b=zeros(1,length);
C=zeros(4,length);
C_max=zeros(1,length);
for i=1:length
C(1,i)=1*r(1,i)+0*r(2,i); %投影（做向量积）
C(2,i)=0*r(1,i)+1*r(2,i);
C(3,i)=(-1)*r(1,i)+0*r(2,i);
C(4,i)=0*r(1,i)+(-1)*r(2,i);
C_max(i)=max([C(1,i),C(2,i),C(3,i),C(4,i)]); %求最大投影
if(C_max(i)==C(1,i))
mj_a(i)=0;mj_b(i)=0; %(1 0)判决为00
elseif(C_max(i)==C(2,i))
mj_a(i)=0;mj_b(i)=1; %(0 1)判决为01
elseif(C_max(i)==C(3,i))
mj_a(i)=1;mj_b(i)=1; %(-1 0)判决为11
else
mj_a(i)=1;mj_b(i)=0; %(0 -1)判决为10
end
end
end最小欧氏距离准则
跟上面的最大投影一样,不过是换成了r向量与各点的距离取最小的时候
min_distance.m
%最小距离准则，分别求r向量与各点的距离，for循环比较，求出最小距离，然后判决
function [mj_a,mj_b]=min_distance(r,length)
mj_a=zeros(1,length);
mj_b=zeros(1,length);
D=zeros(4,length);
D_min=zeros(1,length);
for i=1:length
D(1,i)=(r(1,i)-1)^2+r(2,i)^2; %求欧式距离的平方
D(2,i)=r(1,i)^2+(r(2,i)-1)^2;
D(3,i)=(r(1,i)+1)^2+r(2,i)^2;
D(4,i)=r(1,i)^2+(r(2,i)+1)^2;
D_min(i)=min([D(1,i),D(2,i),D(3,i),D(4,i)]); %求最小欧式距离
if(D_min(i)==D(1,i))
mj_a(i)=0;mj_b(i)=0; %(1 0)判决为00
elseif(D_min(i)==D(2,i))
mj_a(i)=0;mj_b(i)=1; %(0 1)判决为01
elseif(D_min(i)==D(3,i))
mj_a(i)=1;mj_b(i)=1; %(-1 0)判决为11
else
mj_a(i)=1;mj_b(i)=0; %(0 -1)判决为10
end
end
end
计算仿真误比特率
把判决后的信号mj_a和mj_b和原二进制序列信号a b逐个比较,计算误比特率
function [pb]=bit_error(mj_a,mj_b,a,b,length)
numbit=0;%numbit为错误比特个数
for i=1:length
if(mj_a(i)~=a(i))
numbit=numbit+1;
end
if(mj_b(i)~=b(i))
numbit=numbit+1;
end
end
pb=numbit/(2*length); %计算误比特率
end主函数:仿真误比特率曲线和理论误比特率曲线
根据Q函数\(\frac{1}{2}erfc\left ( \frac{a}{\sqrt{2}} \right )\)来计算理论误比特率
画出仿真误码率和理论误码率的曲线
length=input('length='); %输入信源长度
Eb=0.5; %发送信号的比特能量为0.5
SNR=0;
for i=1:18
SNR(i)=i-7;
N0=Eb/(10^(SNR(i)/10)); %给定信噪比求噪声的单边功率谱密度
sigma(i)=sqrt(N0/2); %求方差
[a,b]=Binary_signal_sequence(length); %信源信号
gQm=gray_QPSK_mapping(a,b,length); %qpsk映射,gray
noise=gaussian_sigma(length,sigma(i));%产生加性高斯白噪声
r=noise+gQm; %信号加噪
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%[mj_a,mj_b]=max_projection(r,length); %最大投影准则判决,判决输出信号
[mj_a,mj_b]=min_distance(r,length); %最小欧式距离准则判决,判决输出信号
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pb(i)=bit_error(mj_a,mj_b,a,b,length); %求仿真误比特率
pe(i)=0.5*erfc(sqrt(10^(SNR(i)/10))); %Q函数计算理论误比特率
end
figure;
semilogy(SNR,pb,'r-*'); %画出仿真误比特率曲线，设为红色
hold on;
semilogy(SNR,pe) %画出理论误比特率曲线，设为蓝色
grid on
xlabel('Eb/N0,dB')
ylabel('误比特率')
hold off; 
暂无评论内容