一种基于循环谱相关的常用数字调制信号识别方法

摘要

本发明公开了一种基于循环谱相关的常用数字调制信号识别方法，利用信号循环谱所具有的抗噪特性来提高信号分析的可靠性；并在信号谱相关函数的计算过程中引入α截面小波去噪和叠加求取平均的环节，有效地减弱了原谱相关估计算法结果中由采样点数受限和外界干扰引起的随机波动，利于调制特征的识别和提取；同时，利用信号谱相关计算所获取谱相关图的α截面和f截面，选取合适的特征和参数(如谱相关函数α截面和f截面最大绝对值比、α截面强谱线数目、α截面波动系数、f截面归一化面积、α截面谱线的显著度比等)构建分类方法对通信信号的调制方式进行识别。

说明书

技术领域

本发明属于通信信号调制识别领域，具体涉及一种基于循环谱相关的常用数字调制信号识别方法。

技术背景

通信信号的调制方式识别是介于信号获取和解调之间的一个重要环节，无论是在军用还是民用通信领域中都扮演着重要角色，尤其是在动态无线电频谱的管理和不明干扰的鉴别工作中。另外，通信信号的调制方式识别也是构建软件无线电或认知无线电应用的基础，也为多体制通信互联应用提供技术支持。

基于傅里叶变换的信号功率谱分析是一种常用的较为经典的针对平稳信号的分析方法。通过信号的功率谱分析，可以较好地获取所分析信号的基本特征，包括载波频率、频带宽度等信号参数。但在通信过程中，很大一部分信号由于人工干预(如调制、编码、扫描等)而转换为非平稳信号，功率谱分析已不能较为全面地揭示信号的调制特征，因此基于功率谱分析的调制识别方法通常要借助较多时域统计特征，局限性较大。Gardner提出的循环谱相关理论是近年来人们研究非平稳信号的一个重要工具，该理论充分考虑了信号经过周期性人工干预后所具有的重要性质——循环平稳性，可以较为全面地揭示通信信号的调制特征，这是经典功率谱分析无法比拟的。

虽然循环谱相关理论在各领域的非平稳信号分析中得到较为广泛的应用，但目前在通信信号调制识别领域并未普及。

发明内容

本发明要解决的技术问题是：克服现有技术的不足，提供一种基于循环谱相关的常用数字调制信号识别方法，利用信号循环谱所具有的抗噪特性来提高信号分析的可靠性；并在信号谱相关函数的计算过程中引入α截面小波去噪和叠加求取平均的环节，有效地减弱了原谱相关估计算法结果中由采样点数受限和外界干扰引起的随机波动，利于调制特征的识别和提取；同时，利用信号谱相关计算所获取谱相关图的α截面和f截面，选取合适的特征和参数(如谱相关函数α截面和f截面最大绝对值比、α截面强谱线数目、α截面波动系数、f截面归一化面积、α截面谱线的显著度比等)构建分类方法对通信信号的调制方式进行识别。

本发明解决其技术问题所采用的技术方案为：一种基于循环谱相关的常用数字调制信号识别方法，包括以下几个步骤：

(1)将截获的信噪比调制信号经采样后分为n等份，作为进行谱相关运算的输入信号；

(2)选取合适的傅立叶变换点数和平滑窗宽度，用基于频域平滑的谱相关估计算法分别对这n份信号作谱相关运算；

所述基于频域平滑的谱相关估计算法如下：

$>S_{x\Delta t}^{\alpha }{(t,f)}_{\Delta f}=\frac{1}{M}\underset{v=-(M-1)/2}{\overset{(M-1)/2}{\Sigma }}\frac{1}{\Delta t}{X}_{\Delta t}(t,f+\frac{\alpha }{2}+{\mathrm{vF}}_S)X_{\Delta t}^{*}(t,f-\frac{\alpha }{2}+{\mathrm{vF}}_S)$>

$>X_{\Delta t}(t,f)=\underset{K=0}{\overset{N-1}{\Sigma }}{a}_{\Delta t}\left({\mathrm{KT}}_S\right)x_k(t-{\mathrm{KT}}_S)\exp (-j2\pi f(t-{\mathrm{KT}}_S\left)\right)$>

