一种多进制频移键控信号的检测和分析方法

专利名称：一种多进制频移键控信号的检测和分析方法
技术领域：
本发明属于通信信号分析领域，特别涉及数字调制方式中多进制频移键控(MFSK)信号的自动识别和分析。
背景技术：
调制方式自动识别是电子战(通信侦察)、无线电监测、自适应通信等应用的重要环节，是非协作通信领域的一个重要课题，通信信号分析的各种理论在这一方向获得了广泛应用。频移键控调制方式(FSK)在军队、政府的短波/超短波通信和业余无线电通信中扮演着重要角色，是一种重要的数字信号调制方式，对该类信号的检测和分析是无线电信号监测和军事通信侦察的重要内容。本发明主要涉及频移键控信号的检测和分析，即频移键控信号的识别和调制阶数的分析。
近年来，针对多进制频移键控(MFSK)信号检测和分析，已经提出了如下几种方法1.基于信号包络分析和瞬时频率分析的方法E.E.Azzouz等人在Automatic modulation recognition(Journal Of The FranklinInstitute-Engineering And Applied Mathematics，Mar 1997)中提出了使用信号包络的的谱峰值γmax对恒包络(频率调制信号)信号和非恒包络信号进行分类，使用归一化瞬时频率绝对值的标准偏差σaf对BFSK信号和QFSK信号分类，步骤如下1).将接收到的信号转化为解析信号，去掉弱信号段；2).计算零中心归一化瞬时幅度acn；3).计算信号包络的谱峰值γmax，将该值同预定义门限τγ相比较，若γmax＜τγ则将接收信号判定为频率调制信号；
4).通过解析信号的瞬时相位计算其零中心归一化瞬时频率fcn；5).计算归一化瞬时频率绝对值的标准偏差σaf，将该值同预定义门限τσaf相比较，若σαf＜τσ则判定为BFSK信号，否则判定为QPSK信号。
σaf=1NsΣfcn2(n)-[1NsΣfcn]2]]>该方法的问题在于●由于频率选择性衰落等因素影响，许多实际接收的MFSK信号样本并非恒包络信号；●该方法只能识别BFSK和QFSK信号。
曹志刚等人在相关专利“无需先验知识的常用数字信号调制方式的自动识别方法”(专利号02123627.5)中也涉及到了MFSK信号的检测和分析。相关的步骤如下1).计算信号样本的R参数，若该参数小于预定义门限，则将该信号归为恒包络信号；2).对恒包络信号，计算信号s(n)＝x(n)+iy(n)的瞬时相位 3).对瞬时相位进行去卷叠，给模2π的相位序列加上相位校正序列θ(n)＝(n)+Ck(n) 4).计算信号的瞬时频率f1(n)f1(n)=f11(n)-1NΣn=1Nf11(n),f11(n)=θ(n)-θ(n-1);]]>5).对f1(n)进行绝对值大于π/2的统计，得到统计次数Cp，满足Cp＜TCp的信号判定为没有相位跳变的信号，即MFSK信号；6).对f1(n)进行插值处理，得到f2(n)f2(n)=f1(n),|f1(n)|<π212[f1(n+1)+f1(n-1)],|f1(n)|≥π2;]]>7).对f2(n)进行归一化处理得到归一化频率f(n)f(n)=f2(n)max{f2(n)};]]>8).对f(n)的分布进行统计，得到f(n)的峰值个数Nf，若Nf＝2，则判定为BFSK信号，否则判定为FM信号。
该方法的问题在于●由于频率选择性衰落等因素影响，许多实际接收的MFSK信号样本并非恒包络信号；●该专利没有给出计算f(n)的统计分布的峰值个数Nf的方法，并且只能识别BFSK信号。
2.基于过零采样的方法S.-Z.Hsue等人在论文Automatic Modulation Classification Using Zero Crossing(Radar andSignal Processing，IEE Proceedings For，vol.37，No.6，pp.459-464，Dec.1990.K.)中提出了使用过零采样器提取信号的瞬时频率的倒数，然后通过分析该倒数的方差来检测MFSK信号，通过分析瞬时频率分布直方图来分析MFSK信号的方法，步骤如下1).接收信号，用过零采样器记录信号的过零的时间点序列x(i)；2).通过信号过零的时间x(i)计算值序列y(i)＝x(i+1)-x(i)和值序列z(i)＝y(i+1)-y(i)；3).设定门限τz，将满足z(i)＞τz的样本点定义为“符号间跳变(IST)样本”，去掉相应的y(i)，得到修正后的值序列ya(i)；4).计算ya(i)的方差G，并且设定门限T，满足G＞T的信号样本被判定为MFSK信号；5).将IST样本点之间的ya(i)取均值，计算其分布直方图；6).通过分析该直方图的峰值个数计算MFSK信号的调制阶数。
该方法的问题在于●过零频率估计的准确性对SNR非常敏感，在相关报道中，作者只给出了CNR≥15dB的仿真结果；●过零法测量瞬时频率需要很高的采样率；●在低信噪比条件下，通过分析过零频率估计所得到的频率分布柱状图的峰值个数计算准确的调制阶数比较困难。
3.基于信号离散傅立叶变换的方法Zaihe Yu等人在论文A practical classification algorithm for M-ary frequency shift keyingsignals(Military Communications Conference，2004.MILCOM 2004.IEEE Volume 2，31 Oct.-3Nov.2004 Page(s)1123-1128 Vol.2)中提出了一种通过直接分析MFSK信号功率谱来得到MFSK信号调制阶数的方法。该报道中没有MFSK信号的检测方法，而是假设待分析信号是MFSK信号。该方法步骤如下1).计算信号样本的离散傅立叶变换(DFT)；2).按照门限规则、局部极值规则、等距离规则、全信息规则、分组规则、自适应门限规则和可信度规则分析信号样本DFT峰值个数，得到调制阶数M。
该方法的问题在于●该方法假设各个码元在功率上的分布是相等的，但是，在实际的测试环境中，由于频率选择性衰落等因素，各个码元功率分布往往不是相等的；●在计算DFT的的时候，作者声明使用尽可能长的FFT长度，这不利于实际系统的实现，并且在信号带宽未知的情况下，估计一个合理的FFT长度也是不容易的。

