CN105897626A - 一种mppm调制的软解调算法及误符号率性能估计方法 - Google Patents

Abstract

本发明公开了一种MPPM调制的软解调算法及误符号率性能估计方法，软解调算法包括：产生随机的0、1序列，进行(N,M)MPPM信号调制；产生指数韦伯信道衰落和高斯白噪声，得到受高斯白噪声和信道衰落干扰的MPPM信号块，抽样，找到其中M个最大值，解调，仿真足够次数得到误符号率。该算法复杂度低且误符号率性能优于硬解调算法。误符号率估计方法包括：建立大气信道的数学模型，计算信号与非信号时隙抽样值的累积分布函数，分别求出M和N‑M个非空时隙随机变量最小值和最大值的概率密度函数，求出其差值ΔY的概率密度函数，计算ΔY的概率密度函数的积分，所得结果即为(N,M)MPPM的误符号率。解决了无法通过蒙特卡洛仿真算法直接对(N,M)MPPM的误符号率进行理论分析的问题。

Description

一种MPPM调制的软解调算法及误符号率性能估计方法

技术领域

本发明涉及无线光通信系统中的(N,M)MPPM调制方式在指数韦伯信道下的软解调计算方法，属于无线光通信技术领域。

背景技术

无线光通信又称自由空间光通信(free-space optical,FSO)，其采用大气信道作为传输媒介，而大气湍流导致的光学折射率随机起伏会使光信号在传输过程中产生光强闪烁、光束漂移和光束扩展等效应。这些效应引起的信道衰落会导致传输的光束质量下降，误码率劣化，进而影响系统的稳定性和可靠性。

为了研究大气湍流对FSO通信的影响，国内外科研人员对大气湍流信道模型进行了深入研究，先后提出了对数正态分布、K分布、gamma-gamma分布等模型，结果表明：对数正态分布模型适合于弱湍流，K分布适合于强湍流，而gamma-gamma分布则适合于强、中、弱三种湍流情况。但上述分布都不能很好地模拟采用有限孔径接收机无线光通信系统的光信号强度波动，于是2012年首次从可靠性工程领域引入指数韦伯(exponentiated Weibull，EW)分布以描述孔径平均效应存在情况下的光强衰落。

基于上述信道模型，各国研究人员对无线光通信的调制解调方式也进行了深入探讨和研究。对于调制方式，开关键控(On-Off Keying,OOK)作为最常用的调制方案，实现简单，但其抗干扰性能较差。而新型的多脉冲位置调制(Multipulse Pulse Position Modulation,MPPM)方式，其抗干扰能力相比于OOK具有很大的优势，并且带宽利用率较高。

对于解调方式，通常分为硬解调和软解调。硬解调通过固定、最佳或动态判决门限等方法把含噪声和衰落的信号判定为相应的调制电平，技术简单、易实现，但性能往往不够理想。软解调对于接收到的信号，不急于做判决，通过计算相应比特的软信息，提高判决的准确性，降低系统误码率，其性能远远优于硬解调。

目前的问题在于：由于MPPM的复杂性，关于其软解调特性的报道较少。目前已知软解调系统算法大多是基于条件信道下的计算，即只考虑存在加性高斯白噪声，不考虑大气导致的衰落。该算法不能用于指导实际的FSO通信系统接收机设计。

因此基于EW分布的无线光通信MPPM调制的软解调算法及其误符号率性能的估计实属当前十分重要的研究方向。

发明内容

本发明的目的是提出一种基于指数韦伯(EW)信道的无线光通信(N,M)MPPM调制的软解调算法并对其误符号率进行性能估计。其中软解调算法的优点在于无需判决门限，且判决的性能优于一般解调算法。之后提出的一种误符号率性能估计方法的数学模型，对上述(N,M)MPPM的软解调算法误符号率进行理论分析，进而正确有效的计算出理论误符号率。

本发明的目的是通过下述技术方案来实现的。

一种MPPM调制的软解调算法，包括下述步骤：

1)由MATLAB中rand函数产生0、1随机序列作为仿真输入比特流，并根据映射表进行(N,M)MPPM信号调制，其中，N为一个MPPM信号块中的时隙个数，M为N个时隙中的脉冲个数；

2)在MATLAB中由逆函数法产生符合指数韦伯分布的衰落信道系数并由randn函数产生高斯白噪声随机数；

3)将步骤2)产生的信道衰落系数和高斯白噪声随机数加入到步骤1)获得的(N,M)MPPM调制信号中，得到受到高斯白噪声和信道衰落干扰的MPPM信号块；

4)对步骤3)得到的(N,M)MPPM信号块中的每个时隙进行抽样；

5)利用MATLAB中的MAX函数找到抽样值中的最大值，将抽样值中M个最大值对应的时隙号判为1，其余均判为0；

6)利用步骤1)中(N,M)MPPM信号调制时的映射表，将步骤5)经过判决的信号块恢复为0、1序列，并将该序列与步骤1)中产生的信号对照，判断是否与其相同，若相同，则为正确解调；若不相同，则错误符号个数加1；

7)重复步骤1)到6)，仿真足够多的次数，得到错误符号个数，计算仿真的误符号率。

相应地，本发明给出了基于上述(N,M)MPPM调制的软解调算法进行误符号率性能估计的方法，包括下述步骤：

