[发明专利]基于频率拼接的异步跳频信号参数盲估计方法

权利要求书

1.基于频率拼接的异步跳频信号参数盲估计方法，包括如下步骤：

(1)接收信号：

雷达天线接收任意一段模拟信号；

(2)信号采样：

对雷达天线接收的任意一段模拟信号，以采样频率f_s进行采样，得到离散信号；

(3)对离散信号进行短时傅里叶变换，得到时频域信号：

(4)提取序号：

(4a)提取时频域信号中高能量信号对应的频域采样点序号和时域窗序号；

(4b)将高能量信号对应的频域采样点序号和时域窗序号组成序号集合；

(5)计算每一个频域采样点序号对应的频率值：

(6)聚类频率：

(6a)采用K均值聚类算法，分别对每一个时域窗对应的所有频域采样点序号l₀对应的频率值f(l₀)进行聚类，得到第q₀个时域窗中第c个聚类中心其中，l₀表示时频域信号中高能量信号对应的频域采样点序号，l₀＝1,2,…,L，L表示短时傅里叶变换窗长，f(l₀)表示第q₀个时域窗对应的频域采样点序号l₀对应的频率值，q₀表示时频域信号中高能量信号对应的时域窗序号，q₀＝1,2,…,Q，Q表示短时傅里叶变换所需的总窗数，表示第q₀个时域窗中第c个聚类中心，c表示聚类中心的序号，表示聚类中心的个数；

(6b)将第q₀个时域窗中第c个聚类中心组成频率集合其中，q₀表示时频域信号中高能量信号对应的时域窗序号，q₀＝1,2,…,Q，Q表示短时傅里叶变换所需的总窗数，表示第q₀个时域窗中第c个聚类中心，c表示聚类中心的序号，表示表示聚类中心的个数；

(7)分析频率集合：

(7a)从聚类中心组成的频率集合中提取出所有不同的频率值f_d，fd∈{fcq0,c=1,2,...,Nq0,q0=1,2,...,Q},]]>其中，d表示不同的频率点序号，d＝1,2,…,N_o，N_o表示估计的不同频率的个数，q₀表示时频域信号中高能量信号对应的时域窗序号，q₀＝1,2,…,Q，Q表示短时傅里叶变换所需的总窗数，表示第q₀个时域窗中第c个聚类中心，c表示聚类中心的序号，表示表示聚类中心的个数；

(7b)统计每个频率f_d的起始窗序号终止窗序号以及持续窗数其中，f_d表示提取出的频率，d表示不同的频率点序号，d＝1,2,…,N_o，N_o表示估计的不同频率的个数，频率f_d、起始窗序号终止窗序号以及持续窗数组成向量T_d，

(7c)提取向量T_d中起始窗序号为1的向量T_j⁽⁰⁾，逐行存入矩阵F⁽⁰⁾，其中，d表示不同的频率点序号，d＝1,2,…,N_o，N_o表示估计的不同频率的个数，j表示矩阵F⁽⁰⁾的行序号，j＝1,2,…,r⁽⁰⁾，r⁽⁰⁾表示矩阵F⁽⁰⁾的行数；

(7d)提取向量T_d中终止窗序号为Q的向量T_g⁽²⁾，逐行存入矩阵F⁽²⁾，其中，d表示不同的频率点序号，d＝1,2,…,N_o，N_o表示估计的不同频率的个数，g表示矩阵F⁽²⁾的行序号，g＝1,2,…,r⁽²⁾，r⁽²⁾表示矩阵F⁽²⁾的行数，Q表示短时傅里叶变换所需的总窗数；

(7e)提取向量T_d中起始窗序号不为1且终止窗序号不为Q的向量T_i⁽¹⁾，逐行存入矩阵F⁽¹⁾，其中，d表示不同的频率点序号，d＝1,2,…,N_o，N_o表示估计的不同频率的个数，i表示矩阵F⁽¹⁾的行序号，i＝1,2,…,r⁽¹⁾，r⁽¹⁾表示矩阵F⁽¹⁾的行数，Q表示短时傅里叶变换所需的总窗数；