其中，为谱相关运算后的结果；X_Δt(t,f)为信号x_k(t)短时傅立叶变换的结果；信号x_k(t)为所截获信号的n等份之一；Δt为x_k(t)持续时间；a_Δt是窗函数；Δf为频域平滑间隔；F_S为频域最小增量单位；T_S是时域采样间隔；N是信号x_k(t)样本数；M为平滑窗宽度系数；α为循环频率；f为频率；t为时间；

(3)步骤(2)中的运算结果经α截面小波去噪后相加取平均值，得到改进的谱相关估计算法的输出结果；

所述改进的谱相关估计算法如下：

$>S_{x\Delta t}^{\alpha }\left(f\right)=\frac{1}{N}\underset{k=1}{\overset{N}{\Sigma }}{\mathrm{WDN}}_{\alpha }\left[S_{x\Delta t}^{\alpha }(t_k,f)\right]$>

其中，WDN_α(·)表示对中的循环频率α进行小波去噪；

(4)依据步骤(3)所得谱相关函数的α截面和f截面提取5种调制识别特征，并结合一种基于时域统计的信号特征构建分类方法对几种数字调制信号进行识别；所述5种调制识别特征为：谱相关函数α截面和f截面最大绝对值比(记为R₁)、α截面强谱线数目(记为R₂)、α截面波动系数(记为R₃)、f截面归一化面积(记为R₄)、α截面谱线的显著度比(记为R₅)；所述一种基于信号时域的统计特征为：零中心归一化瞬时幅度绝对值的标准偏差(记为R₆)；所述几种数字调制信号包括2ASK、4ASK、BPSK、QPSK、8PSK、MSK、2FSK、4FSK、2FSK*、4FSK*，其中FSK*表示码元初始相位不相关的频移键控；

所述对几种数字调制信号的识别过程为：

若则判定该信号的调制方式为MSK；

若$>\{{\hat{R}}_1>r_{12},{\hat{R}}_2>r_{21}\}$>或$>\{{\hat{R}}_1<r_{12},{\hat{R}}_3<r_{31}\},$>则判定该信号的调制方为4FSK；

若则判定该信号的调制方式为2FSK；

若则判定该信号的调制方式为BPSK；

若$>\{{\hat{R}}_1>r_{11},{\hat{R}}_2<r_{22},{\hat{R}}_5<r_{51},{\hat{R}}_6<r_{61}\},$>则判定该信号的调制方式为2ASK；

若$>\{{\hat{R}}_1>r_{11},{\hat{R}}_2<r_{22},{\hat{R}}_5<r_{51},{\hat{R}}_6>r_{61}\},$>则判定该信号的调制方式为4ASK；

若则判定该信号的调制方式2FSK*；

若则判定该信号的调制方式为4FSK*；

若$>\{{\hat{R}}_1<r_{12},{\hat{R}}_3>r_{31},{\hat{R}}_4<r_{42},{\hat{R}}_5>r_{52}\},$>则判定该信号的调制方式QPSK；

若$>\{{\hat{R}}_1<r_{12},{\hat{R}}_3>r_{31},{\hat{R}}_4<r_{42},{\hat{R}}_5<r_{52}\},$>则判定该信号的调制方式8PSK；

其中是特征R₁的估计值，r₁₁和r₁₂为对应特征R₁的阈值；其余类同。

所述的步骤(1)中信号的采样频率满足

所述的步骤(2)中每份信号的采样点数取为N＝1024，平滑窗宽度系数为M＝63。

所述的步骤(3)中小波去噪的原则是：使用sym8小波对含噪序列进行分解，在所分解的第五层上，使用软sure域值选择方法对序列进行去噪，并且该域值随第一层小波分解的噪声方差作调整。

所述的步骤(4)中的5种基于调制信号谱相关函数的调制识别特征的定义如下：

①谱相关函数α截面和f截面最大绝对值比R₁定义为：

$>R_1=\underset{N\to \infty }{\lim }\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}\frac{\max \left(\right|S_{xn}^{\alpha }\left(0\right)\left|\right)}{\max \left(\right|S_{xn}^0\left(f\right)\left|\right)}$>

其中和分别为谱相关函数的α截面和f截面函数；

②α截面强谱线数目R₂定义为：

$>R_2=\underset{N\to \infty }{\lim }\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}count({\rho }_n\ge {\rho }_{th}*{\rho }_{max})$>

其中count(·)表示求谱相关图α截面的显著峰值数，ρ_n为截面波峰的显著度，ρ_max为截面波峰的最大显著度，ρ_th为显著度阈值，取定值0.27；