1)建立点对点无线光通信大气信道的数学模型；

2)计算信号时隙与非信号时隙抽样值的概率密度函数；

3)将接收到的(N,M)MPPM调制信号的每个时隙抽样值作为互相独立的随机变量，求出各自的累积分布函数，在此，N为一个MPPM信号块中的时隙个数，M为N个时隙中的脉冲个数；；

4)求出M个非空时隙随机变量最小值Y_smin的概率密度函数，再求出N-M个空时隙随机变量中最大值Y_nmax的概率密度函数；

5)得到非空时隙随机变量的最小值与空时隙随机变量的最大值的差值ΔY，求出差值的概率密度函数；

6)计算ΔY的概率密度函数从负无穷到零处的积分，所得结果即为(N,M)MPPM的误符号率。

本发明具有以下优点：

1)本发明中首先通过蒙特卡洛仿真实现了一种EW信道下MPPM调制的软解调算法，该方法主要通过在接收端对(N,M)MPPM块的每个时隙进行抽样得到前M个较大值，完成对脉冲所在时隙的判断。该方法的算法复杂度低，误符号率系统性能上优于硬解调算法；

2)提出一种适于理论计算软解调MPPM无线光通信系统误符号率性能的方法，解决了无法通过上述蒙特卡洛仿真算法直接对(N,M)MPPM的误符号率进行理论分析的问题。

附图说明

图1是本发明总体的流程步骤；

图2是弱湍流情况下误符号率拟合值、理论值与仿真值的对比图；

图3是中湍流情况下误符号率拟合值、理论值与仿真值的对比图；

图4是强湍流情况下误符号率拟合值、理论值与仿真值的对比图；

图5是中湍流情况下采用动态门限硬解调和软解调误符号率的对比图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合附图和具体实施方式进一步详细说明。本实施例仅表示对本发明的原理性说明，不代表对本发明的任何限制。

如图1所示，本发明以(5,2)MPPM为例，利用MATLAB软件实现蒙特卡洛仿真。详细介绍指数韦伯(EW)信道下的无线光通信(N,M)MPPM调制的软解调算法。

步骤1利用MATLAB中rand函数产生随机的0,1序列，并进行(5,2)MPPM调制

1a)首先需要设定仿真次数T，T的取值需要足够大以确保仿真的一般性，在本实例中我们取T＝10⁸。

根据MPPM调制原理可知，由于所以(5,2)MPPM调制可以对3比特(bit)为一组的比特流进行调制，其二进制比特流与MPPM符号之间的映射关系如表1所示，其中(l₁,l₂)表示有脉冲时隙所在位置，此表格表示的映射关系仅为一种实例。

1b)利用MATLAB中rand函数产生3个0,1随机数，作为一组对应映射表格进行调制，得到有5个时隙的MPPM信号块[x₁,x₂,…,x_i]，其中i＝5。

表1(5,2)MPPM调制映射表

3bit	(5,2)MPPM	(l1,l2)
			000	11000	(1,2)
001	10100	(1,3)
			010	10010	(1,4)
011	10001	(1,5)
			100	01100	(2,3)
101	01010	(2,4)
			110	01001	(2,5)
111	00110	(3,4)

步骤2产生符合EW的信道衰落和高斯白噪声

2a)已知EW型衰落下光强波动I的概率密度函数和累积分布函数分别为：

$f_{EW}\left(I\right)=\frac{\mathrm{\α}\mathrm{\β}}{\mathrm{\η}}{\left(\frac{I}{\mathrm{\η}}\right)}^{\mathrm{\β}-1}\exp \[-{\left(\frac{I}{\mathrm{\η}}\right)}^{\mathrm{\β}}\]\×{\{1-\exp \[-{\left(\frac{I}{\mathrm{\η}}\right)}^{\mathrm{\β}}\]\}}^{\mathrm{\α}-1},\≥0$

和

$F_{EW}\left(I\right)={\{1-\exp \[-{\left(\frac{I}{\mathrm{\η}}\right)}^{\mathrm{\β}}\]\}}^{\mathrm{\α}}$

其中，形状参数β＞0，其取值与闪烁指数相关。尺度参数η＞0，其取值由β决定并与辐照度均值相关。定义α＞0为在给定观测空间内，准直传播的光束与非准直传播的光束被成功接收的平均量。

根据以上公式，利用逆函数法，可以产生符合EW分布的随机数EW_rand。

2b)关于加性高斯白噪声，已知平均信噪比为：

$snr={\mathrm{\μ}}_n^2/2{\mathrm{\σ}}_n^2$

式中，μ为系数，μ²＝(RNP_t/M)²，R表示光电转换效率，P_t表示发射端的平均发射功率。本实例中设定R，P_t的值均为1，N与M分别等于5、2。则加性高斯白噪声的方差可以表示为：

${\mathrm{\σ}}_n^2=25/8snr$

根据上述公式得出的噪声方差利用randn函数即可产生符号高斯白噪声分布的随机数AWGN。

步骤3将信号通过信道，得到受噪声和信道衰落干扰的MPPM信号块

3a)对于接收到的MPPM信号而言，信号被抽样之后的输出模型可用以下式子表示：

