文章和代码曾经归档至【Github仓库：communication-system-simulation】或者公众号【AIShareLab】回复 通信零碎仿真 也可获取。

一、8PSK背景

略

二、原理概述

2.1 PSK调制

发送端发送的是一连串离散而随机的二进制比特流，应用PSK载波相位调制的办法，这样发送端发送的音讯便蕴含在了相位中，此种调制办法能够非常无效地节约带宽。

$$u_{m}(t)=A g_{T}(t) \cos (2 \pi f_{c} t+\frac{2 \pi m}{M}), m=0,1, \ldots, M-1 $$

其中, $g_{T}(t)$ 是发送滤波器的脉冲形态, 传输信号的频谱个性由它决定。A则是信号的幅度。在 $\mathrm{psk}$ 调制中, 所有的 $\mathrm{psk}$ 信号对于所有的 $\mathrm{m}$ 都具备雷同的能量。
能量为:

$$ \varepsilon_{m}=\int_{-\infty}^{+\infty} u_{m}^{2}(t) d t=\varepsilon_{s} \varepsilon_{s}$$
代表每个传输符号的能量。

在本次试验中, 为了不便剖析, 咱们令 $\mathrm{A}=1, g_{T}(t)=\sqrt{\frac{2 \varepsilon_{s}}{T}}, 0 \leq t \leq T$ , 那么, 相应的 $\mathrm{gsk}$ 调制信号的波形为

$$ \begin{array}{l} u_{m}(t)=\sqrt{\frac{2 \varepsilon_{s}}{T}} \cos (2 \pi f_{c} t+\frac{2 \pi m}{M})=\sqrt{\frac{2 \varepsilon_{s}}{T}}(A_{m c} \cos 2 \pi f_{c} t-A_{m s} \cos 2 \pi f_{c} t), \ m=0,1, \ldots, M-1,0 \leq t \leq T \ \end{array} $$

咱们规定:

$$\begin{array}{l}A_{m c}=\cos \frac{2 \pi m}{M} \A_{m s}=\sin \frac{2 \pi m}{M}\end{array}, m=0,1, \ldots, M-1.$$

通过上述剖析, 咱们不难得出, 这样一个相位调制信号能够看作两个正交载波, 因

此, 数字相位调制信号能够在几何上可用二维向量的模式来示意, 即 $$ \vec{s} *{m}=(\sqrt{\varepsilon*{s}} \cos \frac{2 \pi m}{M}, \sqrt{\varepsilon_{s}} \sin \frac{2 \pi m}{M}) $$

正交基函数为:

$$\begin{array}{l}\psi_{1}(t)=g_{T}(t) \cos 2 \pi f_{c} t \\psi_{2}(t)=-g_{T}(t) \sin 2 \pi f_{c} t\end{array}.$$

2.2 信号传输

调制信号在 AWGN 信道传输的时候, 会有噪声混淆进来, 此时输入信号变为: $$ r(t)=u_{m}(t)+n_{c}(t) \cos (2 \pi f_{c} t)-n_{s}(t) \sin (2 \pi f_{c} t) $$

其中, $n_{c}(t)$ 和 $n_{s}(t)$ 别离是加性噪声的同相重量和正交重量, 之后, 咱们将输入信号和 给出的基函数作相干, 则两个相关器的输入为: $r=s_{m}+n=(\sqrt{\varepsilon_{s}} \cos \frac{2 \pi m}{M}+n_{c}, \sqrt{\varepsilon_{s}} \sin \frac{2 \pi m}{M}+n_{s})$ 须要留神的是 $n_{c}(t)$ 和 $n_{s}(t)$ 这两个正交噪声的重量是零均值, 互不相干的高斯随机过程。

2.3 解调形式

（1）最小欧式间隔准则裁决

最小欧式间隔准则裁决: 求出接管到的信号向量与 M 个传输向量的欧式间隔, 选取 对应的最小欧式间隔的向量, 该向量对应的符号即为裁决输入符号。此种办法须要把握 间隔度量的概念并纯熟使用, 上面给出对于间隔度量具体的

实践剖析:

