张华娣，楼华勋,2

（1.中国电子科技集团公司第三十六研究所，浙江 嘉兴 314033；2.通信信息控制和安全技术重点实验室，浙江 嘉兴 314033）

1 引言

QAM（quadrature amplitude modulation）主要是通过两路正交载波的多种幅度来携带符号信息，因正交幅度调制具有多元的特性，又可记为MQAM。在正交幅度调制技术中，2 个支路的幅度具有多种取值，合成信号的相位以及幅度具有多种组合，星座图映射时以信号星座点之间的最小距离最大化为原则[1]。QAM 是一种高性能数字调制技术，近几十年来一直在发展创新，并因频谱利用率高和抗干扰能力强等优点广泛应用于各种有线通信、无线通信场合[2]，因此对其各阶调制的自动识别具有很高的实用意义，但是通过阅读大量文献，发现专门针对QAM 信号特别是高阶QAM 信号调制方式识别的研究比较少。文献[3]提出了一种基于通带QAM 信号的调制方式识别方法，该方法在通带内利用Hilbert 变换提取瞬时包络，然后统计瞬时能量分布向量作为特征向量来实现调制方式分类，该方法需要用到遗传算法对划分向量进行优化，并且需要用参考样本进行训练，计算量较大。文献[4]将特征参数和信号降阶算法相结合，该方法对恢复信号的幅度和相位要求很高，但是在盲识别的情况下，预处理得到的星座图很难和原星座图完全匹配，出现相偏的情况是很普遍的。文献[5]使用各阶QAM 信号包络平方的方差的理论值作为特征参数进行MQAM 信号类内识别，该算法需要对信噪比进行精确估计，不利于工程实现。文献[6]利用QAM 信号矢量图中最小环带的方差，完成5 种QAM 信号调制方式的识别，没有对高阶128QAM 和256QAM 信号进行讨论。文献[7]提出了利用信号星座点平均幅度半径的统计，实现QAM 类内识别，但是只针对16QAM和64QAM 信号。文献[8]使用四阶矩和减法聚类算法对8QAM、16QAM、64QAM 和128QAM 信号进行识别，但是没有对32QAM 和256QAM 进行识别。

文献[9]用减法聚类的方法恢复星座图及确定星座图的点数，从而完成QAM 信号的识别，但该方法只适用于进制数低于64 的QAM 信号。本文受到了文献[9]减法聚类方法的启发，提出通过对32QAM 和128QAM，以及64QAM 和256QAM 信号设置不同的聚类半径，计算聚类点的密度值之和的差异来进行类型的判别。同时结合文献[10]对方形QAM 和十字形QAM 分类的方法以及利用零中心归一化瞬时幅度紧致性对16QAM 识别的方法，完成了16QAM、32QAM、64QAM、128QAM、256QAM 信号的自动识别。

2 QAM 信号模型

接收端幅相调制信号模型表示为

其中，r(t)为接收信号；g(t)为成形滤波器的冲击响应；T0为码元周期；fc为载波频率；θc为载波相位；ε为定时偏差；N为观测符号数目；αi为零均值的平稳复随机序列，即发送码元序列；和ϕi分别为αi的幅度和相位；ω(t)为均值为零且单边功率谱密度为N0的平稳加性高斯噪声。

3 QAM 信号分类

QAM 信号根据星座图的形状可以分为方形QAM 和十字形QAM。图1 为16QAM、32QAM、64QAM、128QAM、256QAM 的星座图。

从图1 可以看出，16QAM、64QAM 和256QAM的星座图是方形的，属于方形QAM 信号；32QAM和128QAM 的星座图是十字形的，属于十字形QAM 信号。

4 QAM 信号识别方法

4.1 四阶累积量