$y=\left\{\begin{array}{c}\mathrm{\μ}h+n\\ n\end{array}\right.$

式中，由上到下依次为MPPM信号块中有脉冲的时隙和无脉冲的时隙，n表示的是均值为0，方差为的加性高斯白噪声，N₀为其单边功率谱密度。h表示信道衰落。

3b)将步骤2所产生的符合EW分布的5个随机数依次与一个(5,2)MPPM信号块中5个时隙的脉冲信号相乘。之后再将产生的5个高斯白噪声随机数依次与脉冲信号相加，最终得到接收的一个MPPM信号块。即，

y_i＝μ·x_i·EW_rand+AWGN

其中μ为系数，EW_rand为步骤2产生的信道随机数，AWGN为步骤2产生的高斯白噪声。

步骤4在接收端对MPPM信号进行抽样

对接收到的信号进行抽样，得到5个时隙的抽样值，即Y＝[Y₁,Y₂,Y₃,Y₄,Y₅]。

步骤5在抽样值中利用MATLAB函数找到最大值，进行判决

在Y＝[Y₁,Y₂,Y₃,Y₄,Y₅]中利用MATLAB中的MAX函数找到前两个最大值(不分次序)Y_i，Y_j，将其对应的时隙号i、j的信号值判为1，其余均判为0。如果有相同的数值而产生两个以上的最大值，则直接判断为解调错误，错误符号个数加1，返回步骤1，进行下一次仿真。

步骤6根据判决得到的脉冲所在时隙号解调信号

由步骤5可知脉冲时隙所在位置(l₁,l₂)为(i,j)(或(j,i)，按序号从小到大排列)对照表1中的映射表，恢复信号为3bit一组的比特流。对照步骤1中产生的信号，判断传输是否相同，若相同则正确解调；若不同或映射表中没有解调出来的情况则判断为解调错误，错误符号数加1。

步骤7得到仿真的误符号率

重复步骤1到6足够的次数T。完成所有解调后，通过错误个数和解调MPPM信号块的总个数，计算出误符号率。

注意，当取不同的平均信噪比snr进行计算时，可以得到不同信噪比下的误符号率。

如图1所示，在上述步骤的基础上，指数韦伯(EW)信道下无线光通信(N,M)MPPM软解调的误符号率性能估计方法，通过下述具体步骤进行计算：

步骤1建立点对点无线光通信大气信道的数学模型

设x为发送信号序列，y为接收信号序列，则信道的数学模型为：

y＝hx+n

n表示均值为0、方差为的加性高斯白噪声。h为服从EW分布的信道衰落。

步骤2计算信号时隙与非信号时隙抽样值的概率密度函数

根据已知EW型衰落下光强波动的概率密度函数和累积分布函数及步骤1，可以计算出信号时隙与非信号时隙抽样值的概率密度函数分别为：

$f\left(y\right|y_0)=f_n\left(y\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_n}\exp (-\frac{y^2}{2{\mathrm{\σ}}_n^2})$

$\begin{array}{c}f\left(y|y_1\right)=\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{y}{\mathrm{\μ}}\right)*f_n\left(y\right)\\ =\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\μ\σ}}_n}{\∫}_0^{\mathrm{\∞}}\frac{\mathrm{\α}\mathrm{\β}}{\mathrm{\η}}{\left(\frac{h}{\mathrm{\μ}\mathrm{\η}}\right)}^{\mathrm{\β}-1}\exp \[-{\left(\frac{h}{\mathrm{\μ}\mathrm{\η}}\right)}^{\mathrm{\β}}-\frac{{\left(y-h\right)}^2}{2{\mathrm{\σ}}_n^2}\]\\ \×{\left\{1-\exp \[-{\left(\frac{h}{\mathrm{\μ}\mathrm{\η}}\right)}^{\mathrm{\β}}\]\right\}}^{\mathrm{\α}-1}dh\end{array}$

其中f(y|y₁)为信号时隙抽样值的概率密度函数，f(y|y₀)为非信号时隙抽样值的概率密度函数。

步骤3将接收端每个时隙抽样值看作互相独立的随机变量，求出各自的累积分布函数

3a)将接收到的(N,M)MPPM调制信号划分为M个非空时隙以及N-M个空时隙，这些时隙之间不仅独立而且同分布，设接收到的(N,M)MPPM信号在N个时隙内的抽样变量为：

Y＝[Y_s1,Y_s2,…,Y_sM,Y_n1,Y_n2,…,Y_n(N-M)]

这里为了方便表示，不考虑二进制比特流与MPPM符号之间的映射关系。其中，Y_si表示的是M个信号中第i个时隙的抽样变量，Y_ni表示的是N-M个非信号时隙中第i个的抽样变量。则：

Y_smin＝min{Y_s1,Y_s2,…,Y_sM}

Y_nmax＝max{Y_n1,Y_n2,…,Y_n(N-M)}；

3b)根据步骤2，可分别求出Y_si的概率密度函数，进而可计算出其对应的累计分布函数为：

$\begin{array}{c}{F}_{Y_{si}}\left(y\right)=F\left(y|y_1\right)\\ ={\∫}_{-\mathrm{\∞}}^{y}f\left(t|y_1\right)dt\\ ={\∫}_{-\mathrm{\∞}}^y{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_n}\exp \left(-\frac{{\left(t-h\right)}^2}{2{\mathrm{\σ}}_n^2}\right)dhdt\\ ={\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)Q\left(\frac{h-y}{{\mathrm{\σ}}_n}\right)dh\end{array}$