在接管音讯尚不确定 (即还没有接管到矢量 $\vec{r}$ ) 的状况下, 要使得先验概率为最大, 最好的裁决办法就是抉择具备最高先验概率 $P(\vec{s}_{m})$ 的信号; 承受到矢量 $\vec{r}$ 后, 其办法与 前者相似, 前者是寻找先验概率的最大值, 此时是寻找后验概率的最大值, 换言之, 抉择使 $P(\vec{s}_{m} \mid \vec{r})$ 最大的 $\vec{s}_{m}$ , 这个裁决准则称为最大后验概率 (MAP) 准则。

依据贝叶斯公式, 后验概率可示意为: $P(\vec{s} *{m} \mid \vec{r})=\frac{f(\vec{r} \mid \vec{s}* {m}) P(\vec{s} *{m})}{f(\vec{r})} $ 当 M 个信号先验概率相等, 因为 $f(\vec{r})$ 和 $P(\vec{s}* {m})$ 均为确定的值 $(P(\vec{s} *{m})=\frac{1}{M})$ 。则寻找 $P(\vec{s}* {m} \mid \vec{r})$ 的最大值就等价于寻找 $f(\vec{r} \mid \vec{s}_{m})$ 的最大值。此时 MAP 准则简化为 ML 准则。

咱们无妨对接管到的矢量 $\vec{r}$ 进行简要的剖析, $\vec{r}=\vec{s}_{m}+\vec{n}$, $\vec{s}_{m}$ 是信号矢量, $\vec{n}$ 是 AWGN 信道中的噪声矢量, 噪声矢量的重量 $n_{k}$ 遵从散布 $N(0, \frac{N_{o}}{2})$ , 则 $r_{k}$ 遵从散布 $N(s_{m k}, \frac{N_{0}}{2})$

因而 $$\begin{aligned} f(\vec{r} \mid \vec{s} *{m}) & =\prod*{k=1}^{N} \frac{1}{\sqrt{\pi N_{0}}} \mathrm{e}^{-\frac{(r_{k}-s_{m k})^{2}}{N_{0}}} \ & =\frac{1}{(\pi N_{0})^{\frac{N}{2}}} e^{\frac{|\vec{r}-\vec{s} *{m}|^{2}}{N*{0}}}, m=1,2, \ldots, M \end{aligned} $$

右端取对数有:

$\ln f(\vec{r} \mid \vec{s} *{m})=-\frac{N}{2} \ln (\pi N*{0})-\frac{1}{N_{0}} \sum_{k=1}^{N}(r_{k}-s_{m k})^{2} $

上式若要获得最大值, 不言而喻 $\sum_{k=1}^{N}(r_{k}-s_{m k})^{2}$ 须要取最小值。这也是合乎咱们直观印象的, 信号空间里两个信号点的欧氏间隔越小, 阐明它们越靠近。

因而, 定义间隔度量 $D(\vec{r} \mid \vec{s}_{m})$

如下:
$$ D(\vec{r} \mid \vec{s} *{m})=\sum*{k=1}^{N}(r_{k}-s_{m k})^{2}, m=1,2, \ldots, M $$

(2) 最佳检测器

最佳检测器将收到的信号向量 r 投射到 M 个可能的传输信号向量 $\{s_{m}\}$ 之一下来, 并 选取对应与最大投影的向量。将上述定义的间隔度量开展:

$$D(\vec{r} \mid \vec{s} *{m})=|\vec{r}|^{2}-2 \vec{r} \cdot \vec{s}* {m}+|\vec{s} *{m}|^{2}, m=1,2, \ldots, M $$

其中, $|\vec{r}|^{2}$项对所有的裁决度量是等价的的, 咱们疏忽这一项, 则失去相干度量:

$$ C(\vec{r} \mid \vec{s}* {m})=2 \vec{r} \cdot \vec{s} *{m}-|\vec{s}* {m}|^{2}, m=1,2, \ldots, M $$

能够看出, 间隔度量越小, 则相干度量越大。

上述剖析也证实了老师要求证实的内容：即相干度量与间隔度量是齐全等价的。

2.4 谬误概率

实践谬误概率：

8PSK: $\operatorname{erfc}(\sqrt{3 \times 10^{\frac{S N R}{10}}} \times \sin (\frac{\pi}{8}))$ ;

