CN104022837B - 一种基于特征提取的信号调制方式识别方法 - Google Patents

Abstract

本发明公开了一种基于特征提取的信号调制方式识别方法，对于接收到的信号，首先计算其归一化功率谱以及归一化功率谱中大于常数B时谱峰的个数K；然后根据K值、归一化零中心瞬时幅度之谱密度的最大值γ_max、圆环特征参数num的大小来综合识别通信信号的调制方式。本发明所述的方法，能够快速、准确地对通信信号的调制方式进行识别。

Description

一种基于特征提取的信号调制方式识别方法

技术领域

本发明涉及通信技术领域，特别涉及一种基于特征提取的信号调制方式识别方法。

背景技术

通信信号调制识别的作用就是对不明信号的调制方式进行分析、判决和归类以及对其调制参数进行估计，可广泛地运用于民用及军用频谱监测领域。调制识别的研究最近二十多年取得了迅速的发展，科研人员研究出了一批有价值的识别方法。文献[1-4]分析了通信信号在自由传播空间当中的传输特性，并提取一些特征参数(如γ_max、δ_dp等参数)，这些特征参数对卫星通信或近距离通信的通信信号有较好的识别效果。由于实际通信信号的传播大多经历多径衰减，接收机接收到的信号的瞬时幅度、瞬时相位短时间内可能产生剧烈的变化，在此环境下利用瞬时幅度、瞬时相位等特征参数对信号调制方式识别可能会产生误导，甚至是完全无效。相对而言，信号的频谱特征仍然保持较好的特性，许多科研人员对信号的频谱特性进行了研究并得到了一系列有意义的性质。文献[5]当中证明了MFSK的频谱在多径条件下仍然能较好地保持其谱峰特征；文献[6]基于数字调制的循环平稳特性，证明了BPSK的平方谱在二倍载频处存在明显的单频分量；OQPSK的平方谱在处存在离散谱线(注：实际上谱线明显与否与发送滤波器的滚降因子有关，f代表频率，T为发送符号的符号周期)；QPSK的四次方谱在其四倍载频处存在明显单频分量；π/4-QPSK的四次方谱在存在明显的离散谱线。理论上MPSK的M次方谱应该在其M倍处存在明显的单频分量，但文献[7]中提到信号M次方后，信噪比会有明显降低。因此在实际运用当中，我们最多考虑信号的4次方谱。此外，运用高阶累积量来区分数字调制的研究也比较多，文献[8]当中所提到的特征参数F_x可以较好地区分MPSK和MQAM，并且具有较好的抗噪声能力。

本专利对上述所提到的方法加以改进并结合自己设计的参数得到一种新的调制识别方法。

参考文献

[1]E.E.Azzouz,A.K.Nandi,.Automatic Modulation recognition ofcommunication signals[M].Kluwer.Acadmic Publishers,1996.

[2]范海波,杨志俊,曹志刚,.卫星通信常用调制方式的自动调制识别[J].通信学报,2004,25(1):140-149.

[3]李杨,李国通,杨根庆,.通信信号数字调制方式自动识别算法研究[J].电子与信息学报,2005,27(2):197-201.

[4]A.Kubankova,D.Kubanek,.Extended method of digital modulationrecognition and its testing[J].RADIOENGINEERING,2011,20(1):25-30.

[5]崔伟亮,.衰弱信道下非合作接受中调制识别技术研究[D].郑州.解放军信息工程大学博士学位论文,2011.

[6]杨琳,.数字通信信号调制方式自动识别技术研究[D].合肥.中国科学技术大学博士学位论文，2008.

[7]廖红舒,.通信侦察信号处理关键技术研究[D].成都.电子科技大学博士学位论文,2011.

[8]陈卫东,.数字通信信号调制识别算法研究[D].西安.西安电子科技大学博士学位论文,2001.

发明内容

本发明的目的在于克服现有技术的缺点与不足，提供一种基于特征提取的信号调制方式识别方法。

本发明的目的通过以下的技术方案实现：

一种基于特征提取的信号调制方式识别方法，包含以下顺序的步骤：

1)对于接收到的信号，计算其归一化功率谱以及归一化功率谱中大于常数B时谱峰的个数K，其中B的范围为(-6,-4)；

2)当K为1时，则判断该信号的调制方式为CW、ASK、AM中的一种：首先计算其归一化零中心瞬时幅度之谱密度的最大值γ_max，若γ_max小于2则为CW，若γ_max大于2则为ASK或AM；然后计算对数幅度谱的峰度系数β，若β小于4则为ASK，若β大于4则为AM；