式中，F(y|y₁)表示信号时隙抽样值的累积分布函数，为Q函数。

3c)同理，Y_ni的累积分布函数为：

$\begin{array}{c}{F}_{Y_{ni}}\left(y\right)=F\left(y|y_0\right)\\ ={\∫}_{-\mathrm{\∞}}^{y}f\left(t|y_0\right)dt\\ =\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_n}{\∫}_{-\mathrm{\∞}}^y\exp \left(-\frac{t^2}{2{\mathrm{\σ}}_n^2}\right)dt\\ =1-Q\left(\frac{y}{{\mathrm{\σ}}_n}\right)\\ =Q\left(-\frac{y}{{\mathrm{\σ}}_n}\right)\end{array}$

式中，F(y|y₀)表示非信号时隙抽样值的累积分布函数。

步骤4求出M个非空时隙随机变量最小值(记为Y_smin)的概率密度函数，再求出N-M个空时隙随机变量中最大值(记为Y_nmax)的概率密度函数

4a)由步骤3知非空时隙与空时隙抽样值的累积概率密度函数，根据概率论中求解最大最小分布的相关知识以及Q函数的性质1-Q(x)＝Q(-x)，Y_smin和Y_nmax的累积分布函数和可以分别表示为：

$\begin{array}{c}{F}_{Y_{s\min }}\left(y\right)=-1{\[1-F_{Y_{si}}\left(y\right)\]}^M\\ =1-{\[1-{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)Q\left(\frac{h-y}{{\mathrm{\σ}}_n}\right)dh\]}^M\\ =1-{\[{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)Q\left(\frac{y-h}{{\mathrm{\σ}}_n}\right)dh\]}^M\end{array}$

$F_{Y_{n\max }}\left(y\right)={\[F_{Y_{ni}}\left(y\right)\]}^{N-M}={\[Q(-\frac{y}{{\mathrm{\σ}}_n})\]}^{N-M};$

4b)对Y_smin和Y_nmax的累积分布函数关于y求导，可得到它们的概率密度函数和分别为：

$\begin{array}{c}{f}_{Y_{s\min }}\left(y\right)={\[F_{Y_{s\min }}\left(y\right)\]}^{\′}\\ ={\left\{{\[{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)Q\left(\frac{y-h}{{\mathrm{\σ}}_n}\right)dh\]}^M\right\}}^{\′}\\ =\frac{M}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_n}{\[{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)Q\left(\frac{y-h}{{\mathrm{\σ}}_n}\right)dh\]}^{M-1}\×\\ \[{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)\exp \left(-\frac{{\left(y-h\right)}^2}{2{\mathrm{\σ}}_n^2}\right)dh\]\end{array}$

$\begin{array}{c}{f}_{Y_{n\max }}\left(y\right)={\[F_{Y_{n\max }}\left(y\right)\]}^{\′}\\ ={\left\{{\[Q\left(-\frac{y}{{\mathrm{\σ}}_n}\right)\]}^{N-M}\right\}}^{\′}\\ =\frac{N-M}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_n}\exp \left(-\frac{y^2}{2{\mathrm{\σ}}_n^2}\right){\[Q\left(-\frac{y}{{\mathrm{\σ}}_n}\right)\]}^{-N-M-1}\end{array}.$

步骤5求出非空时隙随机变量最小值与空时隙随机变量最大值的差值ΔY，即ΔY＝Y_smin-Y_nmax，求出ΔY的概率密度函数。

5a)将ΔY看成两个随机变量相加Y_smin+(-Y_nmax)，则ΔY的概率密度函数等于Y_smin与-Y_nmax概率密度函数的卷积，其中-Y_nmax的概率密度函数为：

$\begin{array}{c}{f}_{-Y_{n\max }}\left(y\right)=f_{Y_{n\max }}\left(-y\right)\\ =\frac{N-M}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_n}\exp \left(-\frac{y^2}{2{\mathrm{\σ}}_n^2}\right){\[Q\left(\frac{y}{{\mathrm{\σ}}_n}\right)\]}^{N-M-1}\end{array}$

5b)再将与进行卷积运算，可求得ΔY的概率密度函数f_ΔY(y)为：

$\begin{array}{c}{f}_{\mathrm{\Δ}Y}\left(y\right)=\frac{M(N-M)}{2{\mathrm{\π\σ}}_n^2}{\∫}_{-\mathrm{\∞}}^{\mathrm{\∞}}\exp \[-\frac{{\left(y-t\right)}^2}{2{\mathrm{\σ}}_n^2}\]{\[Q\left(\frac{y-t}{{\mathrm{\σ}}_n}\right)\]}^{N-M-1}\\ {\[{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)Q\left(\frac{t-h}{{\mathrm{\σ}}_n}\right)dh\]}^{M-1}\[{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)\exp \left(-\frac{{\left(t-h\right)}^2}{2{\mathrm{\σ}}_n^2}\right)dh\]dt\end{array}$

由于ΔY为非空时隙随机变量最小值与空时隙随机变量最大值的差值，可知当ΔY为正，则非空时隙随机变量最小值大于空时隙随机变量最大值，信号可以正确解调；反之，当ΔY为负表示非空时隙随机变量最小值小于空时隙随机变量最大值，信号解调将发生错误。