信号的瞬时特征统计量反映的是信号的二阶统计特性，而信号的调制特点还反映在信号的高阶统计特性上，因此，在信号调制识别中还经常用到信号的高阶统计量作为特征参数[11]。又由于高斯噪声大于二阶的累积量恒为零，把接收到的含有高斯噪声的非高斯信号变换到累积量域处理，就可以剔除噪声的影响，而信号的各阶累积量又取决于信号的调制方式，因此可以利用接收到的带有高斯白噪声的信号的高阶累积量来建立识别参数，不但可以很好地抑制高斯噪声，而且可以识别其调制方式，具体要使用什么阶数的统计量依赖于具体问题。

假设s*(k)表示信号的复共轭，M pq表示信号的各阶矩，定义[12]

其中，E(⋅)表示取平均值，则有

图1 QAM 信号星座图

表1 列出了16QAM、32QAM、64QAM、128QAM、256QAM 信号对应的各阶累积量的理论值。

表1 信号的高阶累积量理论值

从表1 可以看出，16QAM 和32QAM 信号的四阶累积量|C40|理论值相差较大，而 64QAM 和256QAM 的四阶累积量|C40|的理论值与16QAM 的四阶累积量|C40|接近，并且128QAM 的四阶累积量|C40|的理论值与32QAM 的四阶累积量|C40|接近，因此理论上用这个参数直接就可以对方形QAM 信号（包括16QAM、64QAM 和256QAM 信号）和十字形QAM 信号（包括32QAM 和128QAM 信号）进行区分，而在实际应用中考虑到使用的特征参数的顽健性和通用性，也就是说，使用的参数对信号星座图伸缩具有不变性，并且与接收信号功率无关，考虑多使用一个参数|C42|，与|C40|进行比值处理，从而构造了特征参数F，其定义如式(7)所示。

用Matlab 软件在5 dB、10 dB、15 dB、20 dB、25 dB 和30 dB 信噪比下各随机产生100 组根升余弦成形（α=0.3）处理后的16QAM、32QAM、64QAM、128QAM、256QAM 信号，分别进行特征参数F统计估计，得到的F随信噪比变化曲线如图2 所示。其中，每组信号4 096 个码元，每个码元10 个采样点。从图2 可以看出，16QAM、64QAM和256QAM 信号的特征参数F值趋近于1，而32QAM 和128QAM 信号的特征参数F值小于1，用此特征参数很容易对方形QAM 信号和十字形QAM 信号进行分类。

图2 特征参数F 随信噪比的变化曲线

4.2 零中心归一化瞬时幅度紧致性

QAM 信号调制信息不仅体现在相位变化上，同时也体现在信号幅度变化上。由于QAM 信号相位变化种类较多，不宜用来进行QAM 信号类内识别。下面从方形QAM 信号的幅度分布入手，对其分布规律进行研究。

在理想条件下，16QAM信号幅度取值个数为3，64QAM信号幅度取值个数为9，256QAM信号幅度取值个数为32，对各阶方形QAM信号幅度进行最大值归一化，得到各方形QAM信号幅度及各幅度概率分布对应关系如表2～表4所示。从表2～表4可知，各阶方形QAM信号的幅度分布存在差异性。

表2 16QAM 信号幅度及各幅度概率分布对应关系

表3 64QAM 信号幅度及各幅度概率分布对应关系

文献[13]提到零中心归一化瞬时幅度紧致性特征参数反映了瞬时幅度分布的密集性，可以用该参数来反映各阶方形QAM 信号的幅度分布的差异性。零中心归一化瞬时幅度紧致性定义如式(8)所示。

其中，αcn(n)为零中心归一化瞬时幅度。

用Matlab 软件在5 dB、10 dB、15 dB、20 dB、25 dB 和30 dB 信噪比下各随机产生100 组根升余弦成形处理（α=0.3）后的各阶QAM 信号，分别进行特征参数统计估计，得到的随信噪比变化曲线，如图3 所示。其中，每组信号4 096 个码元，每个码元10 个采样点。