其中波峰显著度ρ定义为：

$>\rho =\frac{h^2}{l*max\left(h\right)}$>

其中h为谱相关图α截面中波峰幅值与相邻两个波谷幅值中较大者之差，l为波峰的宽度值，max(h)为α截面中的最大h值；

③α截面波动系数R₃定义为：

$>R_3=\underset{N\to \infty }{\lim }\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}count({\beta }_n\ge {\beta }_{th}*{\beta }_{max})$>

其中count(·)表示统计谱相关函数α截面起伏度大于β_th*β_max的波峰数目，β_n为截面波峰的起伏度，β_max为截面波峰的最大起伏度，β_th为起伏度阈值，取定值0.1；

其中波峰起伏度β定义为：

$>\beta =\frac{h}{l}$>

其中h和l的含义同②中所述；

④f截面归一化面积R₄定义为：

$>R_4=\underset{N\to \infty }{\lim }\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}\left\{\frac{1}{max\left(S_{xn}^0\right(f\left)\right)}{\int }_{-f_0}^{f_0}\right|S_{xn}^0\left(f\right)\left|df\right\}$>

其中表示求取谱相关函数f截面函数的最大值，表示求取f截面的面积；

⑤α截面谱线的显著度比R₅定义为：

$>R_5=\underset{N\to \infty }{\lim }\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}\frac{\sec \left({\rho }_n\right)}{max\left({\rho }_n\right)}$>

其中max(ρ_n)为α截面谱线的最大显著度值，sec(ρ_n)为α截面谱线的次大显著度值；

所述的步骤(4)中的基于信号时域的统计特征零中心归一化瞬时幅度绝对值的标准偏差R₆定义为：

$>\delta =\sqrt{\frac{1}{N}\left[\underset{i=1}{\overset{N}{\Sigma }}{A}_{cn}^2\left(i\right)\right]-{\left[\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}\right|A_{cn}\left(i\right)\left|\right]}^2}$>

其中N为观测样本数，A_cn(i)为零中心归一化瞬时幅度，定义为A_cn(i)＝A_n(i)-1；

其中除数而A(i)为样本点瞬时幅度；

所述的步骤(4)中的各调试识别特征的阈值取值如下：

特征R₁的阈值：r₁₁＝0.67，r₁₂＝0.36；特征R₂的阈值：r₂₁＝4.8，r₂₂＝3.2；特征R₃的阈值：r₃₁＝60.46；；特征R₄的阈值：r₄₁＝903.51，r₄₂＝498.54；特征R₅的阈值：r₅₁＝0.073，r₅₂＝7.95×10^-4；特征R₆的阈值：r₆₁＝0.29。

本发明的原理是：循环谱相关理论是一种非平稳信号分析方法，适用于具有循环平稳性的通信调制信号的分析。针对原谱相关估计算法输出结果中由外界干扰和样本点数受限引起的随机波动现象，引入小波去噪和叠加求取平均的环节对谱相关运算过程进行优化，有效减弱了输出结果的随机波动，并使谱特征显著加强。在此基础上，可以得到调制信号的5种基于循环谱的识别特征，通过结合一种基于信号时域的统计特征，对多种数字调制信号进行识别。

本发明与现有技术的优点在于：与传统的功率谱分析相比，本发明继承循环谱分析的优点，可以充分发掘通信信号的调制特征，使调制识别过程不过分依赖于噪声敏感性较强的信号时域统计特征，从而提高低信噪比时的识别精度；与原谱相关估计算法相比，本发明引入小波去噪和叠加求取平均的环节对谱相关运算过程进行优化，从而减弱输出结果中由外界干扰和样本点数受限引起的随机波动，使谱特征得以加强，利于调制识别特征的提取；本发明提出的5个基于谱相关函数的调制识别特征较基于信号时域统计的识别特征形式简明，噪声敏感性较弱。

附图说明

图1为本发明的基于循环谱相关的调制识别方法流程图；

图2为谱相关运算之后的调制识别流程图。

具体实施方式

下面结合附图及具体实施方式详细介绍本发明。

如图1所示，本发明的基于循环谱相关的常用数字调制信号识别方法的具体实施步骤如下：

(1)由于本发明在原谱相关估计算法中引入了累加求取平均的环节，需将原始信号经采样后进行等份，采样频率设为T_S是时域采样间隔，每等份信号一般取1024个样本点；

