数字调制信号码速率估计方法和装置与流程

本发明涉及一种数字调制信号码速率估计方法和装置，属于通信技术领域。

背景技术：

通信信号调制识别目的就是对截获的一段通信信号，依据较少的先验信息，判断出信号的调制方式。随着通信技术的飞速发展，通信信号的体制和调制样式变得更加复杂多样，信号环境日趋密集，这使得通信信号的识别变得越加困难。通信信号中包含两大类调制信号：数字调制信号和模拟调制信号，利用判决理论时，模拟调制信号和数字调制信号的区分门限往往难以设置。由于模拟调制信号没有码速率之说，其估计结果为任意无规律值，因此利用码速率估计的方法可以有效区分模拟和数字调制，码速率是区分两种调制方式的重要特征。

码速率对于数字调制信号特征提取和分析以及识别和解调有着重要的意义，在电子侦察中，数字调制信号通过正交变频和低通滤波后得到系带信号，为了分析与识别信号，需要知道数字调制信号的码速率。因此正确估计数字调制信号的码速率是十分必要的。

估计数字调制信号码速率的方法有多种方法，其核心思想是根据调制信号码元之间存在突变，根据突变之间的间隔来估计码速率。对于数字调制基带信号，可以运用延迟相乘的方法来获得信号的码速率。经过延迟相乘之后，输出信号的双边功率谱密度主要包含三项，第一项是直流分量，第二项是码速率及高次谐波，第三项为连续谱。可以根据第二项对调制信号的码速率进行估计。但是时延长度必须大于码元宽度的一半，这在实际的应用中很难保证。

论文《无线电信号的调制识别研究》(作者王生兵)公开了利用承载数字信息的幅度调制波形p(n)进行码速率估计。对于ask信号，p(n)就是信号的瞬时幅度；对于psk信号，p(n)就是信号的瞬时相位；对于fsk信号，p(n)就是信号的瞬时频率。该方法的流程可概括如下：

1.对p(n)求l阶差分序列，并求取绝对值，得到p′(n)序列，l选为15；

2.计算p′(n)序列的均值mp，以3mp为门限对p′(n)序列进行判决；

3.当p′(n)＞3mp，判决器输出1；否则，判决器输出0；得到阶跃位置序列l(n)

4.以阶跃带的中心作为符号转换点，得到符号转换位置序列lc(n)，在对lc(n)取一阶差分，得到差分序列lcd(n)，其长度为nd；

5.令nm＝min[lcd(n)]，计算序列ni(n)＝round(lcd(n)/nm)，得到平均符号间距的估计

6.采集信号的码速率估计为：其中fb为输入数据的速率。

该方法的不足之处在于码速率估计方法中的门限通过经验取值，缺乏理论依据，门限的取值对不同参数的调制信号不具有适应性；

除此之外，因小波变换对瞬变信号有着很好的检测性能，基于小波变换的码速率估计是另外一种思路。文献【基于vssr法的数字与模拟调制信号的区分-戚玉鹏】对该方法进行了详细的描述，其码速率估计方法可概括为如下：

1.对采集的信号s(t)进行小波变换，得到a(t)；

2.对a(t)小波变换后取包络值，得到b(t)；

3.对b(t)求其傅里叶变换，得到c(f)；

4.求得|c(f)|中的相邻两尖峰之间的距离d，该距离d与|c(f)|的点个数的比值，即为符号率与信号的采样率的比值，从而可求得信号的符号率。

基于小波变换的方法，对于小波变换的尺度选择没有给出说明，如果对采集信号的码速率完全未知，很难选择合适的小波变换的尺度，因而导致码速率估计错误；

此外现有码速率估计方法对各类数字调制信号分析的不够全面，一般是对某种特定调制方式的信号进行码速率估计，缺乏一种统一的码速率估计方法且可能依赖于先验知识，如需要预先知道调制方式(幅度、相位和频率)、信号的载频和载波相位等，实际情况很难准确获知这些参数；