步骤6计算ΔY的概率密度函数从负无穷到零处的积分，所得结果即为(N,M)MPPM的误符号率。

根据步骤5，对f_ΔY(y)从负无穷到零对y进行积分，即

$\begin{array}{c}{P}_s={\∫}_{-\mathrm{\∞}}^0{f}_{\mathrm{\Δ}Y}\left(y\right)dy\\ =\frac{M\left(N-M\right)}{2{\mathrm{\π\σ}}_n^2}{\∫}_{-\mathrm{\∞}}^0{\∫}_{-\mathrm{\∞}}^{\mathrm{\∞}}\exp \[-\frac{{\left(y-t\right)}^2}{2{\mathrm{\σ}}_n^2}\]{\[Q\left(\frac{y-t}{{\mathrm{\σ}}_n}\right)\]}^{N-M-1}\\ {\[{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)Q\left(\frac{t-h}{{\mathrm{\σ}}_n}\right)dh\]}^{M-1}\[{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)\exp \left(-\frac{{\left(t-h\right)}^2}{2{\mathrm{\σ}}_n^2}\right)dh\]dtdy\end{array}$

式中，P_s为误符号率。

将M＝2，N＝5，代入上式即可获得(5,2)MPPM在软解调时的误符号率。

进一步，为了化简上述误符号率计算过程，可以通过拟合近似代替f_Δy(y)的表达式，缩小计算量，进而求得P_s。

由于多重积分和高幂次数以及Q函数的存在，f_Δy(y)的表达式非常复杂，导致求解P_s的解析闭合表达式非常困难。为了方便计算误符号率表达式，引入一种函数拟合的方法去近似代替f_ΔY(y)的表达式，从而避免了多次积分和高幂次计算，大大缩小了计算量。由于P_s的表达式中是对f_ΔY(y)从负无穷到零这一范围进行积分，因此拟合时只需考虑这一部分的数据，对其进行拟合，再对所获得拟合表达式进行积分，最终得到误符号率的近似值。

借助于MATLAB软件中的曲线拟合工具箱(curve fitting tools)完成拟合。这里，使用高斯函数作为拟合函数对其进行拟合，高斯函数表达式为：

$G\left(y\right)=a\exp \[-{\left(\frac{y-b}{c}\right)}^2\]$

其中a,b,c为需要通过拟合获得的参数。为了提高拟合的准确性与稳定性，通常使用多个高斯函数相加进行拟合，用叠加后的表达式去表示所需拟合的曲线，其表达式为：

$G_n\left(y\right)=\underset{i=1}{\overset{n}{\Σ}}{a}_i\exp \[-{\left(\frac{y-b_i}{c_i}\right)}^2\]$

式中，a_i,b_i,c_i为需要通过拟合获得的参数；i表示序号。

对两项叠加的高斯函数进行积分运算，即：

$\begin{array}{c}{P}_s={\∫}_{-\mathrm{\∞}}^0\underset{i=1}{\overset{2}{\Σ}}{a}_i\exp \[-{\left(\frac{y-b_i}{c_i}\right)}^2\]dy\\ =\underset{i=1}{\overset{2}{\Σ}}{\∫}_{-\mathrm{\∞}}^0{a}_i\exp \[-{\left(\frac{y-b_i}{c_i}\right)}^2\]dy\end{array}$

进行变量替换，令z＝(y-b_i)/c_i，得：

$\begin{array}{c}{P}_s=\underset{i=1}{\overset{2}{\Σ}}{a}_i{c}_i{\∫}_{-\mathrm{\∞}}^{-b_i/c_i}\exp \left(-z^2\right)dz\\ =\underset{i=1}{\overset{2}{\Σ}}{a}_i{c}_i{\∫}_{b_i/c_i}^{\mathrm{\∞}}\exp \left(-z^2\right)dz\\ =\underset{i=1}{\overset{2}{\Σ}}\frac{a_i{c}_i\sqrt{\mathrm{\π}}}{2}erfc\left(\frac{b_i}{c_i}\right)\end{array}$

本发明的正确性和优点可通过以下仿真对比进一步说明：

通过MATLAB进行解析计算和蒙卡仿真，首先，确定拟合值的准确性；其次，通过理论值和仿真值的对比证明软解调算法的正确性；最后，再与中湍条件下的硬解调误符号率特性进行对比，凸显软解调算法的优越性。

仿真结果

图2、3和4分别给出了弱湍、中湍和强湍条件下(5,2)MPPM采用软解调时，误符号率与平均信噪比之间的关系曲线。计算用到的指数韦伯信道参数均来自已有文献。图中包含理论值、拟合值和仿真值。可以看出，三种湍流情况下，拟合值与理论值吻合很好，而且具有较高的准确性。另外，理论值也与仿真值吻合的很好，证明了理论值推导的正确。