(7f)将矩阵F⁽⁰⁾的行序号j初始化为1；

(7g)将矩阵F⁽⁰⁾第j行的向量T_j⁽⁰⁾对应的终止窗序号依次与矩阵F⁽¹⁾的每一行的向量T_i⁽¹⁾对应的起始窗序号相减，将差值依次放入集合Ψ中，其中，j表示矩阵F⁽⁰⁾的行序号，i表示矩阵F⁽¹⁾的行序号，i＝1,2,…,r⁽¹⁾，r⁽¹⁾表示矩阵F⁽¹⁾的行数；

(7h)从集合Ψ中搜索大于等于-1的最小差值η，若搜索到，执行步骤(7i)；若未搜索到，将矩阵F⁽⁰⁾的行序号j加1，执行步骤(7g)；

(7i)按照下式，计算同时包含跳频信号相连两跳频率的时域窗的窗数；

w＝η+1

其中，w表示同时包含跳频信号相连两跳频率的时域窗的窗数，η表示从集合Ψ中搜索的大于等于-1的最小差值；

(8)初步估计跳周期：

(8a)采用K均值聚类算法，对矩阵F⁽¹⁾中的第4列数据进行聚类，得到聚类中心o_m和每一类持续窗数的个数Z_m，每一类持续窗数的个数等于该类持续窗数对应的频率个数，其中，矩阵F⁽¹⁾中的第4列数据表示持续窗数，m表示聚类中心的序号，第m个聚类中心对应第m类信号，m＝1,2,…，N，N表示聚类中心的个数；

(8b)将每一类持续窗数对应的频率组成集合Ω_m；

(8c)按照下式，初步估计跳周期：

T~m=om×B/fs]]>

其中，表示第m类信号的跳周期，m表示聚类中心的序号，第m个聚类中心对应第m类信号，m＝1,2,…,N，N表示聚类中心的个数，o_m表示所得聚类中心，B表示短时傅里叶变换时，时域窗间隔的采样点数，f_s表示采样频率；

(9)检测碰撞频率：

(9a)提取步骤(8a)得到的持续窗数的个数为1的那类持续窗数对应的频率f^p；

(9b)判断频率f^p是否满足第一类拼接条件，若是，执行步骤(11)，否则，执行步骤(9c)；

(9c)判断频率f^p是否满足第二类拼接条件，若是，则将该频率f^p作为碰撞频率，执行步骤(10)，否则，执行步骤(11)；

(10)分解碰撞频率：

(10a)将碰撞频率f^p作为频率所在类对应的跳频信号在窗内的频率，其中，表示频率的起始窗序号，W₁表示频率对应的持续窗数，表示能够拼接在碰撞频率f^p之前的起始窗序号不为1且终止窗序号不为Q的频率，Q表示短时傅里叶变换所需的总窗数；

(10b)将碰撞频率f^p作为频率所在类对应的跳频信号在窗内的频率，其中，表示频率的终止窗序号，W₂表示频率对应的持续窗数，表示能够拼接在碰撞频率f^p之后的起始窗序号不为1且终止窗序号不为Q的频率，Q表示短时傅里叶变换所需的总窗数；

(11)构造矩阵：

将步骤(8a)得到的第m个聚类中心对应的所有频率对应的向量按起始窗序号从小到大进行排序，将排序结果逐行存入矩阵其中，m表示聚类中心的序号，第m个聚类中心对应第m类信号，m＝1,2,…,N，N表示聚类中心的个数，i₀表示矩阵的行序号，表示向量的第一个元素，Ω_m表示步骤(8b)得到的每一类持续窗数对应的频率组成的集合；

(12)拼接每一个起始窗序号不为1且终止窗序号不为Q的跳频信号的第一跳频率以及最后一跳频率，其中，Q表示短时傅里叶变换所需的总窗数：

