本发明涉及通信技术领域，尤其涉及一种无需先验知识的常用通信信号调制方式的自动识别方法。
背景技术：
通信信号调制方式自动识别，指的是在缺少先验条件时对接收机截取信号的调制模式做出判断，并进一步估计信号的调制参数。自上世纪90年代以来，无线通信技术迅速发展，通信调制模式复杂多样使得通信信号调制模式识别成为通信系统及信号处理领域的热点问题。调制模式自动识别作为信号接收和解调的中间步骤，被广泛应用于军事和商用的非协作通信领域之中。
在民用方面，主要是为频谱管理进行信号身份确认、干扰确认等，正确识别收发双方的调制样式，是保证信息无误转发的基本条件。商业应用中，通信信号的自动调制识别技术可以应用在信号证实、干扰识别、频率管理等方面。军用领域中，调制方式的自动识别技术主要应用在电子战争、威胁评估等方面。电子战争作为现代化战争的一部分要求识别系统在密集的电磁环境下做出快速的反应，信号调制模式的识别成为不可缺少的一部分。电子对抗战需要截获敌方信号并进行侦破，从而获取敌方信息。由于非协作通信，接收机截取的通常是全盲信号，首先需要识别出信号的调制类型，在此基础上才能进一步完成信号解调和恢复工作，获取信息并破坏敌方的通信。在早期的通信情报系统中，采用手动调制方式识别，即依靠操作员解释测量到的参数，判断信号调制方式。随着通信技术的飞速发展，通信信号的体制和调制样式变得更加复杂多样，信号环境日趋密集，这使得通信信号的识别变得越来越复杂、也越来越困难，尤其是在无任何先验知识的情况下。调制方式的差异性也很小，因此，如何能够在信噪比较低的情况下进行信号特征参数的有效提取和多类调制信号的调制识别是至关重要的问题。
技术实现要素：
根据现有技术存在的问题，本发明公开了一种无需先验知识的常用通信信号调制方式的自动识别方法，具体包括以下方案包括以下步骤：
S1：接收待处理信号数据s(n)；
S2：提取信号瞬时特征参数信息，其中信号瞬时特征参数包括瞬时幅度a(n)、瞬时相位和瞬时频率f(n)；
S3：将信号截取成多段子信号，采用小波变换方法对每段子信号提取码速率、并计算码速率的方差判断方差的值是否大于门限如果大于门限那么该数据为模拟调制信号，进入S4；如果小于门限那么判断为数字调制信号，并进入S7；
S4：对信号瞬时幅度a(n)计算γmax1参数，判断γmax1参数的值是否大于门限t(γmax1)，如果大于门限值t(γmax1)，进入S5，否则进入S6；
S5：对瞬时非线性相位计算参数σdp，判断σdp参数是否大于门限值t(σdp)，如果大于门限值t(σdp)，判断为DSB信号，否则判断为AM信号；
S6：计算信号s(n)的对称性参数P，判断参数P是否大于门限t(P)，如果大于门限值，判断为SSB信号，否则判断为FM信号；
S7：对信号瞬时幅度a(n)计算γmax2参数，判断γmax2参数的值是否大于门限t(γmax2)，如果大于门限t(γmax2)，进入S8，否则进入S11；
S8：对信号s(n)计算参数，判断参数的值是否大于门限如果大于门限进入S9，否则进入S10；
S9：计算信号瞬时幅度a(n)的幅度种类个数M，如果M≤3，判决为2ASK；如果4≤M，判决为4ASK；
S10：计算信号瞬时幅度a(n)的幅度种类个数L，如果3≤L≤6，判决为16QAM；7≤L，判决为64QAM；
S11：对信号瞬时频率f(n)计算参数，判断参数的值是否大于门限如果大于门限进入S12，否则进入S13；
S12：计算信号的二相特征参数判断参数的值是否大于门限如果大于门限判断为2PSK信号，否则判断为4PSK信号；
S13：计算信号功率谱的峰值个数N，如果N<2，判决为CW，如果2≤N<3判决为2FSK，其他情况下判决为4FSK。
所述码速率的方差采用如下方式计算：
(1)对接收信号进行采样并截取多段子信号，得到s1(n),s2(n)…sM(n)；
(2)对截取的子信号进行非线性变换得到Q1(n),Q2(n)…QM(n)；
Q(i)＝s2(i)-s(i-1)s(i+1)，
(3)对Q1(n),Q2(n)…QM(n)进行小波变换，其中小波变换采用小波函数为haar基函数，其数学表达式为：
其傅里叶变换为满足允许条件：
且ψ(t)满足一般函数的约束条件：
(4)将母函数ψ(t)经过平移和伸缩后得：
称其为一个小波序列，其中，a为伸缩因子，b为平移因子；
(5)对ψ(t)进行离散小波变换，对于任意的函数f(t)∈L2(R)的连续小波变换为：
离散小波变换即将变换中尺度参数a和平移参数b分别取做离散形式：
(6)重复步骤(3)～(5)，进行两次小波变换，得到WQ1(n),WQ2(n)…WQM(n)；
(7)对WQ1(n),WQ2(n)…WQM(n)进行傅里叶变换，得到Y1(ω),Y2(ω)…YM(ω)；
(8)考虑Y1(ω),Y2(ω)…YM(ω)正频率，Ym(ω),m∈[1,2...M]的第一个尖峰所在位置即为码速率Rsm,m∈[1,2...M]；
(9)计算Rs1,Rs2…RsM的方差
所述信号瞬时特征参数信息采用如下方式获取：
(1)对于连续波调制，已调实信号的表达式为：
式中，ωc表示载波的角频率，调制信号存放在已调信号的振幅a(t)、频率ωc和相位中；
(2)对信号进行希尔伯特变换，得到信号的正交分量：
获得信号s(t)的解析表达式z(t)＝s(t)+jv(t)；
(3)提取瞬时幅度，信号包络瞬时幅度估计为：
(4)提取瞬时相位，瞬时相位估计为：
(5)对瞬时相位进行去相位折叠，计算修正相位序列C(i)：
在求得瞬时相位后对相位去折叠、恢复出没有折叠的相位，得到无折叠的相位为：
(6)对去相位折叠后的相位序列去载频线性相位，得到非线性相位分量
信号的相位包括：载频引起的线性相位分量和调制相位的非线性相位分量非线性相位分量由下式计算：
其中在载频fc不能精确已知的情况下，采用相性规划法估计出线性相位分量并使误差最小，可求出C1和C2两个常数；
(7)提取瞬时频率，瞬时频率可以由瞬时相位的导数计算求得，对数字信号可直接求差分，即：
由于采用了上述技术方案，本发明提供的一种无需先验知识的常用通信信号调制方式的自动识别方法，方法中提取信号瞬时特征参数信息对信号进行自动识别，该方法中的信号既有数字调制信号、又有模拟调制信号，算法包含的调制信号种类丰富，适用范围更广，其中处理数据和算法的判决规则简单并且收敛，能够快速得到信号的识别结果，可在线分析，并且该算法对特征参数，如载频、瞬时幅度、相位、频率不敏感，所提的判决参数具有较高的稳定性。
附图说明
为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明中方法的流程图。
具体实施方式
为使本发明的技术方案和优点更加清楚，下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚完整的描述：
通信领域中我们把信号分为两种类型即：
模拟调制：AM，FM，DSB，SSB
数字调制：2ASK，4ASK，2PSK，4PSK，2FSK，4FSK，16QAM，64QAM；
设多段信号的码速率的方差为数字调制信号具有恒定的码速率，而模拟调制信号则无此特征。因此可以截取多段信号，提取每段信号的码速率。数字调制信号的码速率序列的方差理论上为而模拟调制信号的方差不为零，因此可设定门限区分两类调制信号。门限用表示，当时，判断为模拟调制信号，否则判断为数字调制信号。
进一步的，模拟调制信号的区分为：
1、零中心归一化瞬时幅度谱密度的最大值γmax1
γmax1＝max|FFT(acn(i))|2/Ns
式中，Ns为取样点数，acn(i)为零中心归一化瞬时幅度，acn(i)＝a(i)/mean(a)-1。该参数主要用来区分信号有无幅度变换。在模拟调制信号中，FM信号和SSB信号无幅度调制，而DSB信号和AM信号有幅度调制。故可以设置门限区别有无幅度调制。
门限用t(γmax1)表示，当γmax1>t(γmax1)时，判断为DSB或者AM信号，否则判断为FM或者SSB信号。
2、谱对称参数的绝对值|P|
式中，其中S(i)＝FFT(s(n))是信号的傅里叶变换。fcn＝fc*Ns/fs-1，其中fc为载频，fs为采样频率，Ns为采样点数。
特征参数P是信号频谱对称性的量度，用来区分FM信号与SSB信号。FM信号其频谱满足对称性，其P值很小，在无噪声的理想情况下P＝0，而SSB信号其频谱满足不对称性，在无噪声情况下，SSB信号只有单边带。该参数用来识别信号的频谱相对载频是否对称。故可通过设置门限区别两类信号。
门限用t(|P|)表示，当|P|>t(|P|)时，判断为SSB信号，否则判断为FM信号。
3、非弱信号段零中心瞬时相位非线性分量的方差σdp
该参数主要用来区别相位无变化的AM类信号和相位有变化的DSB类信号。故可通过设置门限区别两类信号。
门限用t(σdp)表示，当σdp>t(σdp)时，判断为DSB信号，否则判断为AM信号。
进一步的，数字调制信号的区分：
1、零中心归一化瞬时幅度谱密度的最大值γmax2
γmax2＝max|FFT(acn(i))|2/Ns
式中，Ns为取样点数，acn(i)为零中心归一化瞬时幅度，acn(i)＝a(i)/mean(a)-1。该参数主要用来区分信号有无幅度变换。在数字调制信号中，MPSK信号和MFSK信号无幅度调制，而MASK信号和MQAM信号有幅度调制。故可以设置门限区别有无幅度调制。
门限用t(γmax2)表示，当γmax2>t(γmax2)时，判断为MASK或者MQAM信号，否则判断为MFSK或者MPSK信号。
2、归一化零中心瞬时频率的紧致性
该参数主要用来区别MPSK信号还是MFSK信号，因为对MPSK信号，其瞬时频率为常值，具有较高的紧致性，即该参数较大，而MFSK信号只有两个或者四个电平，其紧致性较差。故可通过设置门限区别两类信号。
门限用表示，当时，判断为MPSK信号，否则判断为MFSK信号。
3、归一化零中心信号的紧致性
该参数主要用来区别MASK信号还是MQAM信号。scn(i)的定义和acn(i)的定义相同。故可通过设置门限区别两类信号。
门限用表示，当时，判断为MASK信号，否则判断为MQAM信号。
MPSK\MFSK\MASK\MQAM信号的各自识别方式
MPSK归一化二相特征参数
定义信号x循环相关函数为：其傅里叶变换为信号的循环谱密度，又称谱相关密度函数。在对MPSK信号进行四次方之后，求得的谱相关函数可以识别不同调制方式的信号。定义归一化二相特征函数：
和均是信号四次方之后的谱相关密度。BPSK信号的归一化二相特征参数为QPSK信号的归一化二相特征参数为因此可根据此识别BPSK和QPSK信号。
门限用表示，当时，判断为BPSK信号，否则判断为QPSK信号。
MFSK功率谱的峰值个数N
2FSK信号的功率谱有两个频率分量，因此有两个峰值；4FSK信号的功率谱有四个频率分量，因此有四个峰值。为了准确提取谱峰的个数，用自回归模型(AR)进行谱密度估计。
首先通过功率谱估计，得到功率谱的光滑曲线，以准确提取正弦分量所形成的尖峰；其次得到功率谱序列的一阶差分序列；最后通过判断差分序列的符号变化，来判断谱峰的个数N。
设AR模型为：
其中p是AR模型的阶数，ak是模型的系数，u(n)是一个均值为0，方差为σ2的白噪声序列。由AR模型得到的功率谱的估计为
有了MFSK信号的功率谱估计序列，通过分析它的一阶差分序列符号变化，就可以检测出谱峰。记一阶差分序列为如果d(n)>0并且d(n+1)<0，则判定在该频率处存在谱峰。并由此得到整个功率谱上的峰值个数N。
设定如下判决规则：N<2判决为CW，2≤N<3判决为2FSK，其他情况下判决为4FSK。
但是上述结论是在功率谱估计结果没有波纹起伏的理想情况下的结论，实际中估计的功率谱即使非常平滑，也势必存在起伏，尤其是在频谱中的弱信号部分。为排除弱信号起伏的影响，我们在对差分序列进行符号判决之前，首先设立门槛值λ(Px)，即功率谱估计序列中，只有信号能量大于门槛值，其相应的差分序列值才能参与谱峰个数判决。
MASK\MQAM信号瞬时幅度的个数M/L
MASK与MQAM信号幅度提取方法一致，此处以MASK信号为例。
对于MASK信号，其时域上的特点是较为明显的，即就L个码元时隙而言，有M种振幅。因此，本文采用的识别MASK信号幅度调制特点的特征参数为：MASK信号在L个码元时隙内的振幅种类。所采用的方法是：计算每个码元内采样值的绝对值和，并相对于最大的绝对值和取归一化，检测这L个码元数值的种类。幅度归一化后的取值范围可以为综上所述，即对于MASK信号而言，可提取出幅值的个数为M。
设立判决：
MASK：M≤3，判决为2ASK；4≤M，判决为4ASK。
MQAM：3≤L≤6，判决为16QAM；7≤L，判决为64QAM。
本发明以实际通信系统中最常用的AM，FM，DSB，SSB，ASK，4ASK，2PSK，4PSK，2FSK，4FSK，16QAM，64QAM等信号为例，进行无需先验知识的常用通信信号调制方式的自动识别计算机仿真，以随机序列作为数字调制信号的信源。用均值为0的正态分布随机数发生器产生加性高斯白噪声AWGN序列。SNR变化范围从6dB到24dB，步进为2dB。每个样本仿真400次。仿真平台操作系统是WINDOWS7，采用MATLAB进行编程。
其中用于仿真的信号调制参数：
数字调制信号：
载频：10kHz，采样率：100kHz，采样时间：200ms，码速率：RB＝500(Baud)；2FSK频率：5kHz，10kHZ；4FSK频率：5kHz，10kHZ，15kHz，20kHZ.
模拟调制信号：
AM：s＝[1+ma cos(2πfmt)]cos(2πf0t)，fm＝100Hz，ma＝0.6
FM：s＝cos[2πf0t+βsin(2πfmt)]，β＝1，fm＝100Hz
DSB：s＝ma cos(2πfmt)cos(2πf0t)，fm＝100Hz，ma＝0.6
SSB：s＝ma[cos(2πfmt)cos(2πf0t)+sin(2πfmt)sin(2πf0t)]，fm＝500Hz，
ma＝0.6
仿真中个判决的门限值是基于大量仿真统计得到的，具体设定如下：
码速率方差判决门限t(Rsv)＝0.05；
瞬时幅度a(n)参数γmax1的门限t(γmax1)＝10dB；
信号s(n)的对称性参数P的门限t(|P|)＝0.6；
瞬时非线性相位参数σdp的门限t(σdp)＝6；
瞬时幅度a(n)参数γmax2的门限t(γmax2)＝10dB；
信号s(n)参数的门限
信号瞬时频率f(n)参数的门限
如图1所示的一种无需先验知识的常用通信信号调制方式的自动识别方法，包括以下步骤：
1)接收待处理信号数据s(n)；
2)提取信号瞬时特征参数：瞬时幅度a(n)、瞬时相位瞬时频率f(n)；
3)将信号截取成多段，利用小波变换方法对每段子信号提取码速率，并计算码速率的方差判断的值是否大于门限如果大于门限，那么该数据为模拟调制信号，进入第4步；如果小于门限，那么判断为数字调制信号，并进入第7步；(由于仅在此步将信号数据截取为多段，其他参数的计算均通过完整的信号数据计算)；
4)对信号瞬时幅度a(n)计算γmax1参数，判断γmax1参数的值是否大于门限t(γmax1)，如果大于门限，进入第5步，否则进入第6步；
5)对瞬时非线性相位计算参数σdp，判断σdp参数是否大于门限t(σdp)，如果大于门限值，判断为DSB信号，否则判断为AM信号；
6)计算信号s(n)的对称性参数P，判断参数P是否大于门限t(P)，如果大于门限值，判断为SSB信号，否则判断为FM信号；
7)对信号瞬时幅度a(n)计算γmax2参数，判断γmax2参数的值是否大于门限t(γmax2)，如果大于门限，进入第8步，否则进入第11步；
8)对信号s(n)计算参数，判断参数的值是否大于门限如果大于门限，进入第9步，否则进入第10步；
9)计算信号瞬时幅度a(n)的幅度种类个数M，如果M≤3，判决为2ASK；如果4≤M，判决为4ASK；
10)计算信号瞬时幅度a(n)的幅度种类个数L，如果3≤L≤6，判决为16QAM；7≤L，判决为64QAM；
11)对信号瞬时频率f(n)计算参数，判断参数的值是否大于门限如果大于门限，进入第12步，否则进入第13步
12)计算信号的二相特征参数判断参数的值是否大于门限如果大于门限，判断为2PSK信号，否则判断为4PSK信号；
13)计算信号功率谱的峰值个数N，如果N<2，判决为CW，如果2≤N<3判决为2FSK，其他情况下判决为4FSK。
所述的码速率的方差计算方法为：
(1)对接收信号进行采样，并截取多段，得到s1(n),s2(n)…sM(n)；
(2)对截取的信号进行非线性变换：
Q(i)＝s2(i)-s(i-1)s(i+1)，
得到Q1(n),Q2(n)…QM(n)，对于单频正弦波信号(CW)输入，其输出为一个恒量；对于已调信号的输入，其输出序列值发生变化，序列取值变化的特征反映了信号调制的信息；对于数字调制信号，该变换能够提取出码元的变化位置；这种非线性变换的特点是，对于单频正弦波信号(CW)输入，其输出为一个恒量。对于已调信号的输入，其输出序列值发生变化，序列取值变化的特征反映了信号调制的信息。对于数字调制信号，该变换能够提取出码元的变化位置。之后利用小波变换，能够有效提取信号的码速率。
(3)对Q1(n),Q2(n)…QM(n)进行小波变换。本发明采用的小波函数为haar基函数，其数学表达式为：
其傅里叶变换为满足允许条件(完全重构条件或恒等分辨条件)：
且ψ(t)满足一般函数的约束条件：
(4)将母函数ψ(t)经过平移和伸缩后得：
称其为一个小波序列。其中，a为伸缩因子，b为平移因子；
(5)对ψ(t)进行离散小波变换，对于任意的函数f(t)∈L2(R)的连续小波变换为：
离散小波变换即将变换中尺度参数a和平移参数b分别取做离散形式：
(6)重复步骤(3)～(5)，进行两次小波变换，得到WQ1(n),WQ2(n)…WQM(n)；
(7)对WQ1(n),WQ2(n)…WQM(n)进行傅里叶变换，得到Y1(ω),Y2(ω)…YM(ω)；
(8)考虑Y1(ω),Y2(ω)…YM(ω)正频率，Ym(ω),m∈[1,2...M]的第一个尖峰所在位置即为码速率Rsm,m∈[1,2...M]；
(9)计算Rs1,Rs2…RsM的方差
所述的瞬时特征参数的计算方法为：
(1)尽管信号的调制样式多种多样，但实质上调制是用调制信号去控制载波的某一个或某几个参数，使这个参数按照调制信号的规律而变化的过程。载波可以是正弦波或者脉冲序列，以正弦波信号作为载波的调制叫做连续波调制。对于连续波调制，已调实信号的一般表达式为：
式中，ωc表示载波的角频率。调制信号可以分别“寄生”在已调信号的振幅a(t)、频率ωc和相位中。
(2)对信号进行希尔伯特变换，希尔伯特变换实际是对原信号进行了π/2相移，使它成为自身的正交对，得到信号的正交分量：
可得信号s(t)的解析表达式z(t)＝s(t)+jv(t)。
(3)提取瞬时幅度，信号包络瞬时幅度估计为：
(4)提取瞬时相位，瞬时相位估计为：
(5)对瞬时相位进行去相位折叠，计算修正相位序列C(i)：
由于实际相位计算是是以模2π来计算的，因此在求得瞬时相位之后需要对相位去折叠，恢复出没有折叠的相位，得到无折叠的相位为：
(6)对去相位折叠后的相位序列去载频线性相位，得到非线性相位分量
信号的相位包括两部分，载频引起的线性相位分量和调制相位的非线性相位分量因此非线性相位分量可由下式计算：
在载频fc不能精确已知的情况下，可以采用相性规划法估计出线性相位分量并使误差最小，则可以求出C1和C2两个常数，所以有：
提取瞬时频率，瞬时频率可以由瞬时相位的导数计算求得，对数字信号可直接求差分，即：
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。