QPSK: $\operatorname{erfc}(\sqrt{10^{\frac{S N R}{10}}}) \times(1-0.25 \times \operatorname{erfc} \sqrt{10^{\frac{S N R}{10}}})$

理论谬误概率:

误码率: 错误码元/传输总码元

误比特率: 谬误比特/传输总比特

三、零碎框图

8PSK:

图3.1 8PSK零碎框图

四、主函数设计

4.1 星座图绘制主函数

1.流程图

图4.1 星座图绘制主函数流程图

2.代码实现

clc,clear,close;% Symbol sequence lengthL=100000;% Generate the original bit sequencesourceSeq=randnum(L);[pI,pQ,sourceSeqCode]=Map(sourceSeq,L);Eb = 1/3;errbit = zeros(1,26);errnum = zeros(1,26);SNR=-5:20;LS = length(SNR);for i=1:LS    % Find the one-sided power spectral density of noise for a given signal-to-noise ratio    N0=Eb/(10^(SNR(i)/10));     % variance    var(i)=N0/2  ;          [rI,rQ]=noise(var(i),pI,pQ);    [result,I]=judgment(rI,rQ);    [errbit(i),errnum(i)]=Count(result,I,sourceSeq,sourceSeqCode);enddraw(sourceSeqCode,rI,rQ);

4.2 QPSK与8PSK误码率比照主函数

1.流程图

图4.2 qpsk与8psk误码率比照主函数流程图

2.代码实现

clc,clear,close;% Symbol sequence lengthL=100000;% Generate the original bit sequencesourceSeq=randnum(L);[pI,pQ,sourceSeqCode]=Map(sourceSeq,L);Eb = 1/3;errbit = zeros(1,26);errnum = zeros(1,26);SNR=-5:20;LS = length(SNR);for i=1:LS    % Find the one-sided power spectral density of noise for a given signal-to-noise ratio    N0=Eb/(10^(SNR(i)/10));     % variance    var(i)=N0/2  ;          [rI,rQ]=noise(var(i),pI,pQ);    [result,I]=judgment(rI,rQ);    [errbit(i),errnum(i)]=Count(result,I,sourceSeq,sourceSeqCode);end% draw(sourceSeqCode,rI,rQ);% Calculate the theoretical bit error rate using the erfc functionTheory8PSKSER = erfc(sqrt(3*10.^(SNR/10)) * sin(pi/8)); % QPSK TheoryQPSKSER = erfc(sqrt(10.^(SNR/10))).*(1-0.25*erfc(sqrt(10.^(SNR/10))));;% Load the SER-SNR curve of QPSKload('qpsk_errnum');figure(2);epskerr = errnum/L;semilogy(SNR,epskerr,'r-o');hold on;semilogy(SNR,qpskerr,'b-o');hold on;semilogy(SNR,Theory8PSKSER,'k-');hold on;semilogy(SNR,TheoryQPSKSER,'k-');hold on;ylabel('SER');xlabel('SNR/dB')legend('8PSK仿真曲线','QPSK仿真曲线','实践曲线','Location', 'northeast' );grid on;axis([-5,20,10e-7,1]);

五、子函数设计

5.1 随机比特序列的产生

代码实现：

function [SourceSeq]=randnum(L)% L is the length of the generated sequence code% Since a code is composed of 3 bits, it is generated here with 3*L.randnum=rand(3*L,1);% Initialize the original sequenceSourceSeq=zeros(3*L,1);% The randomly generated sequence is judged% if the random number is greater than 0.5, it is judged as 1, otherwise it is judged as 0.for i=1:3*L    if(randnum(i)>=0.5)        SourceSeq(i)=1;    else        SourceSeq(i)=0;    endend

5.2格雷编码序列

代码实现