技术实现要素：

本发明的目的是为了克服目前码速率估计技术的不足，提出一种无需先验知识的基于小波变换的数字调制信号码速率估计方法，其自适应程度高。

为解决上述技术问题，本发明采用了以下的技术方案：

数字调制信号码速率估计方法，包括：

步骤1：对以预设采样率fs及总采样点数l进行采样后的信号提取瞬时幅度a(l)和瞬时频率f(l)；

步骤2：按照设定的小波变换尺度分别对瞬时幅度a(l)和瞬时频率f(l)进行小波变换，小波变换母函数选用haar小波，得到ca(l)和cf(l)；

步骤3：对ca(l)进行离散傅里叶变换，得到ya(fk)；对频率f(l)进行离散傅里叶变换，得到yf(fk)；

步骤4：提取ya(fk)和yf(fk)的脉冲尖峰频率位置并对所得位置进行排序，分别得到序列xa(n1)和序列xf(n2)，n1和n2为正整数；

步骤5：对序列xa(n1)做差分，得到dxa(n1)序列，对dxa(n1)中的数值进行统计：选出dxa(n1)中重复次数最多的数值，记为xai,相应的重复次数为nai；对序列xf(n2)做差分，得到dxf(n2)序列，对dxf(n2)中的数值进行统计：选出dxf(n2)中重复次数最多的数值，记为xfj,相应的重复次数为nfj；其中下角标i和j都是整数且初始值为零；

步骤6：按照预设条件改变小波变换尺度因子m，并令i等于i加1，j等于j加1，并重复步骤2～5，直到小波变换的变换尺度a大于数据长度，并记录此时的i和j均等于m；

步骤7：挑选出na0,na1,na2,…nam的最大值，并记录其下标为ma，ma为0到m中的某一整数，把xama记为rsa，nama记为票数pa；挑选出nf0,nf1,nf2,…nfm的最大值，并记录其下标为mf，mf为0到m中的某一整数，把xfmf记为rsf，nfmf记为票数pf；

步骤8：若pa大于pf，则码速率rs为rsa与采样率fs的乘积与总采样点数l的比值；否则，则码速率rs为rsf与采样率fs的乘积与总采样点数l的比值。

进一步地，提取瞬时幅度a(l)和瞬时频率f(l)包括：

对于连续波调制，已调实信号的表达式为：

式中，fc表示载波的频率，a(t)表示已调信号的瞬时幅度；f(t)表示频率；为相位；

对信号进行希尔伯特变换，得到信号的正交分量，表达式如下：

得到信号s(t)的解析表达式z(t)＝s(t)+jv(t)。

再进一步地，提取提取瞬时幅度a(l)，方法如下：

对信号进行采样后，得到离散形式的瞬时幅度a(l)，l＝1,2,…l。

进一步地，提取瞬时频率f(l)的方法：瞬时相位提取，瞬时相位表达式如下：

对信号采样后，得到离散形式的瞬时相位序列，l＝1,2,…l；

去相位折叠，首先计算修正相位序列c(l),l＝1,2,…l：

则无折叠的相位为：

去载频线性相位：采用相性规划法估计出线性相位分量并使误差最小，求出c1和c2两个常数，表示为：

其中为非线性相位分量，为无折叠的相位；

瞬时频率提取，表达式为

其中fs为采样率。

进一步地，小波变换尺度因子为m，m取不小于零的整数，则小波变换的变换尺度a的表达式为：a＝2^m。

再进一步地，设定初始小波变换尺度因子m为零。

进一步地，改变小波变换尺度因子时令小波变换尺度因子m等于m加1。

进一步地，脉冲尖峰频率位置的提取方法如下：

|y(fk)|在码速率对应的频率及码速率倍频位置会出现脉冲尖峰，如果满足以下表达式：

|y(fk)|＞|y(fk-1)|&|y(fk)|＞|y(fk+1)|(9)，