3)当K不为1且γ_max小于2、圆环特征参数num小于8时，则该信号的调制方式为MFSK或FM：计算其凹凸参数P，若P小于门限t₁则为MFSK,若P大于门限t₁则为FM，其中t₁的范围为(6，9)；

4)当K不为1且不满足γ_max小于2、num小于8，则该信号的调制方式被判断为BPSK、QPSK、OQPSK、Π/4-QPSK、MQAM中的一种：首先计算其归一化平方功率谱，若大于门限t₂时只有一条谱线则为BPSK，若大于门限t₂时有两条谱线则为OQPSK；然后计算其归一化四次方功率谱，若大于门限t₃时只有一条谱线则为Π/4-QPSK；接着计算其高阶累计值F，若F小于16.5则为QPSK，否则为MQAM，其中t2、t3的范围为(0.3，0.45)。

所述的步骤1)，具体包含以下步骤：

a、令 $\left\{\begin{array}{c}{e}_0^f\left(n\right)=x\left(n\right)\\ e_0^b\left(n\right)=x\left(n\right)\end{array}\right.,$ 利用下式计算：

${\hat{k}}_m=\frac{-2\underset{n=m}{\overset{N-1}{\Σ}}{e}_{m-1}^f\left(n\right)e_{m-1}^{b*}(n-1)}{\underset{n=m}{\overset{N-1}{\Σ}}|e_{m-1}^f\left(n\right){|}^2+\underset{n=m}{\overset{N-1}{\Σ}}|e_{m-1}^b(n-1){|}^2},m=1,2,...,p$

这里p我们取10，x(n)＝I(n)+iQ(n)；

b、利用

${\hat{r}}_x\left(0\right)=\frac{1}{N}\underset{n=0}{\overset{N-1}{\Σ}}|x\left(n\right){|}^2$

得到m＝1时的参数：

${\hat{a}}_1\left(1\right)={\hat{k}}_1,{\mathrm{\ρ}}_1=(1-|{\hat{k}}_1{|}^2){\hat{r}}_x\left(0\right)$

c、由和下式求出和再依据式(1)估计其中

$\left\{\begin{array}{c}{e}_m^f\left(n\right)=e_{m-1}^f\left(n\right)+k_m{e}_{m-1}^b(n-1)\\ e_m^b\left(n\right)=e_{m-1}^b\left(n\right)+k_m^{*}{e}_{m-1}^f\left(n\right),m=1,2,...,p\\ e_0^f\left(n\right)=e_0^b\left(n\right)=x\left(n\right)\end{array}\right.$

d、依据下式的递推关系，求出m＝2时的及

$\left\{\begin{array}{c}{\hat{a}}_m\left(k\right)={\hat{a}}_{m-1}\left(k\right)+{\hat{k}}_m{\hat{a}}_{m-1}^{*}(m-k),k=1,2,...,m-1\\ {\hat{a}}_m\left(m\right)={\hat{k}}_m\\ {\widehat{\mathrm{\ρ}}}_m=(1-|{\hat{k}}_m{|}^2){\widehat{\mathrm{\ρ}}}_{m-1}\end{array}\right.$

e、重复上面的过程，直到m＝p，从而求得

${\hat{a}}_p\left(k\right),k=1,2,...,p$

f、将上式所求得的代入到下式当中：

$P\left(\mathrm{\ω}\right)=\frac{1}{|1+\underset{k=1}{\overset{p}{\Σ}}{a}_k\left(cos\right(\mathrm{\ω}k)-i\·sin\left(\mathrm{\ω}k\right)){|}^2}$

上式a_k分别对应着前面我们所得到的

g、对上式的ω在[-π,π)取512个等间隔点，注意我们的取值是以下的离散值：

h、对所得到的每一个值P(ω)先归一化然后再进行dB转换，即

PW(ω)＝10log₁₀|P(ω)/max(P(ω))|

然后在这些点当中统计谱峰个数；此时所得的功率谱是利用自回归模型估计信号的平滑功率谱，即AR功率谱。

步骤2)中，所述的γ_max计算式如下：

${\mathrm{\γ}}_{max}=\frac{max|fft\[a_{cn}\left(i\right)\]|}{N_s}$

其中N_s为每一帧数据的个数，a_cn为零中心归一化瞬时幅度；

a_cn由下式计算得到：

a_cn(i)＝a_n(i)-1

其中 $a_n\left(i\right)=\frac{a\left(i\right)}{m_a},m_a=\frac{1}{N_s}\underset{i=1}{\overset{N_s}{\Σ}}a\left(i\right).$

步骤2)中，所述的β由下式计算：

$\mathrm{\β}=\frac{\frac{1}{n}\underset{i}{\Σ}{(v_i-\stackrel{\‾}{v})}^4}{{\left(\frac{1}{n}\underset{i}{\Σ}{(v_i-\stackrel{\‾}{v})}^2\right)}^2},$

其中v_i代表信号频谱第i个点的电平值，单位为dB。

步骤3)中，所述的P由下式计算：