% Define 8psk mapping functionfunction [pI,pQ,SourceCode] = Map(SourceSeq,L)% pI - in-phase component% pQ - quadrature component% SourceCode - The size of the binary number of each digit of the sequence% initializationpI = zeros(L,1);pQ = zeros(L,1);% In order to facilitate subsequent expressions, sqrt(2)/2 is represented hereroot =sqrt(2)/2;% Constructing the mapping matrix according to the Gray code of 8PSKMappingMat = [[1,0];[root,root];[-root,root];[0,1];[root,-root];[0,-1];[-1,0];[-root,-root]];SourceCode =zeros(L,1);% mapping processfor i=1:L    % Since a source symbol is composed of three bits    % the low and high bits are read in reverse order here, and expressed in decimal    SourceCode(i)=SourceSeq(3*i-2)*4+SourceSeq(3*i-1)*2+SourceSeq(3*i)+1;    % Find the corresponding code through the decimal representation and map it    pI(i) = MappingMat(SourceCode(i),1);    pQ(i) = MappingMat(SourceCode(i),2);endend% Octal Gray Code Conversionfunction[a1,a2]=Map_other(N)% a1 is a sequence of random bits in binary% a2 is a sequence of octal symbols% Random bit sequence of length 3La1=bit(3*N);% a2 is used to store the code element sequence of length La2=zeros(1,N);for i=1:3:3*N-2   a2((i+2)/3)=abs(a1(i)*7-abs(a1(i+1)*3-a1(i+2)));% Converts binary bit sequence Gray-encoded to octal sequenceend end

5.3 映射函数

代码实现

% 8PSK coordinate mappingfunction [y3]=coordinate(x1,bit)% x1 is the encoded octal sequence, bit is the originally generated binary random sequenceN=length(x1);Es=bit*bit'/N;% The first line of y3 is used to store the abscissa% the second line is used to store the ordinatey3=zeros(2,N);% Coordinate mappingfor i=1:N    y3(1,i)=sqrt(Es)*cos(pi/4*x1(i)+pi/8);    y3(2,i)=sqrt(Es)*sin(pi/4*x1(i)+pi/8);endend

5.4 噪声生成与叠加输入

代码实现

% Generate Gaussian random noise sub-function , var is the variancefunction [rI,rQ] = noise(var,pI,pQ)L = length(pI);nc=zeros(L,1);ns=zeros(L,1);for k=1:Lu=rand;z=sqrt(var*2*log(1/(1-u)));nc(k)=z*cos(2*pi*u);ns(k)=z*sin(2*pi*u);end% Output two mutually orthogonal Gaussian signalsrI = pI+nc;rQ = pQ+ns;end

5.5 裁决函数

代码实现

% Judgment Criterion: Minimum Euclidean Distance Criterionfunction [result,I]=judgment(rI,rQ)root=sqrt(2)/2;L = length(rI);I=zeros(L,1);mapl = [[1,0];[root,root];[-root,root];[0,1];[root,-root];[0,-1];[-1,0];[-root,-root]];for i=1:L    index=0;    minp=100;    for j=1:8        % Traverse each coordinate        % find the point with the smallest Euclidean distance from the point, and record it.        if((mapl(j,1)-rI(i))^2+(mapl(j,2)-rQ(i))^2<minp)            minp=(mapl(j,1)-rI(i))^2+(mapl(j,2)-rQ(i))^2;            index=j;        end    end    I(i)=index;end% According to the judgment point, make the corresponding assignmentresult = zeros(3*L,1);for i=1:L    switch(I(i))        case 1            [result(3*i-2),result(3*i-1),result(3*i)]=SetValue(0,0,0);        case 2            [result(3*i-2),result(3*i-1),result(3*i)]=SetValue(0,0,1);        case 3            [result(3*i-2),result(3*i-1),result(3*i)]=SetValue(0,1,0);        case 4            [result(3*i-2),result(3*i-1),result(3*i)]=SetValue(0,1,1);        case 5            [result(3*i-2),result(3*i-1),result(3*i)]=SetValue(1,0,0);        case 6            [result(3*i-2),result(3*i-1),result(3*i)]=SetValue(1,0,1);        case 7            [result(3*i-2),result(3*i-1),result(3*i)]=SetValue(1,1,0);        case 8            [result(3*i-2),result(3*i-1),result(3*i)]=SetValue(1,1,1);    endendend

SetValue将对应的值进行赋值即可，因为裁决函数中调用次数过多，因而形象封装成一个函数，便于应用。

% SetValue functionfunction [a,b,c]=SetValue(d,e,f)a=d;b=e;c=f;end