从图3 可以看出，64QAM 和256QAM 的零中心归一化瞬时幅度紧致性的统计值比较接近，不宜用该特征参数进行两者的区分，但是，16QAM 信号的零中心归一化瞬时幅度紧致性的统计值明显比64QAM 和256QAM 的小，因此，可以通过此参数将16QAM 从方形QAM 中分离出来。

表4 256QAM 信号幅度及各幅度概率分布对应关系

图3 特征参数随信噪比的变化曲线

4.3 聚类点密度值差异判别法

4.3.1 减法聚类算法

假设有N个数据{x1,x2,…,xN},将每一个数据点作为一个潜在的聚类中心，计算每个点的密度，密度定义为

其中，正常数rb定义了一个使密度指标显著减小的邻域，这样使靠近第一个聚类中心的数据点的密度指标明显下降，减小其作为下一个聚类中心的可能性，k2表示不同类型调制信号的调整权值。选取更新后具有最大密度值的点作为第二个聚类中心为其密度值。依次类推，可以得到第k阶聚类中心和密度值。

4.3.2 邻域半径设置

减法聚类算法中，有2 个常数（ra和rb）是?未知的，这2 个常数的大小如何设定是该算法的难题。在调制样本已知的情况下，可以通过计算机仿真方法选取，但在调制阶数值未知的情况下，无法找到一个合适的常数ra，满足高阶和低阶调制信号都可使用减法聚类算法。当以低阶调制为准选择ra值时，减法聚类对高阶调制分离性比较差；当以高阶调制为准选择ra值时，减法聚类对低阶调制的聚合性比较差。由此可见，选用单一的聚类半径并不能很好地解决2 种信号的分类识别。由此，针对不同的调制方式对应的理论最大聚类半径对ra和rb进行设置，假设MQAM 信号的水平方向的星座点对应的最大距离为L，由图1 可知，32QAM 信号相邻2个星座点的距离为，64QAM 信号相邻2 个星座点的距离为，128QAM 信号相邻2 个星座点的距离为，256QAM 信号相邻2 个星座点的距离为，由2 个星座点之间的距离来确定最大聚类半径，则32QAM、64QAM、128QAM 和256QAM 对应的ra和rb可以分别设置为

4.3.3 密度差异判别法

文献[14]提出，假设当前聚类星座图为待识别信号中最简单的星座结构，以相应的邻域半径进行减法聚类。如果聚类正确，聚类点数必定不大于星座点数；反之，说明当前假设不正确，继续假设当前聚类星座图为较复杂的星座图结构，进行下一轮减法聚类，依次类推，直到聚类点数小于或等于当前假设的星座图中的点数。

但是，对于不同信噪比情况下，用文献[14]提到的判断聚类点数小于当前假定的星座图中的点数的方法时，需要根据信噪比来调整ra和rb的大小，这显然增加了算法的难度，不利于工程实现。由此，本文使用普适性更强的密度差异判别法来进行信号类型的判断。假设要识别的信号属于十字形QAM 信号（32QAM 或128QAM），分别设置聚类半径为，用减法聚类算法进行聚类，分别得到2 种聚类中心点，根据式(9)每个聚类中心点对应一个密度值，在假设正确的情况下，得到的聚类点的密度值之和必然是最大的，由此，根据算出的2 个聚类点的密度值之和的差异就可以对32QAM或128QAM 信号进行判断。同理，对于64QAM 或256QAM 信号，分别设置聚类半径为，用同样的方法可以进行类型判别。

5 自动识别流程

识别步骤如下。

1)对采样信号进行频率粗估计，然后根据估计的频率进行正交下变频。

2)对正交下变频后的信号进行频率的精确估计，可用四次方谱和三点法得到精确估计的载频值，具体见文献[15]，用估计到的载频值来消除载波的影响。

3)估计信号带宽，并进行低通滤波，以减小噪声引起的相位抖动，具体见文献[9]。

4)计算特征参数F，设置门限参数0.635，计算出的特征参数F的值大于该门限的信号判断为方形QAM 信号，否则为十字形QAM 信号。