$P=\frac{P_{21}+P_{22}}{P_1},$

假设在区域[f₁,f₂]中，f₁、f₂表示频率值，P₁是指信号载波附近区域中的电平最大值，P₂₁和P₂₂分别表示f₁右边区域的最大值和f₂左边区域的最大值。

步骤3)中，所述的圆环特征参数num为统计归一化IQ平面内以原点为中心，0.1为半径的点的个数。

所述的基于特征提取的信号调制方式识别方法，在步骤1)之前还包括对接收到的信号进行预处理，具体如下：首先对接收到的信号进行中值滤波，利用中值滤波所得信号进行带宽估计，并根据所得带宽对信号进行线性滤波处理。预处理的目的是提高接收信号的信噪比，方便后续处理。

本发明与现有技术相比，具有如下优点和有益效果：

本发明对上述方法进行了改进，如二、四次方谱谱线的提取方法以及MFSK信号的谱峰特征的提取等，同时结合自己所提出的特征参数(如圆环特征参数等)构造了一个新的调制识别判别体系。能够快速、准确地对通信信号的调制方式进行识别。

附图说明

图1为本发明所述的一种基于特征提取的通信信号调制识别方法的流程图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

所涉及到的英文缩写说明如下：

AM 幅度调制

FM 频率调制

CW 单载波

ASK 振幅键控

MFSK M进制频移键控

BPSK 二相相移键控

QPSK 四相相移键控

π/4-QPSK π/4相移QPSK

OQPSK 偏移四相相移键控

MQAM M进制正交振幅调制

一种基于特征提取的信号调制方式识别方法，包含以下顺序的步骤：

1)首先对接收到的信号进行中值滤波，利用中值滤波所得信号进行带宽估计，并根据所得带宽对信号进行线性滤波处理；

2)对于接收到的信号，计算其归一化功率谱以及归一化功率谱中大于常数B时谱峰的个数K，其中B的范围为(-6,-4)；具体包含以下步骤：