5.6星座图绘制函数

代码实现

function draw(I,rI,rQ)figure;root = sqrt(2)/2;%映射矩阵MappingMat = [[1,0];[root,root];[-root,root];[0,1];[root,-root];[0,-1];[-1,0];[-root,-root]];L=length(I)for i=1:L    % Draw points with different colors according to different values    switch(I(i))        case 1            plot(rI(i),rQ(i),'*','color','r');hold on;        case 2            plot(rI(i),rQ(i),'*','color','g');hold on;        case 3            plot(rI(i),rQ(i),'*','color','b');hold on;        case 4            plot(rI(i),rQ(i),'*','color','c');hold on;        case 5            plot(rI(i),rQ(i),'*','color','m');hold on;        case 6            plot(rI(i),rQ(i),'*','color','y');hold on;        case 7            plot(rI(i),rQ(i),'*','color','k');hold on;        case 8            plot(rI(i),rQ(i),'*','color','[0.5,0.5,0.5]');hold on;    endendx = -4:0.1:4;y = -2:0.1:2;% plot(x,zeros(length(x)),'k');hold on;% plot(zeros(length(y)),y,'k');hold on;axis equalaxis([-2,2,-2,2]);end

5.7误码率及误比特率计算函数

代码实现

function [errbit,errsymbol]=Count(result,I,sourceSeq,sourceSym)% result: the sequence after the decision% I: symbol after decision% sourceSeq: original bit sequence% sourceSym: original symbol sequenceerrbit=0;errsymbol=0;for i=1:length(sourceSeq)    if(sourceSeq(i)~=result(i))        errbit=errbit+1;    endendfor i=1:length(I)    if(I(i)~=sourceSym(i))        errsymbol=errsymbol+1;    endendend

六、性能剖析与试验后果

6.1 比拟8PSK与QPSK的Monte Carlo仿真误符号率曲线、实践误符号率曲线差异

在AWGN信道下，未加信道纠错码的8PSK调制通信零碎检测器的裁决准则选为最小间隔法（星座图上符号间的间隔），格雷码映射，比拟数据点为100000时8PSK与QPSK的Monte Carlo仿真误符号率曲线，实践误符号率曲线，比拟差异（横坐标是SNR=Eb/N0）。(一张图上出现4条曲线)

图6. 1 QPSK与8PSK性能比拟

通过Monte Carlo仿真误符号率曲线能够看出，整体仿真后果根本合乎实践计算曲线，并且在不同的信噪比下，对应的QPSK的误码率显著小于8PSK的误码率，8PSK的星座图如下：

图6. 2 8PSK星座图

6.2 实践剖析8PSK性能比QPSK

实践证实如下：

$$ \begin{array}{l} r_{o 1}=\sqrt{\frac{\varepsilon_{g}}{2}} \cos (\frac{2 \pi m}{M})+n=\sqrt{\varepsilon_{s}} \cos (\frac{2 \pi m}{M})+n_{c} \ r_{o 2}=\sqrt{\frac{\varepsilon_{g}}{2}} \sin (\frac{2 \pi m}{M})+n=\sqrt{\varepsilon_{s}} \sin (\frac{2 \pi m}{M})+n_{s} \end{array} $$

假如发送信号的相位为 $\theta=0$ ，那么

$$ r_{o 1}=\sqrt{\varepsilon_{s}}+n_{c} \quad r_{o 2}=n_{s} $$

因为 $n_{c}$ 和 $n_{s}$ 是联结高斯随机过程
$$ f_{r}(r_{o 1}, r_{o 2})=\frac{1}{2 \pi \sigma_{r}^{2}} \exp (-\frac{(r_{o 1}-\sqrt{\varepsilon_{s}})^{2}+r_{o 2}^{2}}{2 \sigma_{r}^{2}}) $$

检测的测度为相位 $\Theta_{r}=\tan ^{-1}(r_{o 2} / r_{o 2})$

$$f_{\Theta_{r}}(\theta_{r})=\frac{1}{2 \pi} e^{-2 \rho_{s} \sin ^{2} \theta_{r}} \int_{0}^{\infty} v e^{-(v-\sqrt{4 \rho_{s}} \cos \theta_{r})^{2} / 2} d v$$