5)若判断为方形QAM 则执行步骤6)，否则直接执行步骤8)。

6)对信号的信噪比进行粗估，当信噪比小于8 dB 时，设置零中心归一化瞬时幅度紧致性特征参数门限为2.43；当信噪比为8～12 dB 时，设置门限为2.40；当信噪比大于12 dB 时，设置门限为2.38。

8)通过包络的平方来求取波特率，估计定时误差，在最佳采样点进行定时取样得到较为理想的星座图，具体插值滤波和解卷方法可见文献[16]。

6 仿真实验

实验1用Matlab 软件分别产生信噪比20 dB的16QAM、32QAM、64QAM、128QAM、256QAM信号各一组，根升余弦成形滤波器α=0.3，信号实际载频为1 MHz，采样率为10 MHz，码元速率为250 kHz。通过Matlab 软件对这5 种信号的自动识别流程进行仿真，仿真内容和结果如下。

图4 为5 种信号的调制波形图。通过对仿真信号求FFT，取幅度最大值的方法进行频率估计。

图5 为5 种信号根据粗估计的频率进行正交下变频后的波形图。通过三点法对正交下变频之后的信号做频谱精确估计。通过对正交下变频之后的信号做功率谱，以低于最大功率3 dB 为准则进行带宽估计。

图6 为5 种信号正交下变频后并进行低通滤波后的波形图。根据式(7)对正交下变频并低通滤波之后的信号进行特征参数F的计算。表5 为5 种信号的频率粗估计值、频率精确估计值、带宽估计值和特征参数F。

图4 5 种信号调制波形

图5 5 种信号正交下变频后波形

图6 5 种信号低通滤波后波形

表5 频率粗估计值、频率精确估计值、带宽估计值和特征参数F

结合表5 和自动识别流程步骤4)，16QAM、64QAM和256QAM的特征参数F的值大于该门限，判断为方形QAM 信号，32QAM 和128QAM 被判断为十字形QAM 信号。

根据式(8)对正交下变频并低通滤波之后的信号进行中心归一化瞬时幅度紧致性特征参数的计算，表6 为16QAM、64QAM 和256QAM 信号的零中心归一化瞬时幅度紧致性特征参数。

表6 零中心归一化瞬时幅度紧致性特征参数

结合表6 和步骤7)，16QAM 信号被正确区分出来，同时，64QAM 和256QAM 信号进入步骤8)继续识别。

通过求取正交下变频并低通滤波之后的信号的包络的平方的功率谱，取相应的最大值对应频率点的方法求取信号的波特率，表 7 为32QAM、64QAM、128QAM 以及256QAM 信号的波特率。

由表7 估计的波特率结合实际信号采样率，可以得到每个符号的实际采样点数，再结合基于过零检测的定时误差算法，可得到每个符号的定时误差估计值，通过定时误差估计值来控制内插滤波器对采样得到的信号样本值进行插值运算，从而得到信号在最佳采样时刻的近似值。

表7 波特率

图7 为 32QAM、64QAM、128QAM 以及256QAM 信号在最佳采样点进行定时取样得到的较为理想的星座图仿真结果。

根据4.3.2 节邻域半径计算方法，首先要得到MQAM 信号的水平方向的星座点对应的最大距离L，表8 为32QAM、64QAM、128QAM 以及256QAM信号的水平方向的星座点对应的最大距离L的估计值。

图7 4 种信号恢复星座图

表8 最大距离L

根据4.3.1 节所提的减法聚类算法，计算得到32QAM信号在邻域半径下的聚类中心及相应的密度值。

图8 32QAM 信号在邻域半径下得到的聚类中心和密度值

通过将图8 和图9 的聚类中心对应的密度值进行相加，可以得到32QAM 信号在2 种邻域半径下聚类点数以及密度值之和。表9 为32QAM 信号在2 种邻域半径下的聚类点数以及密度值之和。