图5给出了中湍条件下，(5,2)MPPM采用硬解调和软解调时误符号率与平均信噪比的关系曲线。硬解调选取的是性能最好的动态判决门限与软解调对比。对比发现，25mm孔径(α＝2.61,β＝1.37,η＝0.68)时，若要达到10^-7数量级的误符号率，动态判决门限和软解调所需的平均信噪比分别约为49.5和46dB；对于60mm(α＝1.52β＝4.02,η＝0.99)的孔径而言，若要达到同等数量级的误符号率，只需31.5和28dB的平均信噪比即可。由此可知，两种不同孔径时软解调比动态判决门限的硬解调有大约3.5dB的信噪比性能增益。因此，软解调方法在系统性能上优于硬解调。

以上所揭露的仅是本发明的较佳实施例而已，然而不能以此来限定本发明的权利范围，本领域技术人员利用上述揭露的技术内容做出些许简单修改、等同变化或修饰，仍属于本发明的保护范围之内。

Claims (10)

1.一种MPPM调制的软解调算法，其特征在于，包括下述步骤：

1)由MATLAB中rand函数产生0、1随机序列作为仿真输入比特流，并根据映射表进行(N,M)MPPM信号调制，其中，N为一个MPPM信号块中的时隙个数，M为N个时隙中的脉冲个数；

2)在MATLAB中由逆函数法产生符合指数韦伯分布的衰落信道系数并由randn函数产生高斯白噪声随机数；

3)将步骤2)产生的信道衰落系数和高斯白噪声随机数加入到步骤1)获得的(N,M)MPPM调制信号中，得到受到高斯白噪声和信道衰落干扰的MPPM信号块；

4)对步骤3)得到的(N,M)MPPM信号块中的每个时隙进行抽样；

5)利用MATLAB中的MAX函数找到抽样值中的最大值，将抽样值中M个最大值对应的时隙号判为1，其余均判为0；

6)利用步骤1)中(N,M)MPPM信号调制时的映射表，将步骤5)经过判决的信号块恢复为0、1序列，并将该序列与步骤1)中产生的信号对照，判断是否与其相同，若相同，则为正确解调；若不相同，则错误符号个数加1；

7)重复步骤1)到6)，仿真足够多的次数，得到错误符号个数，计算仿真的误符号率。

2.根据权利要求1所述的方法，其特征在于，所述步骤1)中，进行(N,M)MPPM信号调制，是通过下述方法来实现的：

1a)已知(N,M)MPPM调制的传信能力为：

$C_N^M=\frac{N!}{M!(N-M)!}$

式中，N为一个MPPM信号块中的时隙个数，M为N个时隙中的脉冲个数；设一个MPPM信号块以b个比特信号进行调制，则其传信能力需满足下述不等式：

$C_N^M\≥2^b$

这里，b取满足不等式的最大自然数；

1b)将产生的二进制0,1随机序列以b个比特为一组，依据调制映射表进行(N,M)MPPM调制。

3.根据权利要求1所述的方法，其特征在于，所述步骤2)中，产生符合指数韦伯分布的信道衰落系数和高斯白噪声随机数，是通过下述方法来实现的：

2a)已知指数韦伯型衰落下光强波动I的概率密度函数f_EW(I)和累积分布函数F_EW(I)分别为：

$f_{EW}\left(I\right)=\frac{\mathrm{\α}\mathrm{\β}}{\mathrm{\η}}{\left(\frac{I}{\mathrm{\η}}\right)}^{\mathrm{\β}-1}\exp \[-{\left(\frac{I}{\mathrm{\η}}\right)}^{\mathrm{\β}}\]\×{\{1-\exp \[-{\left(\frac{I}{\mathrm{\η}}\right)}^{\mathrm{\β}}\]\}}^{\mathrm{\α}-1},I\≥0$

和

$F_{EW}\left(I\right)={\{1-\exp \[-{\left(\frac{I}{\mathrm{\η}}\right)}^{\mathrm{\β}}\]\}}^{\mathrm{\α}}$

式中,形状参数β＞0，其取值与闪烁指数相关；尺度参数η＞0，其取值由β决定并与辐照度均值相关；定义α＞0为给定观测空间内，准直传播的光束与非准直传播的光束被成功接收的平均量；

利用逆函数法，产生符合指数韦伯分布的随机数；

2b)关于加性高斯白噪声，已知平均信噪比为：

$snr={\mathrm{\μ}}^2/2{\mathrm{\σ}}_n^2$

式中，μ为系数，μ²＝(RNP_t/M)²，R表示光电转换效率，P_t表示发射端的平均发射功率；

则加性高斯白噪声的方差可以表示为：

${\mathrm{\σ}}_n^2={\mathrm{\μ}}^2/2snr$

根据得出的噪声方差，利用MATLAB中的randn函数即可产生符合高斯白噪声分布的随机数。

4.根据权利要求1所述的方法，其特征在于，所述步骤3)中，得到受到高斯白噪声和信道衰落干扰的MPPM信号块，是通过下述方法来实现的：

3a)根据信号被抽样之后的输出模型：

$y=\left\{\begin{array}{c}\mathrm{\μ}h+n\\ n\end{array}\right.$

式中，由上到下依次为MPPM信号块中有脉冲的时隙和没有脉冲的时隙，n表示的是均值为0，方差为的加性高斯白噪声，N₀为其单边功率谱密度，h表示的是信道衰落；

3b)将步骤2)所产生的符合指数韦伯分布的N个随机数依次与一个(N,M)MPPM信号块中N个时隙的信号相乘，之后再将产生的N个高斯白噪声随机数依次与信号相加，最终得到接收的一个MPPM信号块。