则表明在频率fk处存在脉冲尖峰。

在另一方面本发明提供了数字调制信号码速率估计装置，包括：

采样单元：用于根据预设采样率fs及总采样点数l对信号进行采样；

瞬时幅度和瞬时频率提取单元：用于对采用后的信号进行瞬时幅度a(l)和瞬时频率f(l)的提取；

小波变换单元：用于按照设定的小波变换尺度分别对瞬时幅度a(l)和瞬时频率f(l)进行小波变换，小波变换母函数选用haar小波，得到ca(l)和cf(l)；

离散傅里叶变换单元：用于对ca(l)进行离散傅里叶变换，得到ya(fk)；对频率f(l)进行离散傅里叶变换，得到yf(fk)；脉冲尖峰频率位置提取单元：用于提取ya(fk)和yf(fk)的脉冲尖峰频率位置并对所得位置进行排序，分别得到序列xa(n1)和序列xf(n1)，n1和n2为正整数；

差分单元：用于对序列xa(n1)做差分，得到dxa(n1)序列以及对序列xf(n2)做差分，得到dxf(n2)序列；

票数统计单元：用于对dxa(n1)中的数值进行统计：选出dxa(n1)中重复次数最多的数值，记为xam,相应的重复次数为nam；并且对dxf(n2)中的数值进行统计：选出dxf(n2)中重复次数最多的数值，记为xfm,相应的重复次数为nfm；

其中下角标m是每次循环都不同的变量；

小波变换尺度因子设定单元：用于设定的小波变换尺度并且按照预设条件改变小波变换尺度因子直到小波变换的变换尺度a大于数据长度，并记录m等于此时的小波变换尺度因子m，m取不小于零的整数；

码速率计算单元：用于挑选出na0,na1,na2,…nam的最大值，并记录其下标为ma，ma为0到m中的某一整数，把xama记为rsa，nama记为票数pa；

挑选出nf0,nf1,nf2,…nfm的最大值，并记录其下标为mf，mf为0到m中的某一整数，把xfmf记为rsf，nfmf记为票数pf；

若pa大于pf，则码速率rs为rsa与采样率fs的乘积与总采样点数l的比值；否则，则码速率rs为rsf与采样率fs的乘积与总采样点数l的比值。

本发明所达到的有益效果：

1.数字调制信号的种类多，对常见信号都进行了全面分析和码速率估计；

2.不依赖于先验知识，包括信号调制方式、信噪比、载频、载波相位和码元同步等，在对信号完全未知的情况下进行码速率估计，并且在较低信噪比的情况下识别正确率高；

3.本发明提出的“票数统计”策略能够有效解决小波变换尺度选择的问题；

4.具有自适应性，对于数字调制信号的不同码速率都能够自动选择合适的尺度；算法规则简单且收敛，码速率的估计可以在很短时间内完成，适合在线分析。

附图说明

图1是本发明方法流程示意图。

具体实施方式

下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

图1为本发明方法流程示意图，图1示出了本发明数字调制信号码速率估计方法，

步骤1：对以预设采样率fs及总采样点数l进行采样后的信号提取瞬时幅度a(l)和瞬时频率f(l)，l为采样后的离散自变量；因为信号未知的情况下，无法确定码元调制在幅度上还是频率上，因此本发明方法分幅度和频率处理：

具体方法包括：

对于连续波调制，已调实信号的一般表达式为：

式中，fc表示载波的频率。调制信号可以分别“寄生”在已调信号的瞬时幅度a(t)、频率f(t)和相位中。

(1).对信号进行希尔伯特变换，得到信号的正交分量。

希尔伯特变换实际是对原信号进行了π/2相移，使它成为自身的正交对。可得信号s(t)的解析表达式z(t)＝s(t)+jv(t)。

(2).瞬时幅度提取

信号包络瞬时幅度估计：

(3).瞬时相位提取

瞬时相位估计：

