[发明专利]一种多进制频移键控信号的检测和分析方法

权利要求书

1.一种多进制频移键控信号的检测和分析方法，其特征在于，所述方法是在计算机中依次按

以下步骤实现的：

步骤(1)接收待处理数据，得到实信号x(n)；

步骤(2)按以下步骤把实信号x(n)转换成解析信号s(n)：

步骤(2.1)把实信号x(n)补零后计算FFT，得到信号的离散傅立叶变换F_x(n)，补零的长度N_z由下式得到：

N_z＝N_FFT-N_x，N_x为实信号x(n)的长度，表示大于“·”的最小整数，N_FFT是进行FFT运算的长度；

步骤(2.2)取所述F_x(n)的前半部分，n＝1，……，N_FFT/2，得到信号s(n)的离散傅立叶变换F_s(n)；

步骤(2.3)把F_s(n)的前半部分，n＝1，………，N_FFT/4，和后半部分n＝N_FFT/4+1，………，N_FFT/2调换；

步骤(2.4)计算步骤(2.3)得到的F_s(n)的IFFT，得到解析信号s(n)，s(n)＝IFFT(F_s(n))；

步骤(2.5)把步骤(2.4)得到的解析信号s(n)的部分截去，剩余部分即为后续分析用的解析信号，其有效长度为

步骤(3)根据预定义门限，分段检测s(n)的幅值，去掉弱信号段，其步骤依次如下：

步骤(3.1)把步骤(2.5)得到的信号s(n)＝x(n)+jy(n)分成等长的N_seg段，用一个序列S_i，i＝1，2，……，N_seg表示，N_seg＝5～20；

步骤(3.2)计算步骤(3.1)中得到的每一段信号幅值的和m_i：

mi=Σs(n)∈Si|s(n)|;]]>

步骤(3.3)设定门限τ_m：

τm=12max{mi,i=1,2,······,Nseg};]]>

步骤(3.4)将m_i<τ_m的信号段判定为弱信号段，去掉所有弱信号段，把剩余的信号段连接成新的解析信号序列s(n)，长度为N；

步骤(4)计算步骤(3.4)所述解析信号序列s(n)的归一化瞬时频率序列f_n(n)，按以下步骤实现：

步骤(4.1)按下式计算解析信号序列s(n)＝x(n)+jy(n)，n＝1，2，……，N的瞬时相位序列

n＝1，2，……，N；

步骤(4.2)对步骤(4.1)得到的瞬时相位序列做差分运算后取余得到信号的瞬时频率序列f(n)：

n＝1，2，……，N-1；

步骤(4.3)把该瞬时频率序列f(n)归一化为零均值、单位方差序列f_n(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)得到的f_n(n)的归一化瞬时频率的变化率序列Δf_n(n)及其绝对值的均值其中：

Δfn(n)=|fn(n+2)-fn(n)|,]]>n＝1，2，……，N-3，mΔfn=1N-3Σn=1N-3Δfn(n)]]>

步骤(6)把所述同预定义门限τ_m比较：若MΔfn<τm,]]>则把所述信号x(n)判定为MFSK信号，并且执行步骤(7)，否则，判断接收序列不是MFSK信号，流程结束；

步骤(7)依次执行以下步骤以利用均值滤波器来减小瞬时频率的抖动：

步骤(7.1)利用均值滤波器对步骤(4.3)得到的瞬时频率序列f_n(n)进行平滑滤波，得到更新的瞬时频率序列f_n(n)，滤波器长度l_f由被分析信号可能的最大码元速率f_dmax和信号的设定采样频率f_s按下式求出：

lf=m.fsfdmax,]]>m=13~12;]]>

步骤(7.2)根据步骤(7.1)得到的更新的瞬时频率序列重新计算信号的归一化瞬时频率变化率序列Δf_n(n)：Δf_n(n)＝|f_n(n+2)-f_n(n)|；

步骤(7.3)根据下式计算步骤(7.2)得到的归一化瞬时频率变化率序列的方差