发明内容
本发明的主要目的在于识别和分析多进制频移键控(MFSK)信号，包括一种MFSK信号的检测方法和一种MFSK信号的分析方法。
本发明的特征在于，所述方法是在计算机中依次按以下步骤实现的步骤(12)接收待处理数据，得到实信号x(n)；步骤(13)按以下步骤把实信号x(n)转换成解析信号s(n)步骤(13.1)把实信号x(n)补零后计算FFT，得到信号的离散傅立叶变换Fx(n)，补零的长度Nz由下式得到Nz＝NFFT-Nx，Nx为实信号x(n)的长度， 表示大于“·”的最小整数，NFFT是进行FFT运算的长度；步骤(13.2)取所述Fx(n)的前半部分，n＝1，……，NFFT/2，得到信号s(n)的离散傅立叶变换Fs(n)；步骤(13.3)把Fs(n)的前半部分，n＝1，……，NFFT/4，和后半部分n＝NFFT/4+1，……，NFFT/2调换；步骤(13.4)计算步骤(2.3)得到的Fs(n)的IFFT，得到解析信号s(n)，s(n)＝IFFT(Fs(n))；步骤(13.5)把步骤(2.4)得到的解析信号s(n)的 部分截去，剩余部分即为后续分析用的解析信号，其有效长度为 步骤(14)根据预定义门限，分段检测s(n)的幅值，去掉弱信号段，其步骤依次如下
步骤(14.1)把步骤(2.5)得到的信号s(n)＝x(n)+jy(n)分成等长的Nseg段，用一个序列Si，i＝1，2，……，Nseg表示，Nseg＝5～20；步骤(14.2)计算步骤(3.1)所述的每一段信号幅值的和mi
步骤(14.3)设定门限τmτm=12max{mi,i=1,2,······,Nseg};]]>步骤(14.4)将mi＜τm的信号段判定为弱信号段，去掉所有弱信号段，把剩余的信号段连接成新的解析信号序列s(n)，长度为N；步骤(15)计算步骤(3.4)所述解析信号s(n)的归一化瞬时频率序列fn(n)，按以下步骤实现步骤(15.1)按下式计算解析信号序列s(n)＝x(n)+jy(n)，n＝1，2，……，N的瞬时相位序列(n) (n)∈(-π，π)，n＝1，2，……，N；步骤(15.2)对步骤(4.1)得到的瞬时相位序列(n)做差分运算后取余得到信号的瞬时频率序列f(n)f(n)＝mod((n+1)-(n)+π，2π)，n＝1，2，……，N-1；步骤(15.3)把该瞬时频率序列f(n)归一化为零均值、单位方差序列fn(n)fn(n)=f(n)-mfσf,]]>其中，mf=1N-1Σn=1N-1f(n),σf=1N-2Σn=1N-1(f(n)-mf)2;]]>步骤(16)计算步骤(4.3)得到的fn(n)的归一化瞬时频率的变化率序列Δfn(n)及其绝对值的均值mΔfn，其中Δfn(n)＝|fn(n+2)-fn(n)|，n＝1，2，……，N-3，mΔfn=1N-3Σn=1N-3Δfn(n)]]>
步骤(17)把所述mΔfn同预定义门限τm比较若mΔfm<τm,]]>则把所述信号x(n)判定为MFSK信号，否则，便非MFSK信号；步骤(18)依次执行以下步骤以利用均值滤波器来减小瞬时频率的抖动步骤(18.1)利用均值滤波器对瞬时频率值序列进行平滑滤波，得到更新的瞬时频率序列fn(n)，滤波器长度lf，由被分析信号可能的最大码元速率fdmax和信号的设定采样频率fs按下式求出lf=mfsfdmax,m=13~12;]]>步骤(18.2)根据步骤(7.1)得到的更新的瞬时频率序列重新计算信号的归一化瞬时频率变化率序列Δfn(n)Δfn(n)＝|fn(n+2)-fn(n)|；步骤(18.3)根据下式计算步骤(7.2)得到的归一化瞬时频率变化率序列的方差τΔfτΔfn=1N-4Σn=1N-3(Δfn(n)-mΔfn)2,mΔfn=1N-3Σn=1N-3Δfn(n);]]>步骤(18.4)搜索Δfn(n)序列，若|Δfn(n)|<τΔfn,]]>则fn(n+1)属符号稳定区间Гsi，i＝1，2，……，Ns(Ns为符号稳定区间的个数)，若|Δfn(n)|>τΔfn]]>则fn(n+1)属符号跳变区间Гti，i＝1，2，……，Nt(Nt为符号跳变区间的个数)；步骤(18.5)把步骤(7.4)中每一个符号稳定区间的瞬时频率值改为该区间所有瞬时频率值的均值，舍弃所有符号跳变区间的瞬时频率值，构成新的瞬时频率值序列fn(n)；步骤(19)根据“熵原则”按以下步骤分析步骤(7.5)所述瞬时频率序列值的直方图步骤(19.1)根据归一化瞬时频率序列fn(n)建立分布直方图 以此作为fn(n)的概率密度分布函数的估计p^(fi)=count{fn(n),fi≤fn(n)<fi+Δf}Nfn,fi=min{fn(n)}+iΔf,i=0,1,······,NP^-1,]]>其中，Nfn为fn(n)的长度，Δf=(max{fn(n)}-min{fn(n)})/Np^,]]> 为分布直方图的长度，Np^=512~4096,]]>count(x，y}表示满足条件y的x值的个数步骤(19.2)确定门限τp的搜索范围为 搜索步长Δτp=max{p^(fi)}/2N,]]>
N＝10～100步骤(19.3)对搜索范围内的每个、值执行一下步骤步骤(19.3.1)搜索满足p^(fi)>τp]]>的fi，把连续的fi归为一组；步骤(19.3.2)计算每一组的概率密度函数估计 之和，得到MFSK信号各个码元的分布概率的估计值 步骤(19.3.3)计算对应该τp的熵值ee=-Σi=1M^P^klog(P^k),]]> 是连续的男值组成的组的个数；步骤(19.4)取对应熵值最大的τp，将相应的fi组的个数作为码元状态个数 的初始估计，每组的最小瞬时频率 和最大瞬时频率 作为每个码元状态的频率范围(f^minkf^maxk),k=1,2,······,M^;]]>步骤(20)根据“等概率分布原则”，按以下步骤精确估计码元状态个数 的精确值步骤(20.1)设定门限0<τPmin<1;]]>步骤(20.2)去掉P^k<τPminmax{P^k,k=1,2,······,M^]]>的状态；步骤(20.3)重新计算码元状态个数 的初始估计和每个码元状态的频率范围(f^mink,f^maxk),k=1,2,······,M^;]]>步骤(21)根据“等频差原则”精确估计码元状态个数 和各状态中心频率，其步骤如下步骤(21.1)计算各个码元状态的中心频率 f^mk=Σf=fminkfmaxkp^(f)fΣf=fminkfmaxkp^(f);]]>步骤(21.2)计算各个状态之间的频差 f^dkl=|f^mk-f^ml|;]]>步骤(21.3)计算每个频差的权重wkl权重等于得到该频差的两个码元状态的概率分布估计的较小值，减去这两个码元状态之间的所有码元状态的概率分布估计之和，如果差值小于0，则权重为0wkl=max{0,min{P^k,P^l}-Σn=k+1l-1P^n},k<l-1max{0,min{P^k,P^l}},k=l-1;]]>步骤(21.4)对各个频差加权平均，得到频差的估计f^d=Σk,lwklf^dklΣk,lwkl]]>步骤(21.5)设定门限τmin、τmaxτmin=0.4f^d~0.8f^d,τmax=1.2f^d~1.6f^d;]]>步骤(21.6)对于f^dk,k+1<τmin]]>的相邻状态，合并码元状态k和k+1；步骤(21.7)对于f^dk,k+1>τmax]]>的相邻状态，在码元状态k和k+1之间插入 个码元状态( 表示大于“·”的最小整数)；步骤(21.8)更新码元状态个数 和各个码元状态的中心频率 步骤(22)后续处理按下式计算得到MFSK信号的调制阶数 M^0=2[log2M^],]]>其中[·]运算表示取离“·”最近的整数。
本发明具有如下优点●不依赖于MFSK信号是恒包络信号的假设模型，允许信号出现频率选择性衰落，适用范围大于其他已有算法；●抓住了MFSK信号存在瞬时频率跳变的本质特征，以瞬时频率作为基本特征来识别和分析MFSK信号，具有较好的鲁棒性；●针对计算机模拟信号的测试和针对实际信号的测试表明，本发明提出的算法具有很好的识别效果和较强的实用性。
本发明达到的效果多进制频移键控(MFSK)是一类重要的数字信号调制方式，本发明主要涉及该类调制方式的通信信号的检测和分析。计算机仿真实验表明，对调制阶数为2～32的MFSK信号，信噪比大于5dB的条件下，本发明所提方法的检测正确率达到了99％以上，信噪比大于7.5dB的条件下，分析正确率达到了95％以上。对实际接收的信号样本的测试表明，本发明所提方法的检测正确率达到了95％以上，分析正确率达到了90％以上。