(2)循环谱相关理论是由Gardner提出的用于非平稳信号分析的一种理论工具，其基于频域平滑的估计算法定义为：

$>S_{x\Delta t}^{\alpha }{(t,f)}_{\Delta f}=\frac{1}{M}\underset{v=-(M-1)/2}{\overset{(M-1)/2}{\Sigma }}\frac{1}{\Delta t}{X}_{\Delta t}(t,f+\frac{\alpha }{2}+{\mathrm{vF}}_S)X_{\Delta t}^{*}(t,f-\frac{\alpha }{2}+{\mathrm{vF}}_S)---\left(1\right)$>

上式中x_k(t)为所截获信号的n等份之一；Δt为x_k(t)持续时间；Δf为频域平滑间隔；F_S为频域最小增量单位；T_S是时域采样间隔；M为平滑窗宽度系数，取定值为63；α为循环频率；f为频率；t为时间；X_Δt(t,f)为步骤(1)中等份信号X_Δt(t)的短时傅里叶变换，定义为：

$>X_{\Delta t}(t,f)=\underset{K=0}{\overset{N-1}{\Sigma }}{a}_{\Delta t}\left({\mathrm{KT}}_S\right)x_k(t-{\mathrm{KT}}_S)\exp (-j2\pi f(t-{\mathrm{KT}}_S\left)\right)---\left(2\right)$>

上式中a_Δt是窗函数，此处取为与等份信号等长的矩形窗函数；短时傅里叶变换的点数取为N＝1024；

(3)由于外界干扰和样本点数受限的影响，步骤(2)中的谱相关计算结果波动较大，不利于弱特征的观测和提取，故先对上述n份等份信号的谱相关函数的循环频率α进行小波去噪，再相加求取平均，如下式所示：

$>S_{x\Delta t}^{\alpha }\left(f\right)=\frac{1}{N}\underset{k=1}{\overset{N}{\Sigma }}{\mathrm{WDN}}_{\alpha }\left[S_{x\Delta t}^{\alpha }(t_k,f)\right]$>

(3)

式中WDN_α(·)表示对中的循环频率α进行小波去噪，如下所示：

y＝wden(x,'heursure','s','sln',5,'sym8')(4)

其中x和y为去噪前后的离散序列，整个式子表示使用sym8小波对含噪信号序列进行分解，在所分解的第五层上，使用软sure域值选择方法对序列进行去噪，并且该域值随第一层小波分解的噪声方差作调整。

(4)通过分析调制信号的谱相关函数可以发现，谱特征主要集中在谱相关图的α截面和f截面上，因此可以通过提取这两个截面的谱特征组合成调试识别特征来对调制信号进行识别，所提出的特征为谱相关函数α截面和f截面最大绝对值比R₁、α截面强谱线数目R₂、α截面波动系数R₃、f截面归一化面积R₄、α截面谱线的显著度比R₅；另外，该调制识别过程还涉及一种信号时域统计特征——零中心归一化瞬时幅度绝对值的标准偏差R₆；为了准确描述调制识别特征，引入截面波峰起伏度β和显著度ρ的定义，这两个定义分别用于描述α截面的背景抖动和谱线强度，

其中波峰起伏度β定义为：

$>\beta =\frac{h}{l}---\left(5\right)$>

其中h为谱相关图α截面中波峰幅值与相邻两个波谷幅值中较大者之差，l为波峰的宽度值；

波峰显著度ρ定义为：

$>\rho =\frac{h^2}{l*max\left(h\right)}---\left(6\right)$>

其中max(h)为α截面中的最大h值；在此基础上对各调制识别特征作描述如下：

①谱相关函数α截面和f截面最大绝对值比R₁定义为：

$>R_1=\underset{N\to \infty }{\lim }\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}\frac{\max \left(\right|S_{xn}^{\alpha }\left(0\right)\left|\right)}{\max \left(\right|S_{xn}^0\left(f\right)\left|\right)}---\left(7\right)$>

其中和分别为谱相关函数的α截面和f截面函数；该特征的阈值取为r₁₁＝0.67，r₁₂＝0.36；

②α截面强谱线数目R₂定义为：

$>R_2=\underset{N\to \infty }{\lim }\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}count({\rho }_n\ge {\rho }_{th}*{\rho }_{max})$>