τΔfn=1N-4Σn=1N-3(Δfn(n)-mΔfn)2,]]>mΔfn=1N-3Σn=1N-3Δfn(n);]]>

步骤(7.4)搜索Δf_n(n)序列，若|Δfn(n)|<τΔfn,]]>则f_n(n+1)属符号稳定区间Г_si，i＝1，2，……，N_s，其中N_s为符号稳定区间的个数，若|Δfn(n)|>τΔfn]]>则f_n(n+1)属符号跳变区间Г_ti，i＝1，2，……，N_t，其中N_t为符号跳变区间的个数；

步骤(7.5)把步骤(7.4)中每一个符号稳定区间的瞬时频率值改为该区间所有瞬时频率值的均值，舍弃所有符号跳变区间的瞬时频率值，构成新的瞬时频率序列f_n(n)；

步骤(8)根据“熵原则”按以下步骤分析步骤(7.5)所述瞬时频率序列值的直方图：

步骤(8.1)根据归一化瞬时频率序列f_n(n)建立分布直方图以此作为f_n(n)的概率密度分布函数的估计：

p^(fi)=count{fn(n),fi≤fn(n)<fi+Δf}Nfn,]]>f_i＝min{f_n(n)}+iΔf，i＝0，1，……，

其中，为f_n(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)}/(2K),]]>K＝10～100；

步骤(8.3)对搜索范围内的每个τ_p值执行以下步骤：

步骤(8.3.1)搜索满足p^(fi)>τp]]>的f_i，把连续的f_i归为一组；

步骤(8.3.2)计算每一组的概率密度分布函数估计之和，得到MFSK信号各个码元的分布概率的估计值

步骤(8.3.3)计算对应该τ_p的熵值e：

e=-Σk=1M^P^klog(P^k),]]>是连续的f_i值组成的组的个数；

步骤(8.4)取对应熵值最大的τ_p，将相应的f_i组的个数的值作为码元状态个数M的初始估计值，每组的最小瞬时频率f_mink和最大瞬时频率f_maxk作为每个码元状态的频率范围(f_mink，f_maxk)，k＝1，2，……，

步骤(9)根据“等概率分布原则”，按以下步骤精确估计码元状态个数M的精确值：

步骤(9.1)设定门限0<τPmin<1;]]>

步骤(9.2)去掉P^k<τPmin×max{P^k,k=1,2,······,M^}]]>的状态，此处，是一个取值范围在0～1的门限，max{P^k,k=1,2,······,M^}]]>是中的最大值，二者相乘得到一个用于筛选合理的码元状态的值；

步骤(9.3)重新计算码元状态个数M的估计值和每个码元状态的频率范围(f_mink，f_maxk)，k＝1，2，……，

步骤(10)根据“等频差原则”精确估计码元状态个数M和各状态中心频率，其步骤如下：

步骤(10.1)计算各个码元状态的中心频率

f^mk=Σi=iminkimaxkp^(fi)fiΣi=iminkimaxkp^(fi)]]>

其中，i_mink和i_maxk是每组码元状态的最小、最大瞬时频率的标号，是步骤(8.1)计算得到的瞬时频率序列的概率密度分布函数的估计；

步骤(10.2)计算各个状态之间的频差f^dk,l=|f^mk-f^ml|;]]>

步骤(10.3)计算每个频差的权重w_kl：

权重等于得到该频差的两个码元状态的概率分布估计的较小值，减去这两个码元状态之间的所有码元状态的概率分布估计之和，如果差值小于0，则权重为0：

wkl=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^dk,lΣ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.2)中得到的相邻码元状态k和k+1之间的频差；

步骤(10.7)对于f^dk,k+1>τmax]]>的相邻状态，在码元状态k和k+1之间插入个码元状态，其中，表示大于“·”的最小整数，表示步骤(10.2)中得到的相邻码元状态k和k+1之间的频差；

步骤(10.8)更新码元状态个数的估计和各个码元状态的中心频率

步骤(11)后续处理：

按下式计算得到MFSK信号调制阶数的最终估计

M^0=2[log2M^],]]>

其中[·]运算表示取离“·”最近的整数。