5.一种权利要求1所述的MPPM调制的软解调算法进行误符号率性能估计方法，其特征在于，包括下述步骤：

1)建立点对点无线光通信大气信道的数学模型；

2)计算信号时隙与非信号时隙抽样值的概率密度函数；

3)将接收到的(N,M)MPPM调制信号的每个时隙抽样值作为互相独立的随机变量，求出各自的累积分布函数，在此，N为一个MPPM信号块中的时隙个数，M为N个时隙中的脉冲个数；

4)求出M个非空时隙随机变量最小值Y_smin的概率密度函数，再求出N-M个空时隙随机变量中最大值Y_nmax的概率密度函数；

5)得到非空时隙随机变量的最小值与空时隙随机变量的最大值的差值ΔY，求出差值的概率密度函数；

6)计算ΔY的概率密度函数从负无穷到零处的积分，所得结果即为(N,M)MPPM的误符号率。

6.根据权利要求5所述的方法，其特征在于，所述步骤2)中，计算信号时隙与非信号时隙抽样值的概率密度函数，是通过下述方法来实现的：

根据已知指数韦伯型衰落下光强波动的概率密度函数f_EW(I)和累积分布函数F_EW(I)及步骤1)中的数学模型，计算出信号时隙与非信号时隙抽样值的概率密度函数分别为：

$f\left(y\right|y_0)=f_n\left(y\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_n}\exp (-\frac{y^2}{2{\mathrm{\σ}}_n^2})$

$\begin{array}{c}f\left(y|y_1\right)=\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{y}{\mathrm{\μ}}\right)*f_n\left(y\right)\\ =\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\μ\σ}}_n}{\∫}_0^{\mathrm{\∞}}\frac{\mathrm{\α}\mathrm{\β}}{\mathrm{\η}}{\left(\frac{h}{\mathrm{\μ}\mathrm{\η}}\right)}^{\mathrm{\β}-1}\exp \[-{\left(\frac{h}{\mathrm{\μ}\mathrm{\η}}\right)}^{\mathrm{\β}}-\frac{{\left(y-h\right)}^2}{2{\mathrm{\σ}}_n^2}\]\\ \×{\left\{1-\exp \[-{\left(\frac{h}{\mathrm{\μ}\mathrm{\η}}\right)}^{\mathrm{\β}}\]\right\}}^{\mathrm{\α}-1}dh\end{array}$

式中，β＞0为形状参数，η＞0为尺度参数，定义α＞0为在给定观测空间内，准直传播的光束与非准直传播的光束被成功接收的平均量，h表示信道衰落，μ为系数，且μ＝RNP_t/M；f(y|y₁)为有信号时隙抽样值的概率密度函数，f(y|y₀)为无信号时隙抽样值的概率密度函数。

7.根据权利要求5所述的方法，其特征在于，所述步骤3)中，求出各自的累积分布函数，是通过下述方法来实现的：

3a)将接收到的(N,M)MPPM调制信号划分为M个非空时隙以及N-M个空时隙，这些时隙之间不仅独立而且同分布，设接收到的(N,M)MPPM信号在N个时隙内的抽样变量为：

Y＝[Y_s1,Y_s2…Y_sM,Y_n1,Y_n2…Y_n(N-M)]

其中，Y_si表示的是M个信号中第i个时隙的抽样变量，Y_ni表示的是N-M个非信号时隙中第i个的抽样变量；

则M个抽样值中的最小值Y_smin和N-M个抽样值中的最大值Y_nmax分别为：

Y_smin＝min{Y_s1,Y_s2,…,Y_sM}

Y_nmax＝max{Y_n1,Y_n2,…,Y_n(N-M)}；

3b)根据抽样值的概率密度函数，分别求出Y_si的概率密度函数，则可计算出其对应的累计分布函数为：

$\begin{array}{c}{F}_{Y_{si}}\left(y\right)=F\left(y|y_1\right)\\ ={\∫}_{-\mathrm{\∞}}^{y}f\left(t|y_1\right)dt\\ ={\∫}_{-\mathrm{\∞}}^y{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_n}\exp \left(-\frac{{\left(t-h\right)}^2}{2{\mathrm{\σ}}_n^2}\right)dhdt\\ ={\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)Q\left(\frac{h-y}{{\mathrm{\σ}}_n}\right)dh\end{array}$

其中，为加性高斯白噪声的方差，σ_n为加性高斯白噪声的标准差，为Q函数；

3c)同理Y_ni的累积分布函数为：

$\begin{array}{c}{F}_{Y_{ni}}\left(y\right)=F\left(y|y_0\right)\\ ={\∫}_{-\mathrm{\∞}}^{y}f\left(t|y_0\right)dt\\ =\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_n}{\∫}_{-\mathrm{\∞}}^y\exp \left(-\frac{t^2}{2{\mathrm{\σ}}_n^2}\right)dt\\ =1-Q\left(\frac{y}{{\mathrm{\σ}}_n}\right)\\ =Q\left(-\frac{y}{{\mathrm{\σ}}_n}\right)\end{array}.$

8.根据权利要求5所述的方法，其特征在于，所述步骤4)中，求解M个非空时隙随机变量最小值Y_smin的概率密度函数和N-M个空时隙随机变量中最大值Y_nmax的概率密度函数，通过下述方法实现：