(4).对信号采样后，得到离散形式的瞬时幅度a(l)序列和瞬时相位序列，l＝1,2,…l

(5).去相位折叠

由于实际相位计算时是以模2π来计算的，因此在求得瞬时相位之后需要对相位去折叠，恢复出没有折叠的相位，为此首先计算修正相位序列c(l),l＝1,2,…l：

则无折叠的相位为：

(6).去载频线性相位

信号的相位包括两部分，载频引起的线性相位分量和调制相位的非线性相位分量因此非线性相位分量可由下式计算：

在载频fc不能精确已知的情况下，可以采用相性规划法估计出线性相位分量并使误差最小，则可以求出c1和c2两个常数，所以有：

(7).瞬时频率提取

瞬时频率可以由瞬时相位的导数计算求得，对数字信号可直接求差分，即：

步骤2：按照设定的小波变换尺度分别对瞬时幅度a(l)和瞬时频率f(l)进行小波变换，由于haar小波简单、计算快和最常用，因此小波变换母函数选用haar小波，得到ca(l)和cf(l)；

小波变换尺度因子为m，m取不小于零的整数，则小波变换的变换尺度a的表达式为：a＝2^m。

从图1可以看出设定初始小波变换尺度时设小波变换尺度因子为m为零。

步骤3：对ca(l)进行离散傅里叶变换，得到ya(fk)；对频率f(l)进行离散傅里叶变换，得到yf(fk)；

步骤4：提取ya(fk)和yf(fk)的脉冲尖峰频率位置并对所得位置进行排序，分别得到序列xa(n1)和序列xf(n2)，n1和n2为正整数；

ya(fk)和yf(fk)在码速率对应的频率及码速率倍频位置会出现脉冲尖峰，记录尖峰对应频率位置；在对所得位置进行排序时升序和降序均可。

步骤5：对序列xa(n1)做差分，得到dxa(n1)序列，对dxa(n1)中的数值进行统计：选出dxa(n1)中重复次数最多的数值，记为xai,相应的重复次数为nai；对序列xf(n2)做差分，得到dxf(n2)序列，对dxf(n2)中的数值进行统计：选出dxf(n2)中重复次数最多的数值，记为xfj,相应的重复次数为nfj；其中下角标i和j都是整数且初始值为零；

步骤6：按照预设条件改变小波变换尺度因子m，并令i等于i加1，j等于j加1，并重复步骤2～5，直到小波变换的变换尺度a大于数据长度，并记录此时的i和j均等于m；

优选地，改变小波变换尺度因子时令小波变换尺度因子m等于m加1。

步骤7：挑选出na0,na1,na2,…nam的最大值，并记录其下标为ma，ma为0到m中的某一整数，把xama记为rsa，nama记为票数pa；挑选出nf0,nf1,nf2,…nfm的最大值，并记录其下标为mf，mf为0到m中的某一整数，把xfmf记为rsf，nfmf记为票数pf；

步骤8：若pa大于pf，则码速率rs为rsa与采样率fs的乘积与总采样点数l的比值；否则，则码速率rs为rsf与采样率fs的乘积与总采样点数l的比值。

小波变换原理

定义：ψ(t)∈l²(r)，其傅里叶变换为当满足允许条件(完全重构条件或恒等分辨条件)：

称ψ(t)为一个基本小波或者母小波。将母函数经过平移和伸缩后得：

称其为一个小波序列。其中a为伸缩因子，b为平移因子。对于任意的函数x(t)∈l²(r)的连续小波变换为：

对应采样后信号的小波变换为：

本发明中，x(l)对应于瞬时幅度a(l)和瞬时频率f(l)。

由于基小波在小波变换中对被分析的信号起着观测窗的作用，所以ψ(t)还应该满足一般函数的约束条件：

通常，把连续小波变换中尺度参数a和平移参数b分别取做离散形式：j∈z称为离散小波变换。

(2).haar母小波

haar母小波的数学表达式为