(12a)将步骤(8a)得到的聚类中心的序号m初始化为1；

(12b)将向量中的终止窗序号依次与矩阵第一行的起始窗序号相减，得到差值Δ_j，其中，j表示矩阵F⁽⁰⁾的行序号，j＝1,2,…,r⁽⁰⁾，r⁽⁰⁾表示矩阵F⁽⁰⁾的行数；

(12c)按照下式，搜索差值Δ_j的最大值的下标：

j*=argmaxj{Δj∈[-1,w]}]]>

其中，j^*表示差值Δ_j的最大值的下标，表示取最大值所对应的下标值操作，w表示同时包含同一跳频信号相连两跳频率的时域窗的窗数；

(12d)将向量对应的频率作为矩阵第一行对应的第一跳频率；

(12e)将矩阵最后一行的终止窗序号与向量T_g⁽²⁾的起始窗序号相减，得到差值其中，g表示矩阵F⁽²⁾的行序号，g＝1,2,…,r⁽²⁾，r⁽²⁾表示矩阵F⁽²⁾的行数；

(12f)按照下式，搜索差值的最大值的下标：

g^=arg maxg{Δ~g∈[-1,w]}]]>

其中，表示差值的最大值的下标，表示取最大值所对应的下标值操作，w表示同时包含同一跳频信号相连两跳频率的时域窗的窗数；

(12g)将向量对应的频率作为矩阵最后一行对应的最后一跳频率；

(12h)按照式得到新的矩阵

(12i)将聚类中心的序号m加1，如果m≤N，返回步骤(12b)，如果m＞N，执行步骤(13)，其中，m表示聚类中心的序号，第m个聚类中心对应第m类信号，m＝1,2,…,N，N表示聚类中心的个数；

(13)按照下式，精确估计跳周期：

T^m=(1I-1Σs=1I-1(T‾m(s+1,2)-T‾m(s,2)))×B/fs]]>

其中，表示第m类信号的跳周期，m表示聚类中心的序号，第m个聚类中心对应第m类信号，m＝1,2,…,N，N表示聚类中心的个数，I表示矩阵的总行数，矩阵是步骤(12h)得到的矩阵，Σ表示求和符号，s表示矩阵行序号，表示矩阵的第s行第2列的元素，B表示短时傅里叶变换时，时域窗间隔的采样点数，f_s表示采样频率；

(14)估计跳变时刻：

(14a)按照下式，计算第m类信号第一跳的跳变时刻：

tm(1)=[T‾m(1,3)-w2]×1fs]]>

其中，表示第m类信号第一跳的跳变时刻，m表示聚类中心的序号，第m个聚类中心对应第m类信号，m＝1,2,…,N，N表示聚类中心的个数，矩阵是步骤(12h)得到的矩阵，表示矩阵第1行第3列的元素，w表示同时包含同一跳频信号相连两跳频率的时域窗的窗数，f_s表示采样频率；

(14b)按照下式，计算第m类信号第k跳的跳变时刻:

tm(k)=tm(1)+(k-1)T^m]]>

其中，表示第m类信号第k跳的跳变时刻，m表示聚类中心的序号，第m个聚类中心对应第m类信号，m＝1,2,…,N，N表示聚类中心的个数，k表示跳数的序号，k＝1,2…Z_m-1，Z_m表示步骤(8a)聚类得到的每一类持续窗数的个数；

(15)按照下式，估计第m类信号的第k跳频率：

f^mk=T‾m(k,1)]]>

其中，表示第m类信号的第k跳频率，m表示聚类中心的序号，第m个聚类中心对应第m类信号，m＝1,2,…,N，N表示聚类中心的个数，k表示跳数的序号，k＝1,2…Z_m-1，Z_m表示步骤(8a)聚类得到的每一类持续窗数的个数，矩阵是步骤(12h)得到的矩阵，表示矩阵第k行第1列的元素。