图1为本发明提出的MFSK信号检测和分析流程图。
图2为实施例中实际信号测试时采用的硬件原理图。
图3为仿真MFSK信号的识别测试结果，其中虚线表示FSK信号的mΔfn值，实线表示非FSK信号的mΔfn值。
图4为仿真MFSK信号的分析测试结果，每条曲线代表一种调制方式在不同信噪比条件下的分析正确率。
图5为实际接收MFSK信号的识别测试结果其中虚线表示FSK信号的mΔfn值，实线表示非FSK信号的mΔfn值。
具体实施例方式
本发明所提出的MFSK信号检测和分析方法总体流程如图1所示，包括以下步骤1)接收待处理数据，得到实信号x(n)；2)将实信号x(n)变成解析信号s(n)；3)根据预定义门限，分段检测s(n)幅值，去掉弱信号段；4)计算s(n)的归一化瞬时频率序列fn(n)；5)计算归一化瞬时频率的变化率序列Δfn(n)及其绝对值的均值mΔfn；6)将mΔfn同预定义门限τm比较，如果mΔfn<τm,]]>则将该信号判定为MFSK信号，否则判定该信号不是MFSK信号；7)利用均值滤波减小瞬时频率的抖动；8)根据“熵原则”分析瞬时频率分布直方图，得到码元状态个数 的初始估计；9)根据“等概率分布原则”精确估计码元状态个数 10)根据“等频差原则”精确估计码元状态个数 和各状态中心频率；11)后续处理。
这里给出针对计算机仿真信号的测试和针对实际接收信号的测试两个实施例。
1.针对计算机仿真信号的测试
a)MFSK信号的识别在不同的信噪比条件下，生成两组信号用于测试。一组是FSK信号，一组是非FSK信号●FSK信号组包含2FSK、4FSK、8FSK、16FSK、32FSK信号，信号长度256个码元，每个码元100个采样点，码元之间频差100Hz，采样率为最高调制频率的4倍；●非FSK信号组包含BPSK、QPSK、8PSK、16QAM、64QAM、OQPSK、CW信号，信号长度256个码元，采样频率为信号波特率的8倍，采用滚降系数0.8的升余弦窗脉冲整形。
实际应用环境中接收的信号，信噪比是未知的，也是不稳定的。为了更符合实际应用的情况，对生成的两组测试信号，随机加入了5dB至30dB的白噪声。最后得到两组仿真信号，各包含500段信号样本，分别计算每段信号的mΔfn值，估计合理的门限τm，得到如附图3的测试结果。附图3中横轴为信号样本编号，纵轴为mΔfn值。可以看出，FSK信号组和非FSK信号组被明显分离，合理的门限值范围是τm＝0.8～1.2。当τm＝0.84时，识别的正确率达到100％。
b)MFSK信号的分析在0～20dB信噪比条件下，生成2FSK、4FSK、8FSK、16FSK、32FSK信号用于测试。信号长度为1024个码元，每个码元100个采样点，码元之间频差200Hz，信号采样率为最高调制频率的4倍。用本发明提出的算法对生成的信号进行分析，对每种调制阶数/信噪比组合重复100次测试，得到了如附图4的测试结果。附图4是运用本算法对不同调制阶数的测试信号在不同的信噪比条件下进行分析得到的分析正确率的曲线，横轴是测试信号的信噪比，纵轴是分析正确率，每条曲线代表一种调制方式。可以看出，当SNR≥7.5dB时，对2FSK、4FSK、8FSK、16FSK、32FSK调制信号的分析正确率均达到了95％以上。
2.针对实际接收信号的测试(测试用硬件原理图见附图2)a)MFSK信号的识别同针对仿真信号的测试类似，将测试数据分成两组●MFSK信号，包含了2FSK(58)、4FSK(16)、6FSK(3)、8FSK(17)、12FSK(3)、13FSK(2)、16FSK(1)、20FSK(2)、32FSK(5)信号(括号内是该调制阶数信号的样本数)，采样率范围为6000～8000Hz，码元速率范围13Bd～1000Bd。
●非MFSK信号，包含了BPSK(8)、QPSK(2)、8PSK(9)信号(括号内是该调制方式信号的样本数)，采样率范围为6000～8000Hz，码元速率范围为30Bd～2400Bd。
从两组样本中，随机提取时间长度为1s(6000～8000采样点)的信号段，计算其mΔfn值，估计合理的门限τm，得到如附图5的测试结果。附图5中，横轴是样本信号段编号，纵轴是mΔfn值。从图中可以看出，两组样本基本上被很好的分离，合理的门限值范围是τm＝0.8～1.2，同仿真结果相符。取τm＝1.17的时候，分离的正确率达到了95.6％。
b)MFSK信号的分析我们利用实际接收信号对本发明提出的MFSK信号的分析算法进行了测试，由于足够长的16FSK信号和32FSK信号很难得到，测试用信号集只包括了2FSK、4FSK、8FSK信号。测试中将每段信号样本分成长度为1s(6000～8000个样本点)的数据段，用本发明提出的MFSK信号分析算法进行分析，得到的测试结果如下