结合表9 和自动识别流程步骤9)，32QAM 信号可以被正确区分。

根据4.3.1 节所提的减法聚类算法，计算得到64QAM 信号在邻域半径下的聚类中心及相应的密度值。

图9 32QAM 信号在邻域半径下得到的聚类中心和密度值

表9 32QAM 信号在2 种邻域半径下的聚类点数以及密度值之和

图10 6 4QAM 信号在邻域半径下得到的聚类中心和密度值

通过将图10 和图11 的聚类中心对应的密度值进行相加，可以得到64QAM 信号在2 种邻域半径下聚类点数以及密度值之和。表10 为64QAM 信号在2 种邻域半径下的聚类点数以及密度值之和。

结合表10 和自动识别流程步骤10)，64QAM信号可以被正确区分。

根据4.3.1 节所提的减法聚类算法，计算得到128QAM 信号在邻域半径下的聚类中心及相应的密度值。

图11 64QAM 信号在邻域半径下得到的聚类中心和密度值

表10 64QAM 信号在2 种邻域半径下的聚类点数以及密度值之和

图12 128QAM 信号在邻域半径下得到的聚类中心和密度值

通过将图12 和图13 的聚类中心对应的密度值进行相加，就得到了128QAM 信号在2 种邻域半径下聚类点数以及密度值之和。表11 为128QAM 信号在2 种邻域半径下的聚类点数以及密度值之和。

图13 128QAM 信号在邻域半径下得到的聚类中心和密度值

表11 128QAM 信号在2 种邻域半径下的聚类点数以及密度值之和

结合表11 和自动识别流程步骤9)，128QAM信号可以被正确区分。

根据4.3.1 节所提的减法聚类算法，计算得到256QAM 信号在邻域半径下的聚类中心及相应的密度值。

图14 256QAM 信号在邻域半径下得到的聚类中心和密度值

图15 256QAM 信号在邻域半径下得到的聚类中心和密度值

通过将图14 和图15 的聚类中心对应的密度值进行相加，可以得到256QAM 信号在2 种邻域半径下聚类点数以及密度值之和。表12 为256QAM 信号在2 种邻域半径下的聚类点数以及密度值之和。

表12 256QAM 信号在2 种邻域半径下的聚类点数以及密度值之和

结合表12 和自动识别流程步骤10)，256QAM信号可以被正确区分。

实验2用Matlab 软件在5 dB、10 dB、15 dB、20 dB、25 dB 和30 dB 信噪比下各随机产生100 组16QAM、32QAM、64QAM、128QAM、256QAM 信号，根升余弦成形滤波器α=0.3，每个码元采样点数10，码元个数为4 096，得到识别正确率如表13 所示。

表13 识别正确率

由实验2 可知，除了256QAM，其他信号在信噪比为15 dB 时信号类型识别正确率都能达到90%以上，说明本文所提的联合识别方法是可行的。

实验3用Matlab 软件在信噪比为15 dB 情况下，分别随机产生100组16QAM、32QAM、64QAM、128QAM、256QAM 信号，根升余弦成形滤波器α=0.3，每个码元对应采样点数10，在码元个数分别为256 个、1024 个、4096 个时，得到的识别正确率如表14 所示。

表14 识别正确率

由实验3 可知，识别正确率与码元个数有关。这是由于QAM 信号本身幅相种类比较多，只有达到一定的码元个数，才能反映各QAM 信号的幅相概率分布情况。当码元个数达到4 096 个时,能得到比较好的识别正确率。

7 结束语

本文针对MQAM 调制方式识别算法进行研究，提出了一种联合四阶累积量、零中心归一化幅度紧致性和减法聚类计算聚类中心密度值的自动识别方法。本文算法不需要预先知道信号的波特率、载波频率，而是通过非线性变换和频谱分析，从接收信号中估计出这些参数，然后逐级进行信号调制类型的识别，该方法运算量适中，易于工程实现。