a、令 $\left\{\begin{array}{c}{e}_0^f\left(n\right)=x\left(n\right)\\ e_0^b\left(n\right)=x\left(n\right)\end{array}\right.,$ 利用下式计算：

${\hat{k}}_m=\frac{-2\underset{n=m}{\overset{N-1}{\Σ}}{e}_{m-1}^f\left(n\right)e_{m-1}^{b*}(n-1)}{\underset{n=m}{\overset{N-1}{\Σ}}|e_{m-1}^f\left(n\right){|}^2+\underset{n=m}{\overset{N-1}{\Σ}}|e_{m-1}^b(n-1){|}^2},m=1,2,...,p$

这里p我们取10，x(n)＝I(n)+iQ(n)；

b、利用

${\hat{r}}_x\left(0\right)=\frac{1}{N}\underset{n=0}{\overset{N-1}{\Σ}}|x\left(n\right){|}^2$

得到m＝1时的参数：

${\hat{a}}_1\left(1\right)={\hat{k}}_1,{\mathrm{\ρ}}_1=(1-|{\hat{k}}_1{|}^2){\hat{r}}_x\left(0\right)$

c、由和下式求出和再依据式(1)估计其中

$\left\{\begin{array}{c}{e}_m^f\left(n\right)=e_{m-1}^f\left(n\right)+k_m{e}_{m-1}^b(n-1)\\ e_m^b\left(n\right)=e_{m-1}^b\left(n\right)+k_m^{*}{e}_{m-1}^f\left(n\right),m=1,2,...,p\\ e_0^f\left(n\right)=e_0^b\left(n\right)=x\left(n\right)\end{array}\right.$

d、依据下式的递推关系，求出m＝2时的及

$\left\{\begin{array}{c}{\hat{a}}_m\left(k\right)={\hat{a}}_{m-1}\left(k\right)+{\hat{k}}_m{\hat{a}}_{m-1}^{*}(m-k),k=1,2,...,m-1\\ {\hat{a}}_m\left(m\right)={\hat{k}}_m\\ {\widehat{\mathrm{\ρ}}}_m=(1-|{\hat{k}}_m{|}^2){\widehat{\mathrm{\ρ}}}_{m-1}\end{array}\right.$

e、重复上面的过程，直到m＝p，从而求得

${\hat{a}}_p\left(k\right),k=1,2,...,p$

f、将上式所求得的代入到下式当中：

$P\left(\mathrm{\ω}\right)=\frac{1}{|1+\underset{k=1}{\overset{p}{\Σ}}{a}_k\left(cos\right(\mathrm{\ω}k)-i\·sin\left(\mathrm{\ω}k\right)){|}^2}$

上式a_k分别对应着前面我们所得到的

g、对上式的ω在[-π,π)取512个等间隔点，注意我们的取值是以下的离散值：

h、对所得到的每一个值P(ω)先归一化然后再进行dB转换，即

PW(ω)＝10log₁₀|P(ω)/max(P(ω))|

然后在这些点当中统计谱峰个数；此时所得的功率谱是利用自回归模型估计信号的平滑功率谱，即AR功率谱；

3)当K为1时，则判断该信号的调制方式为CW、ASK、AM中的一种：首先计算其归一化零中心瞬时幅度之谱密度的最大值γ_max，若γ_max小于2则为CW，若γ_max大于2则为ASK或AM；然后计算对数幅度谱的峰度系数β，若β小于4则为ASK，若β大于4则为AM；

所述的γ_max计算式如下：

${\mathrm{\γ}}_{max}=\frac{max|fft\[a_{cn}\left(i\right)\]|}{N_s}$

其中N_s为每一帧数据的个数，a_cn为零中心归一化瞬时幅度；

a_cn由下式计算得到：

a_cn(i)＝a_n(i)-1

其中 $a_n\left(i\right)=\frac{a\left(i\right)}{m_a},m_a=\frac{1}{N_s}\underset{i=1}{\overset{N_s}{\Σ}}a\left(i\right);$

所述的β由下式计算：

$\mathrm{\β}=\frac{\frac{1}{n}\underset{i}{\Σ}{(v_i-\stackrel{\‾}{v})}^4}{{\left(\frac{1}{n}\underset{i}{\Σ}{(v_i-\stackrel{\‾}{v})}^2\right)}^2},$

其中v_i代表信号频谱第i个点的电平值，单位为dB；

4)当K不为1且γ_max小于2、圆环特征参数num小于8时，则该信号的调制方式为MFSK或FM：计算其凹凸参数P，若P小于门限t₁则为MFSK,若P大于门限t₁则为FM，其中t₁的范围为(6，9)；所述的圆环特征参数num为统计归一化IQ平面内以原点为中心，0.1为半径的点的个数；

所述的P由下式计算：

$P=\frac{P_{21}+P_{22}}{P_1},$

假设在区域[f₁,f₂]中，f₁、f₂表示频率值，P₁是指信号载波附近区域中的电平最大值，P₂₁和P₂₂分别表示f₁右边区域的最大值和f₂左边区域的最大值；

5)当K不为1且不满足γ_max小于2、num小于8，则该信号的调制方式被判断为BPSK、QPSK、OQPSK、Π/4-QPSK、MQAM中的一种：首先计算其归一化平方功率谱，若大于门限t₂时只有一条谱线则为BPSK，若大于门限t₂时有两条谱线则为OQPSK；然后计算其归一化四次方功率谱，若大于门限t₃时只有一条谱线则为Π/4-QPSK；接着计算其高阶累计值F，若F小于16.5则为QPSK，否则为MQAM，其中t2、t3的范围为(0.3，0.45)。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (4)

1.一种基于特征提取的信号调制方式识别方法，其特征在于：包含以下顺序的步骤：

1)对于接收到的信号，计算其归一化功率谱以及归一化功率谱中大于常数B时谱峰的个数K，其中B的范围为(-6,-4)；

2)当K为1时，则判断该信号的调制方式为CW、ASK、AM中的一种：首先计算其归一化零中心瞬时幅度之谱密度的最大值γ_max，若γ_max小于2则为CW，若γ_max大于2则为ASK或AM；然后计算对数幅度谱的峰度系数β，若β小于4则为ASK，若β大于4则为AM；

3)当K不为1且γ_max小于2同时圆环特征参数num小于8时，则该信号的调制方式为MFSK或FM：计算其凹凸参数P，若P小于门限t₁则为MFSK,若P大于门限t₁则为FM，其中t₁的范围为(6，9)；