4a)由得到的非空时隙与空时隙抽样值的累积概率密度函数，利用最大最小分布以及Q函数的性质1-Q(x)＝Q(-x)，Y_smin和Y_nmax的累积分布函数和可以分别表示为：

$\begin{array}{c}{F}_{Y_{s\min }}\left(y\right)=1-{\[1-F_{Y_{si}}\left(y\right)\]}^M\\ =1-{\[1-{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)Q\left(\frac{h-y}{{\mathrm{\σ}}_n}\right)dh\]}^M\\ =1-{\[{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)Q\left(\frac{y-h}{{\mathrm{\σ}}_n}\right)dh\]}^M\end{array}$

$F_{Y_{nmax}}\left(y\right)={\[F_{Y_{ni}}\left(y\right)\]}^{N-M}={\[Q(-\frac{y}{{\mathrm{\σ}}_n})\]}^{N-M};$

4b)对Y_smin和Y_nmax的累积分布函数关于y求导，可得到其概率密度函数和分别为：

$\begin{array}{c}{f}_{Y_{s\min }}\left(y\right)={\[F_{Y_{s\min }}\left(y\right)\]}^{\′}\\ ={\left\{{\[{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)Q\left(\frac{y-h}{{\mathrm{\σ}}_n}\right)dh\]}^M\right\}}^{\′}\\ =\frac{M}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_n}{\[{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)Q\left(\frac{y-h}{{\mathrm{\σ}}_n}\right)dh\]}^{M-1}\×\\ \[{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)\exp \left(-\frac{{\left(y-h\right)}^2}{2{\mathrm{\σ}}_n^2}\right)dh\]\end{array}$

$\begin{array}{c}{f}_{Y_{n\max }}\left(y\right)={\[F_{Y_{n\max }}\left(y\right)\]}^{\′}\\ ={\left\{{\[Q\left(-\frac{y}{{\mathrm{\σ}}_n}\right)\]}^{N-M}\right\}}^{\′}\\ =\frac{N-M}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_n}\exp \left(-\frac{y^2}{2{\mathrm{\σ}}_n^2}\right){\[Q\left(-\frac{y}{{\mathrm{\σ}}_n}\right)\]}^{N-M-1}\end{array}.$

9.根据权利要求5所述的方法，其特征在于，所述步骤5)中，求出非空时隙随机变量的最小值与空时隙随机变量的最大值的差值ΔY的概率密度函数，通过下述方法实现：

5a)将ΔY看成两个随机变量相加Y_smin＋(-Y_nmax)，则ΔY的概率密度函数等于Y_smin与-Y_nmax概率密度函数的卷积，其中-Y_nmax的概率密度函数为：

$\begin{array}{c}{f}_{-Y_{n\max }}\left(y\right)=f_{Y_{n\max }}\left(-y\right)\\ =\frac{N-M}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_n}\exp \left(-\frac{y^2}{2{\mathrm{\σ}}_n^2}\right){\[Q\left(\frac{y}{{\mathrm{\σ}}_n}\right)\]}^{N-M-1}\end{array};$

5b)再将与进行卷积运算，可求得ΔY的概率密度函数为：

$\begin{array}{c}{f}_{\mathrm{\Δ}Y}\left(y\right)=\frac{M(N-M)}{2{\mathrm{\π\σ}}_n^2}{\∫}_{-\mathrm{\∞}}^{\mathrm{\∞}}\exp \[-\frac{{\left(y-t\right)}^2}{2{\mathrm{\σ}}_n^2}\]{\[Q\left(\frac{y-t}{{\mathrm{\σ}}_n}\right)\]}^{N-M-1}\\ {\[{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)Q\left(\frac{t-h}{{\mathrm{\σ}}_n}\right)dh\]}^{M-1}\[{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)\exp \left(-\frac{{\left(t-h\right)}^2}{2{\mathrm{\σ}}_n^2}\right)dh\]dt\end{array}.$

10.根据权利要求5所述的方法，其特征在于，所述步骤6)中，得到(N,M)MPPM的误符号率，通过下述方法实现：

根据ΔY的概率密度函数求解误符号率P_S，对f_ΔY(y)从负无穷到零对y进行积分，即：

$\begin{array}{c}{P}_s={\∫}_{-\mathrm{\∞}}^0{f}_{\mathrm{\Δ}Y}\left(y\right)dy\\ =\frac{M\left(N-M\right)}{2{\mathrm{\π\σ}}_n^2}{\∫}_{-\mathrm{\∞}}^0{\∫}_{-\mathrm{\∞}}^{\mathrm{\∞}}\exp \[-\frac{{\left(y-t\right)}^2}{2{\mathrm{\σ}}_n^2}\]{\[Q\left(\frac{y-t}{{\mathrm{\σ}}_n}\right)\]}^{N-M-1}\\ {\[{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)Q\left(\frac{t-h}{{\mathrm{\σ}}_n}\right)dh\]}^{M-1}\[{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\μ}}{f}_{EW}\left(\frac{h}{\mathrm{\μ}}\right)\exp \left(-\frac{{\left(t-h\right)}^2}{2{\mathrm{\σ}}_n^2}\right)dh\]dtdy\end{array}.$