其中count(·)表示求谱相关图α截面的显著峰值数，ρ_n为截面波峰的显著度，ρ_max为截面波峰的最大显著度，ρ_th为显著度阈值，取定值0.27；该特征的阈值取为r₂₁＝4.8，r₂₂＝3.2；

③α截面波动系数R₃定义为：

$>R_3=\underset{N\to \infty }{\lim }\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}count({\beta }_n\ge {\beta }_{th}*{\beta }_{max})$>

其中count(·)表示统计谱相关函数α截面起伏度大于β_th*β_max的波峰数目，β_n为截面波峰的起伏度，β_max为截面波峰的最大起伏度，β_th为起伏度阈值，取定值0.1；该特征的阈值取为r₃₁＝60.46；

④f截面归一化面积R₄定义为：

$>R_4=\underset{N\to \infty }{\lim }\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}\left\{\frac{1}{max\left(S_{xn}^0\right(f\left)\right)}{\int }_{-f_0}^{f_0}\right|S_{xn}^0\left(f\right)\left|df\right\}$>

其中表示求取谱相关函数f截面函数的最大值，表示求取f截面的面积；该特征的阈值取为r₄₁＝903.51，r₄₂＝498.54；

⑤α截面谱线的显著度比R₅定义为：

$>R_5=\underset{N\to \infty }{\lim }\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}\frac{\sec \left({\rho }_n\right)}{max\left({\rho }_n\right)}$>

其中max(ρ_n)为α截面谱线的最大显著度值，sec(ρ_n)为α截面谱线的次大显著度值；该特征的阈值取为r₅₁＝0.073，r₅₂＝7.95×10^-4

⑥零中心归一化瞬时幅度绝对值的标准偏差R₆定义为：

$>\delta =\sqrt{\frac{1}{N}\left[\underset{i=1}{\overset{N}{\Sigma }}{A}_{cn}^2\left(i\right)\right]-{\left[\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}\right|A_{cn}\left(i\right)\left|\right]}^2}$>

其中N为观测样本数，A_cn(i)为零中心归一化瞬时幅度，定义为A_cn(i)＝A_n(i)-1；其中除数而A(i)为样本点瞬时幅度；该特征的阈值取为r₆₁＝0.29。

结合上述6种调制识别特征对几种数字调制信号进行识别，具体的调制识别流程如图2所示，

若则判定该信号的调制方式为MSK；

若$>\{{\hat{R}}_1>r_{12},{\hat{R}}_2>r_{21}\}$>或$>\{{\hat{R}}_1<r_{12},{\hat{R}}_3<r_{31}\},$>则判定该信号的调制方式为4FSK；

若则判定该信号的调制方式为2FSK；

若则判定该信号的调制方式为BPSK；

若$>\{{\hat{R}}_1>r_{11},{\hat{R}}_2<r_{22},{\hat{R}}_5<r_{51},{\hat{R}}_6<r_{61}\},$>则判定该信号的调制方式为2ASK；

若$>\{{\hat{R}}_1>r_{11},{\hat{R}}_2<r_{22},{\hat{R}}_5<r_{51},{\hat{R}}_6>r_{61}\},$>则判定该信号的调制方式为4ASK；

若则判定该信号的调制方式为2FSK*；

若则判定该信号的调制方式为4FSK*；

若$>\{{\hat{R}}_1<r_{12},{\hat{R}}_3>r_{31},{\hat{R}}_4<r_{42},{\hat{R}}_5>r_{52}\},$>则判定该信号的调制方式QPSK；

若$>\{{\hat{R}}_1<r_{12},{\hat{R}}_3>r_{31},{\hat{R}}_4<r_{42},{\hat{R}}_5<r_{52}\},$>则判定该信号的调制方式8PSK；

其中为调制识别过程中特征R的估计值。

根据这样的识别准则，在信噪比为10dB时各数字调制方式的识别率均在88.8％以上，识别效果比较理想，所识别调制方式包括2ASK、4ASK、BPSK、QPSK、8PSK、MSK、2FSK、4FSK、2FSK*、4FSK*，其中FSK*表示码元初始相位不相关的频移键控。

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

尽管为说明目的公开了本发明的最佳实例和附图，但是本领域的技术人员可以理解：在不脱离本发明及所附的权利要求的精神和范围内，各种替换、变化和修改都是可能的。因此，本发明所保护的技术方案不应局限于最佳实例和附图所公开的内容。