4)当K不为1且不满足γ_max小于2或num小于8，则该信号的调制方式被判断为BPSK、QPSK、OQPSK、Π/4-QPSK、MQAM中的一种：首先计算其归一化平方功率谱，若大于门限t₂时只有一条谱线则为BPSK，若大于门限t₂时有两条谱线则为OQPSK；然后计算其归一化四次方功率谱，若大于门限t₃时只有一条谱线则为Π/4-QPSK；接着计算其高阶累计值F，若F小于16.5则为QPSK，否则为MQAM，其中t2、t3的范围为(0.3，0.45)。

2.根据权利要求1所述的基于特征提取的信号调制方式识别方法，其特征在于：所述的步骤1)，具体包含以下步骤：

a、令 $\left\{\begin{array}{c}{e}_0^f\left(n\right)=x\left(n\right)\\ e_0^b\left(n\right)=x\left(n\right)\end{array}\right.,$ 利用下式计算：

${\hat{k}}_m=\frac{-2\underset{n=m}{\overset{N-1}{\mathrm{\Σ}}}{e}_{m-1}^f\left(n\right)e_{m-1}^{b^{*}}(n-1)}{\underset{n=m}{\overset{N-1}{\mathrm{\Σ}}}|e_{m-1}^f\left(n\right){|}^2+\underset{n=m}{\overset{N-1}{\mathrm{\Σ}}}|e_{m-1}^b(n-1){|}^2},m=1,2,...,p$

这里p我们取10，x(n)＝I(n)+iQ(n)，其中I(n)、Q(n)表示输入的I、Q数据序列；

b、利用

${\hat{r}}_x\left(0\right)=\frac{1}{N}\underset{n=0}{\overset{N-1}{\Σ}}|x\left(n\right){|}^2$

得到m＝1时的参数：

${\hat{a}}_1\left(1\right)={\hat{k}}_1,{\mathrm{\ρ}}_1=(1-|{\hat{k}}_1{|}^2){\hat{r}}_x\left(0\right)$

c、由和下式求出和再依据式(1)估计其中

$\left\{\begin{array}{c}{e}_m^f\left(n\right)=e_{m-1}^f\left(n\right)+k_m{e}_{m-1}^b(n-1)\\ e_m^b\left(n\right)=e_{m-1}^b(n-1)+k_m^{*}{e}_{m-1}^f\left(n\right),m=1,2,...,p\\ e_0^f\left(n\right)=e_0^b\left(n\right)=x\left(n\right)\end{array}\right.$

d、依据下式的递推关系，求出m＝2时的及

$\left\{\begin{array}{c}{\hat{a}}_m\left(k\right)={\hat{a}}_{m-1}\left(k\right)+{\hat{k}}_m{\hat{a}}_{m-1}^{*}(m-k),k=1,2,...,m-1\\ {\hat{a}}_m\left(m\right)={\hat{k}}_m\\ {\widehat{\mathrm{\ρ}}}_m=(1-|{\hat{k}}_m{|}^2){\widehat{\mathrm{\ρ}}}_{m-1}\end{array}\right.$

e、重复上面的过程，直到m＝p，从而求得

${\hat{a}}_p\left(k\right),k=1,2,...,p$

f、将上式所求得的代入到下式当中：

$P\left(\mathrm{\ω}\right)=\frac{1}{|1+\underset{k=1}{\overset{p}{\Σ}}{a}_k\left(cos\right(\mathrm{\ω}k)-i\·sin(\mathrm{\ω}k\left)\right){|}^2}$

上式a_k分别对应着前面我们所得到的

g、对上式的ω在[-π,π)取512个等间隔点，注意我们的取值是以下的离散值： $-\mathrm{\π},-\mathrm{\π}+\frac{2\mathrm{\π}}{512},-\mathrm{\π}+2\·\frac{2\mathrm{\π}}{512},-\mathrm{\π}+3\·\frac{2\mathrm{\π}}{512},...,\mathrm{\π}-\frac{2\mathrm{\π}}{512};$

h、对所得到的每一个值P(ω)先归一化然后再进行dB转换，即

PW(ω)＝10log₁₀|P(ω)/max(P(ω))|

然后在这些点当中统计谱峰个数；此时所得的功率谱是利用自回归模型估计信号的平滑功率谱，即AR功率谱。

3.根据权利要求1所述的基于特征提取的信号调制方式识别方法，其特征在于：步骤3)中，所述的圆环特征参数num为统计归一化IQ平面内以原点为中心，0.1为半径的点的个数。

4.根据权利要求1所述的基于特征提取的信号调制方式识别方法，其特征在于：在步骤1)之前还包括对接收到的信号进行预处理，具体如下：首先对接收到的信号进行中值滤波，利用中值滤波所得信号进行带宽估计，并根据所得带宽对信号进行线性滤波处理。