其中 $\rho_{s}=\varepsilon_{s} / N_{0}$ . $\mathrm{v}$ 是接管矢量 $\mathrm{r}$ 的包络. 若 $\rho_{s}>>1$ 且 $|\Theta_{r}|<=\pi / 2$

$f_{\Theta_{r}}(\theta_{r}) \approx \sqrt{\frac{2 \rho_{s}}{\pi}} \cos \theta_{r} e^{-2 \rho_{s} \sin ^{2} \theta_{r}}$ 若发送相位 0 , 当噪声使接管矢量的相位落在区域 $-\pi / M<\Theta_{r}<\pi / M$ 之外时会产生判 决谬误，即

$$\begin{array}{c}P_{e}=1-\int_{-\pi / M}^{\pi / M} f_{\Theta_{r}}(\theta_{r}) d \Theta \approx 1-\int_{-\pi / M}^{\pi / M} \sqrt{\frac{2 \rho_{s}}{\pi}} \cos \theta_{r} e^{-2 \rho_{s} \sin ^{2} \theta_{r}} d \theta_{r} \\approx 2 Q(\sqrt{2 \rho_{s}} \sin \frac{\pi}{M})=2 Q(\sqrt{2 \log _{2} M \rho_{b}} \sin \frac{\pi}{M})\end{array}$$

对于 M=4 ， 码元谬误概率应为:

$$P_{4}=1-(1-P_{2})^{2}=2 Q(\sqrt{\frac{2 \varepsilon_{b}}{N_{0}}})[1-\frac{1}{2} Q(\sqrt{\frac{2 \varepsilon_{b}}{N_{0}}})]$$

对于 M=8 , 码元谬误概率为:

$$P_{8}=2 Q(\sqrt{2 \frac{(\log _{2} 8) E_{b}}{N_{0}}} \sin \frac{\pi}{8})$$

很显著, 随着 $\mathrm{M}$ 的一直增长, $\mathrm{Pe}$ 也在一直减少。合乎试验后果。

七、问题回顾与总结

1.对二进制序列格雷编码的问题

针对二进制序列格雷编码，次要有两种思路，别离是间接法和间接法，间接法是先依据8PSK的格雷码结构映射矩阵，依据该3bit数示意码元的十进制值寻找其在格雷码矩阵中的对应地位，并且进行映射，须要留神的局部是因为一个源码元是由三个bit组成的，因而理论读取中，以3位单位进行遍历，并且通过倒序的形式读取低高位。

2.对于MPSK误符号率的问题

最开始计算误符号率时, 对于 QPSK 的实践误码率, 我最开始采纳的是 MPSK 的对立公式:

$$\begin{array}{c}P_{s, M P S K}=2 Q(\sqrt{2 \frac{E_{s}}{N_{0}}} \sin \frac{\pi}{M}) \=2 Q(\sqrt{2 \frac{k E_{b}}{N_{0}}} \sin \frac{\pi}{M})=2 Q(\sqrt{2 \frac{(\log _{2} M) E_{b}}{N_{0}}} \sin \frac{\pi}{M})\end{array}$$

然而在后续的学习中, 才发现上述公式尽能够在 M>4 的状况下才能够应用, 而对于 \mathrm{M}=4 时的零碎误码率, 应该采纳公式:

$$\begin{array}{cc}P_{4} & =1-P_{e} =2 Q(\sqrt{\frac{2 \varepsilon_{s}}{N_{0}}})[1-\frac{1}{2} Q(\sqrt{\frac{2 \varepsilon_{b}}{N_{0}}})]\end{array}$$

3.对于星座图绘制的问题

在绘制星座图时，初步想法是对于每一个裁决分类的样本点采纳不同的色彩绘制，然而对于如何针对点进行色彩，线性的绘制，我的初步想法是建设一个色彩-线性的向量，而后对于每个点裁决的具体情况找到对应的样本色彩线型，采纳数组援用的模式进行属性的赋值，然而随后发现看似简化了绘制过程，理论却在援用时产生很大的工作量，还可能产生谬误绘制，因而综合比拟下我抉择用switch的办法进行状况判断，并对相应的裁决点进行绘制。