脉冲尖峰频率位置的提取方法

|y(fk)|在码速率对应的频率及码速率倍频位置会出现脉冲尖峰，如果满足以下规则：

|y(fk)|＞|y(fk-1)|&|y(fk)|＞|y(fk+1)|(17)

则表明在频率fk处存在脉冲尖峰，‘&’表示“与门”。为了消除噪声起伏带来的影响，在算法实现中，仅考虑3db带宽对应频率内的脉冲尖峰。

在另一方面本发明提供数字调制信号码速率估计装置，包括：

采样单元：用于根据预设采样率fs及总采样点数l对信号进行采样；

瞬时幅度和瞬时频率提取单元：用于对采用后的信号进行瞬时幅度a(l)和瞬时频率f(l)的提取；

小波变换单元：用于按照设定的小波变换尺度分别对瞬时幅度a(l)和瞬时频率f(l)进行小波变换，小波变换母函数选用haar小波，得到ca(l)和cf(l)；

离散傅里叶变换单元：用于对ca(l)进行离散傅里叶变换，得到ya(fk)；对频率f(l)进行离散傅里叶变换，得到yf(fk)；脉冲尖峰频率位置提取单元：用于提取ya(fk)和yf(fk)的脉冲尖峰频率位置并对所得位置进行排序，分别得到序列xa(n1)和序列xf(n1)，n1和n2为正整数；

差分单元：用于对序列xa(n1)做差分，得到dxa(n1)序列以及对序列xf(n2)做差分，得到dxf(n2)序列；

票数统计单元：用于对dxa(n1)中的数值进行统计：选出dxa(n1)中重复次数最多的数值，记为xam,相应的重复次数为nam；并且对dxf(n2)中的数值进行统计：选出dxf(n2)中重复次数最多的数值，记为xfm,相应的重复次数为nfm；

其中下角标m是每次循环都不同的变量；

小波变换尺度因子设定单元：用于设定的小波变换尺度并且按照预设条件改变小波变换尺度因子直到小波变换的变换尺度a大于数据长度，并记录m等于此时的小波变换尺度因子m，m取不小于零的整数；

码速率计算单元：用于挑选出na0,na1,na2,…nam的最大值，并记录其下标为ma，ma为0到m中的某一整数，把xama记为rsa，nama记为票数pa；

挑选出nf0,nf1,nf2,…nfm的最大值，并记录其下标为mf，mf为0到m中的某一整数，把xfmf记为rsf，nfmf记为票数pf；

若pa大于pf，则码速率rs为rsa与采样率fs的乘积与总采样点数l的比值；否则，则码速率rs为rsf与采样率fs的乘积与总采样点数l的比值。

进一步地，所述小波变换单元设定初始小波变换尺度因子m为零，小波变换的变换尺度a的表达式为：a＝2^m。

在具体实施例中实际通信系统中常见的数字调制信号2ask,4ask,bpsk,qpsk,8psk,oqpsk,pi/4-dqpsk,2fsk,4fsk,16qam,64qam为例，进行码速率的自动估计计算机仿真，使用随机序列作为数字调制信号的信源。用均值为0的正态分布的随机数发生器产生加性高斯白噪声awgn序列。信噪比snr变化范围从5db到15db，步进为1db。每个识别样本用20000个点，每个识别样本仿真400次。仿真平台操作系统是win7，采用软件matlab进行编程。仿真表明，当信噪比达到6db时，码速率的估计准确率高达99％。

现代通信系统中最常用的是数字调制技术，因此，数字调制方式的识别以及关键参数码速率的估计对通信信号的分析以及特征提取有着至关重要的作用。本发明提出的基于小波变换的数字调制信号码速率估计方法无需任何先验知识，并且根据“投票”策略能够自适应确定小波变换的尺度，算法规则简单且收敛。对常见数字调制信号的仿真结果表明，在信噪比较低的情况下，码速率的估计正确率仍能达到99％以上

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。