从上表中的测试结果可以看出，该算法对实际信号的分析达到了比较高的正确率。
权利要求
1.一种多进制频移键控信号的检测和分析方法，其特征在于，所述方法是在计算机中依次按以下步骤实现的步骤(1)接收待处理数据，得到实信号x(n)；步骤(2)按以下步骤把实信号x(n)转换成解析信号s(n)步骤(2.1)把实信号x(n)补零后计算FFT，得到信号的离散傅立叶变换Fx(n)，补零的长度Nz由下式得到Nz＝NFFT-Nx，Nx为实信号x(n)的长度， 表示大于“·”的最小整数，NFFT是进行FFT运算的长度；步骤(2.2)取所述Fx(n)的前半部分，n＝1，……，NFFT/2，得到信号s(n)的离散傅立叶变换Fs(n)；步骤(2.3)把Fs(n)的前半部分，n＝1，……，NFFT/4，和后半部分n＝NFFT/4+1，……，NFFT/2调换；步骤(2.4)计算步骤(2.3)得到的Fs(n)的IFFT，得到解析信号s(n)，s(n)＝IFFT(Fs(n))；步骤(2.5)把步骤(2.4)得到的解析信号s(n)的 部分截去，剩余部分即为后续分析用的解析信号，其有效长度为 步骤(3)根据预定义门限，分段检测s(n)的幅值，去掉弱信号段，其步骤依次如下步骤(3.1)把步骤(2.5)得到的信号s(n)＝x(n)+jy(n)分成等长的Nseg段，用一个序列Si，i＝1，2，……，Nseg表示，Nseg＝5～20；步骤(3.2)计算步骤(3.1)所述的每一段信号幅值的和mimi=Σs(n)∈Si|s(n)|;]]>步骤(3.3)设定门限τmτm=12max{mi,i=1,2,······,Nseg};]]>步骤(3.4)将mi＜τm的信号段判定为弱信号段，去掉所有弱信号段，把剩余的信号段连接成新的解析信号序列s(n)，长度为N；步骤(4)计算步骤(3.4)所述解析信号s(n)的归一化瞬时频率序列fn(n)，按以下步骤实现步骤(4.1)按下式计算解析信号序列s(n)＝x(n)+jy(n)，n＝1，2，……，N的瞬时相位序列(n) 步骤(4.2)对步骤(4.1)得到的瞬时相位序列(n)做差分运算后取余得到信号的瞬时频率序列f(n)f(n)＝mod((n+1)-(n)+π，2π)，n＝1，2，……，N-1；步骤(4.3)把该瞬时频率序列f(n)归一化为零均值、单位方差序列fn(n)fn(n)=f(n)-mfσf,]]>其中，mf=1N-1Σn=1N-1f(n),]]>σf=1N-2Σn=1N-1(f(n)-mf)2;]]>步骤(5)计算步骤(4.3)得到的fn(n)的归一化瞬时频率的变化率序列Δfn(n)及其绝对值的均值mΔfn，其中Δfn(n)=|fn(n+2)-fn(n)|,n=1,2,······,N-3,mΔfn=1N-3Σn=1N-3Δfn(n)]]>步骤(6)把所述mΔfn同预定义门限τm比较若mΔfn<τm,]]>则把所述信号x(n)判定为MFSK信号，否则，便非MFSK信号；步骤(7)依次执行以下步骤以利用均值滤波器来减小瞬时频率的抖动步骤(7.1)利用均值滤波器对瞬时频率值序列进行平滑滤波，得到更新的瞬时频率序列fn(n)，滤波器长度lf由被分析信号可能的最大码元速率fdmax和信号的设定采样频率fs按下式求出lf=mfsfdmax,]]>m=13~12;]]>步骤(7.2)根据步骤(7.1)得到的更新的瞬时频率序列重新计算信号的归一化瞬时频率变化率序列Δfn(n)Δfn(n)＝|fn(n+2)-fn(n)|；步骤(7.3)根据下式计算步骤(7.2)得到的归一化瞬时频率变化率序列的方差τΔfτΔfn=1N-4Σn=1N-3(Δfn(n)-mΔfn)2,]]>mΔfn=1N-3Σn=1N-3Δfn(n);]]>步骤(7.4)搜索Δfn(n)序列，若|Δfn(n)|<τΔfn,]]>则fn(n+1)属符号稳定区间Γsi，i＝1，2，……，Ns(Ns为符号稳定区间的个数)，若|Δfn(n)|>τΔfn]]>则fn(n+1)属符号跳变区间Γti，i＝1，2，……，Nt(Nt为符号跳变区间的个数)；步骤(7.5)把步骤(7.4)中每一个符号稳定区间的瞬时频率值改为该区间所有瞬时频率值的均值，舍弃所有符号跳变区间的瞬时频率值，构成新的瞬时频率值序列fn(n)；步骤(8)根据“熵原则”按以下步骤分析步骤(7.5)所述瞬时频率序列值的直方图步骤(8.1)根据归一化瞬时频率序列fn(n)建立分布直方图 以此作为fn(n)的概率密度分布函数的估计p^(fi)=count{fn(n),fi≤fn(n)<fi+Δf}Nfn,fi=min{fn(n)}+iΔf,i=0,1,······,Np^-1,]]>其中，Nfn为fn(n)的长度，Δf=(max{fn(n)}-min{fn(n)})/Np^,]]> 为分布直方图的长度，Np^=512~4096,]]>count{x，y}表示满足条件y的x值的个数；步骤(8.2)确定门限τp的搜索范围为0~max{p^(fi)}/2,]]>搜索步长Δτp=max{p^(fi)}/2N,]]>N＝10～100；步骤(8.3)对搜索范围内的每个τp值执行一下步骤步骤(8.3.1)搜索满足p^(fi)>τp]]>的fi，把连续的fi归为一组；步骤(8.3.2)计算每一组的概率密度函数估计 之和，得到MFSK信号各个码元的分布概率的估计值 步骤(8.3.3)计算对应该τp的熵值ee=-Σi=1M^P^klog(P^k),]]> 是连续的fi值组成的组的个数；步骤(8.4)取对应熵值最大的τp，将相应的fi组的个数作为码元状态个数 的初始估计，每组的最小瞬时频率 和最大瞬时频率 作为每个码元状态的频率范围(f^minkf^maxk),k=1,2,······,M^;]]>步骤(9)根据“等概率分布原则”，按以下步骤精确估计码元状态个数 的精确值步骤(9.1)设定门限0<τPmin<1;]]>步骤(9.2)去掉P^k<τPminmax{P^k,k=1,2,······,M^}]]>的状态；步骤(9.3)重新计算码元状态个数 的初始估计和每个码元状态的频率范围(f^mink,f^maxk),k=1,2,······,M^;]]>步骤(10)根据“等频差原则”精确估计码元状态个数 和各状态中心频率，其步骤如下步骤(10.1)计算各个码元状态的中心频率 f^mk=Σf=fminkfmaxkp^(f)fΣf=fminkfmaxkp^(f);]]>步骤(10.2)计算各个状态之间的频差 f^dkl=|f^mk-f^ml|;]]>步骤(10.3)计算每个频差的权重wkl权重等于得到该频差的两个码元状态的概率分布估计的较小值，减去这两个码元状态之间的所有码元状态的概率分布估计之和，如果差值小于0，则权重为0wkl=max{0,min{P^k,P^l}-Σn=k+1l-1P^n},k<l-1max{0,min{P^k,P^l}},k=l-1;]]>步骤(10.4)对各个频差加权平均，得到频差的估计f^d=Σk,lwklf^dldΣk,lwkl]]>步骤(10.5)设定门限τmin、τmaxτmin=0.4f^d~0.8f^d,]]>τmax=1.2f^d~1.6f^d;]]>步骤(10.6)对于f^dk,k+1<τmin]]>的相邻状态，合并码元状态k和k+1；步骤(10.7)对于f^dk,k+1>τmax]]>的相邻状态，在码元状态k和k+1之间插入 个码元状态( 表示大于“·”的最小整数)；步骤(10.8)更新码元状态个数 和各个码元状态的中心频率 步骤(11)后续处理按下式计算得到MFSK信号的调制阶数 M^0=2[log2M^],]]>其中[·]运算表示取离“·”最近的整数。
全文摘要
本发明属于通信信号自动识别和分析技术领域，目的在于检测和分析多进制频移键控信号。本发明特征在于，在检测时用归一化后的瞬时频率的变化率序列的绝对值的平均值与设定的门限做比较，来判定MFSK信号；在分析时，用平滑滤波来减少MFSK信号瞬时频率的抖动，再用“熵原则”分析瞬时频率分布直方图，据此用等概率分布原则精确估计码元状态个数和频率范围，在此基础上，再运用“等频差”原则得到更精确的码元状态个数、频率范围和中心频率的估计。通过计算机仿真和实际测试证明本发明的检测正确率在95％以上，分析正确率在90％以上。
文档编号H04L27/156GK1946069SQ20061011408
公开日2007年4月11日 申请日期2006年10月27日 优先权日2006年10月27日
发明者邹红星, 俞力杰, 王普 申请人:清华大学