本发明属于信号处理技术领域，尤其涉及一种基于沿特征提取的脉冲信号到达时间与脉宽估计方法。
背景技术：
脉冲信号的到达时间与脉宽是表征脉冲信号时域特征的重要参数，对含噪声脉冲信号的到达时间与脉宽进行估计是信号处理领域的一个重要研究课题，在声纳、雷达以及电子战等领域具有重要的理论和应用价值，尤其在雷达和声纳信号处理中所占的地位更为突出。
目前国内外学者提出了许多脉冲信号到达时间与脉宽估计方法，主要有时域能量检测法，短时傅里叶变换法(shorttimefouriertransform，stft)，分段fft法，haar小波变换边缘检测法，自相关算法和脉冲包络边缘检测法。其中时域能量检测法无需任何先验条件适用范围最广，在高信噪比条件下可接近参数估计的克拉美罗下界(cramer-raolowbound，crlb)，但该方法对信噪比要求高，且属于非线性的方法，存在小信号抑制效应；短时傅里叶变换法实现简单，稳健性高，但该方法估计精度取决于stft的时间步进，适用于对到达时间和脉宽估计精度要求不高的场合；分段stft法，先将接收信号分段作fft，检测信号并确定信号起点所在的分段，然后将该段的fft数据取平方，再作fft，取模的极大值作为到达时间估计值，但平方运算引入了非线性运算，低信噪比时性能不佳；haar小波变换边缘检测法利用相关检测算法，对信号的起止时间进行粗估计，在检测到信号的条件下，估计出信号频率并将其变换至基带，在一定尺度下对基带信号作haar小波变换，检测出小波变换模值的峰点位置作为到达时间精估计，该方法充分利用了haar小波变换具有边缘检测及定位信号突变点的能力，但脉冲信号到达时间与脉宽估计是一种典型的双阶跃型边缘检测问题，参数估计精度对小波尺度选择的合理性十分敏感；自相关算法利用脉冲信号与宽带噪声的时间相关性差异，通过相关运算，达到抑制噪声、检测信号的目的，适用于调频脉冲信号的到达时间与脉宽的估计，但由于单频脉冲信号自相关峰不够尖锐，自相关峰检测困难，该方法对单频脉冲信号的到达时间与脉宽估计适用性差。
脉冲包络边缘检测法充分利用了脉冲包络沿的瞬变特征，近年来得到了广泛的关注与研究，该类方法在高信噪比条件下估计精度较高，但传统的脉冲包络边缘检测法往往存在两方面的不足：(1)直接利用希尔伯特变换提取脉冲信号包络，几乎没有任何处理增益，对脉冲信号的信噪比要求高；(2)直接利用脉冲信号包络的一阶或二阶差分作为脉冲信号包络前后沿瞬变特征的评价因子，提取脉冲信号前后沿，当脉冲信号经过信道传播，导致接收信号存在多途或混响干扰，从而使脉冲信号包络发生畸变时，该类方法不再适用。
技术实现要素：
发明目的：本发明提供一种基于沿特征提取的脉冲信号到达时间与脉宽估计方法，在满足雷达和声纳信号处理需求的前提下，可起到抗多途和混响干扰的作用，克服统脉冲包络边缘检测法的不足，在运算量相当的条件下，脉冲信号到达时间与脉宽估计性能较传统的脉冲包括边缘检测法有明显改进。
技术方案：为了实现上述发明目的，本发明采用以下技术方案：一种基于沿特征提取的脉冲信号到达时间与脉宽估计方法，包括如下步骤：
(1)获取待处理的脉冲信号采样数据序列：从传感器接收n个采样点的实时采集数据作为待处理的数据序列x(n)，n＝0，1，…，n-1，或从存储器中提取包含整个脉冲信号的n个采样点数据作为待处理的数据序列x(n)，n＝0，1，…，n-1，所述的n为包含整个脉冲信号的采样点数，取值为2的整数次幂；
(2)提取脉冲信号包络：对所述数据序列x(n)依次做离散傅里叶变换，频域带通滤波，离散傅里叶逆变换和取模运算得到其脉冲信号包络数据序列y(n)；
(3)脉冲包络归一化与差值序列计算：对y(n)进行归一化处理得到归一化的脉冲信号包络数据序列z(n)和差值序列σ(n)；
(4)综合评价因子及其各子参数初始化：对脉冲包络上升沿和下降沿综合评价因子及其各子参数进行初始化设置；
(5)计算脉冲包络上升沿综合评价因子：计算归一化包络数据序列z(n)的上升沿综合评价因子u(k)；
(6)计算脉冲包络下降沿综合评价因子：计算归一化包络数据序列z(n)的下降沿综合评价因子d(k)；
(7)估计出脉冲信号到达时间和脉宽：分别搜索u(k)和d(k)的最大值所对应的离散采样时间索引ku和kd，依据ku和kd估计出脉冲信号到达时间ts和脉宽τ。
优选的，步骤(2)中，对所述数据序列x(n)依次做离散傅里叶变换，频域带通滤波，离散傅里叶逆变换和取模运算得到其脉冲信号包络数据序列y(n)，具体包括如下五步：
(2-1)对数据序列x(n)做离散傅里叶变换，x(n)的离散傅里叶变换x(l)为：
其中l为x(l)的离散频率索引，j表示虚数单位，即
(2-2)依据已知的脉冲信号下限频率fl和上限频率fh产生滤波器系数h(l)：
其中kl和kh分别为频域带通滤波器下限频率bl＝fl-ξfl和上限频bh＝fh+ξfh所对应的离散频率索引，ξ为滤波器带宽系数，0＜ξ＜1，kl和kh的值分别为：
其中δf为长度为n的离散傅里叶变换的频率分辨率，δf＝fs/n，fs为采样频率，int()代表四舍五入取整运算；
(2-3)对x(l)做频域带通滤波得到滤波后的离散频谱y(l)：
y(l)＝x(l)h(l)，l＝0，1…，n-1
(2-4)对y(l)做离散傅里叶逆变换得脉冲信号滤波后的时域复数信号yc(n)：
(2-5)取时域复数信号yc(n)的模值得到脉冲信号包络数据序列y(n)：
y(n)＝|yc(n)|，n＝0，1…，n-1
其中||代表取模值运算。优选的，步骤(2)中，步骤(2-1)计算x(n)的离散傅里叶变换x(l)和步骤(2-4)对y(l)做离散傅里叶逆变换得时域复数信号yc(n)均可通过快速傅里叶变换实现；步骤(2-2)中滤波器带宽系数ξ的优选值为0.1。
优选的，步骤(3)中，对包络y(n)进行归一化处理得到归一化的脉冲信号包络数据序列z(n)和差值序列σ(n)，z(n)和σ(n)分别为：
其中min[y(n)]和max[y(n)]分别为y(n)在n满足0≤n≤n-1范围内的最小值和最大值；
优选的，步骤(4)中，对脉冲包络上升沿和下降沿综合评价因子及其各子参数进行初始化设置，具体包括如下参数的初始化：
脉冲包络上升沿综合评价因子u(k)初始化为：u(k)＝0，k＝0，1…，n-1
脉冲包络上升沿瞬变中心权重u1(k)初始化为：u1(k)＝0，k＝0，1…，n-1
脉冲包络上升沿连续上升权重u2(k)初始化为：u2(k)＝0，k＝0，1…，n-1
脉冲包络下降沿综合评价因子d(k)初始化为：d(k)＝0，k＝0，1…，n-1
脉冲包络下降沿瞬变中心权重d1(k)初始化为：d1(k)＝0，k＝0，1…，n-1
脉冲包络下降沿连续下降权重d2(k)初始化为：d2(k)＝0，k＝0，1…，n-1
其中k为上述各参量的离散时间索引；
脉冲沿特征搜索窗长下限和上限wl和wh，wl和wh需同时满足以下三个条件：
1≤wl≤n-1，1≤wh≤n-1，wl≤wh
脉冲沿瞬变中心权重迭代搜索截止门限η1，0≤η1≤1；
脉冲沿连续单调变化权重迭代搜索截止门限η2，0≤η2≤1；
优选的，步骤(4)中，为了兼顾本发明的运算量和估计精度，脉冲沿特征搜索窗长下限和上限wl和wh的优选值分别为10和50，脉冲沿瞬变中心权重迭代搜索截止门限η1优选值为0.5，脉冲沿连续单调变化权重迭代搜索截止门限η2的优选值为0.75。
优选的，步骤(5)中，计算归一化包络z(n)的上升沿综合评价因子u(k)，具体包括如下六步：
(5-1)初始化上升沿搜索离散采样时间索引ks＝wl-1和上升沿搜索窗长wu＝2wl+1；
(5-2)令ks＝ks+1，u1(ks)＝1，u2(ks)＝1，并重新计算u1(ks)：
u1(ks)＝1-2|z(ks)-0.5|
(5-3)重新计算u2(ks)：
u2(ks)＝cu/wu
其中，cu为σ(n)在n满足
(5-4)令wu＝wu+1，并判断是否满足以下条件：
若满足返回步骤(5-3)，否则转入步骤(5-5)；
(5-5)判断ks≤n-wl是否成立，若成立返回(5-2)；否则进入步骤(5-6)；
(5-6)计算上升沿综合评价因子u(k)：
u(k)＝u1(k)u2(k)，k＝0，1…，n-1。
优选的，步骤(6)中，计算归一化包络z(n)的下降沿综合评价因子d(k)，具体包括如下六步：
(6-1)初始化下降沿搜索离散采样时间索引ke＝wl-1和下降沿搜索窗长wd＝2wl+1
(6-2)令ke＝ke+1，d1(ke)＝1，d2(ke)＝1，并重计算d1(ke)
d1(ke)＝1-2|z(ke)-0.5|
(6-3)重计算d2(ke)
d2(ke)＝cd/wd
其中cd为σ(n)在n满足
(6-4)令wd＝wd+1，并判断是否满足以下条件：
若满足返回步骤(6-3)，否则转入步骤(6-5)；
(6-5)判断ke≤n-wl是否成立，若成立返回(6-2)；否则进入步骤(6-6)；
(6-6)计算下降沿综合评价因子d(k)：
d(k)＝d1(k)d2(k)，k＝0，1…，n-1。
优选的，步骤(7)中，分别搜索u(k)和d(k)的最大值所对应的离散采样时间索引ku和kd，依据ku和kd估计出脉冲信号到达时间ts和脉宽τ，具体包括如下两步：
(7-1)分别搜索u(k)和d(k)的最大值所对应的离散采样时间索引ku和kd
其中
(7-2)依据ku和kd估计出脉冲信号到达时间τ0和脉宽τ：
τ0＝ku/fs，τ＝(kd-ku)/fs。
有益效果：本发明与现有的方法相比，有以下几点有益效果：
(1)本发明的估计方法在提取脉冲信号时，利用频域带通滤波技术对传统希尔伯特变换进行了改进，充分利用了脉冲信号的频域特征，可获得近似匹配滤波的处理增益，较常规希尔伯特变换包络提取法，扩展了可适用信噪比的下限；
(2)本发明的估计方法在对脉冲包络边缘进行检测时，设计了脉冲沿综合评价因子，该评价因子通过脉冲包络沿瞬变中心权重和单调变化权重的综合加权，较传统脉冲包络边缘检测法，本发明方法具备更强的抗多途和混响干扰的能力；
(3)本发明的估计方法在对脉冲包络上升沿和下降沿综合评价因子及其各子参数进行初始化设置时，可依据实际使用中对运算量和估计精度的不同要求，动态调整脉冲沿特征搜索窗长下限和上限wl和wh：wl越小，wh越大，估计精度越高，但相应的运算量也会增大；反之，wl越大，wh越小，运算量降低，但估计精度也会变差。
附图说明
图1为本发明方法的流程图；
图2为实施例1的仿真脉冲信号包络图；
图3为实施例1的仿真脉冲信号归一化包络图；
图4为实施例1的仿真脉冲信号归一化包络差值序列图；
图5为实施例1的上升沿综合评价因子图；
图6为实施例1的下降沿综合评价因子图；
图7为实施例2的仿真脉冲信号包络图；
图8为实施例2的仿真脉冲信号归一化包络图；
图9为实施例2的仿真脉冲信号归一化包络差值序列图；
图10为实施例2的上升沿综合评价因子图；
图11为实施例2的下降沿综合评价因子图。
具体实施方式
下面结合附图和实施例对本发明做进一步的说明：
如图1所示，一种基于沿特征提取的脉冲信号到达时间与脉宽估计方法，包括如下步骤：
(1)获取待处理的脉冲信号采样数据序列：从传感器接收n个采样点的实时采集数据作为待处理的数据序列x(n)，n＝0，1，…，n-1，或从存储器中提取包含整个脉冲信号的n个采样点数据作为待处理的数据序列x(n)，n＝0，1，…，n-1，所述的n为包含整个脉冲信号的采样点数，取值为2的整数次幂。
(2)提取脉冲信号包络：对所述数据序列x(n)依次做离散傅里叶变换，频域带通滤波，离散傅里叶逆变换和取模运算得到其脉冲信号包络数据序列y(n)，具体包括如下五步：
(2-1)对数据序列x(n)做离散傅里叶变换，x(n)的离散傅里叶变换x(l)为
其中l为x(l)的离散频率索引，j表示虚数单位，即
(2-2)依据已知的脉冲信号下限频率fl和上限频率fh产生滤波器系数h(l)：
其中kl和kh分别为频域带通滤波器下限频率bl＝fl-ξfl和上限频bh＝fh+ξfh所对应的离散频率索引，ξ为滤波器带宽系数，0＜ξ＜1，kl和kh的值分别为：
其中δf为长度为n的离散傅里叶变换的频率分辨率，δf＝fs/n，fs为采样频率，int()代表四舍五入取整运算；
(2-3)对x(l)做频域带通滤波得到滤波后的离散频谱y(l)：
y(l)＝x(l)h(l)，l＝0，1…，n-1
(2-4)对y(l)做离散傅里叶逆变换得脉冲信号滤波后的时域复数信号yc(n)：
(2-5)取复数yc(n)的模值得到脉冲信号包络数据序列y(n)：
y(n)＝|yc(n)|，n＝0，1...，n-1
其中||代表取模值运算。
在第(2)步中，步骤(2-1)计算x(n)的离散傅里叶变换x(l)和步骤(2-4)对y(l)做离散傅里叶逆变换得时域复数信号yc(n)均可通过快速傅里叶变换实现，提高算法的计算效率；步骤(2-2)中滤波器带宽系数ξ的优选值为0.1。
(3)脉冲包络归一化与差值序列计算：对y(n)进行归一化处理得到归一化的脉冲信号包络数据序列z(n)和差值序列σ(n)，计算过程如下：
(4)综合评价因子及其各子参数初始化：对脉冲包络上升沿和下降沿综合评价因子及其各子参数进行初始化设置，具体包括如下参数的初始化：
脉冲包络上升沿综合评价因子u(k)初始化为：u(k)＝0，k＝0，1…，n-1
脉冲包络上升沿瞬变中心权重u1(k)初始化为：u1(k)＝0，k＝0，1…，n-1
脉冲包络上升沿连续上升权重u2(k)初始化为：u2(k)＝0，k＝0，1…，n-1
脉冲包络下降沿综合评价因子d(k)初始化为：d(k)＝0，k＝0，1…，n-1
脉冲包络下降沿瞬变中心权重d1(k)初始化为：d1(k)＝0，k＝0，1…，n-1
脉冲包络下降沿连续下降权重d2(k)初始化为：d2(k)＝0，k＝0，1…，n-1
其中k为上述各参量的离散时间索引；
脉冲沿特征搜索窗长下限和上限wl和wh，wl和wh需同时满足以下三个条件：
1≤wl≤n-1，1≤wh≤n-1，wl≤wh
脉冲沿瞬变中心权重迭代搜索截止门限η1，0≤η1≤1；
脉冲沿连续单调变化权重迭代搜索截止门限η2，0≤η2≤1；
在第(4)步中，为了兼顾本发明的运算量和估计精度，脉冲沿特征搜索窗长下限和上限wl和wh的优选值分别为10和50，脉冲沿瞬变中心权重迭代搜索截止门限η1优选值为0.5，脉冲沿连续单调变化权重迭代搜索截止门限η2的优选值为0.75。
(5)计算脉冲包络上升沿综合评价因子：计算归一化包络z(n)的上升沿综合评价因子u(k)，具体包括如下六步：
(5-1)初始化上升沿搜索离散采样时间索引ks＝wl-1和上升沿搜索窗长wu＝2wl+1；
(5-2)令ks＝ks+1，u1(ks)＝1，u2(ks)＝1，并重新计算u1(ks)：
u1(ks)＝1-2|z(ks)-0.5|
(5-3)重新计算u2(ks)：
u2(ks)＝cu/wu
其中cu为σ(n)在n满足
(5-4)：令wu＝wu+1，并判断是否满足以下条件：
若满足返回步骤(5-3)，否则转入步骤(5-5)；
(5-5)判断ks≤n-wl是否成立，若成立返回(5-2)；否则进入步骤(5-6)；
(5-6)计算上升沿综合评价因子u(k)：
u(k)＝u1(k)u2(k)，k＝0，1…，n-1。
(6)计算脉冲包络下降沿综合评价因子：计算归一化包络z(n)的下降沿综合评价因子d(k)，具体包括如下六步：
(6-1)初始化下降沿搜索离散采样时间索引ke＝wl-1和下降沿搜索窗长wd＝2wl+1
(6-2)令ke＝ke+1，d1(ke)＝1，d2(ke)＝1，并重计算d1(ke)
d1(ke)＝1-2|z(ke)-0.5|
(6-3)重计算d2(ke)
d2(ke)＝cd/wd
其中cd为σ(n)在n满足
若满足返回步骤(6-3)，否则转入步骤(6-5)；
(6-5)判断ke≤n-wl是否成立，若成立返回(6-2)；否则进入步骤(6-6)；
(6-6)计算下降沿综合评价因子d(k)：
d(k)＝d1(k)d2(k)，k＝0，1…，n-1。
(7)估计出脉冲信号到达时间和脉宽：分别搜索u(k)和d(k)的最大值所对应的离散采样时间索引ku和kd，依据ku和kd估计出脉冲信号到达时间ts和脉宽τ，具体包括如下两步：
(7-1)分别搜索u(k)和d(k)的最大值所对应的离散采样时间索引ku和kd
其中
(7-2)依据ku和kd估计出脉冲信号到达时间τ0和脉宽τ
τ0＝ku/fs，τ＝(kd-ku)/fs。
本发明的实施例中，仿真接收脉冲信号模型为：
其中a为信号幅度，
以采样频率fs对上述脉冲信号进行离散采样可得到脉冲信号采样数据序列：
其中n0＝int(τ0fs)，m＝int(τfs)，n＝int(tfs)。
实施例1：
仿真信号参数分别设置为：信号幅度a＝2，初始相位
下面对仿真脉冲进行到达时间与脉宽估计：
在第(1)步中，设置滤波器带宽系数ξ＝0.1，由第(2)步提取的仿真脉冲信号包络如图2所示；
依据第(3)步，对y(n)进行归一化处理得到归一化的脉冲信号包络数据序列z(n)和差值序列σ(n)分别如图3和图4所示；
在第(4)步中，初始化wl＝10，wh＝50，η1＝0.5，η2＝0.75；
依据第(5)步，计算归一化包络z(n)的上升沿综合评价因子u(k)，如图5所示；
依据第(6)步，计算归一化包络z(n)的下降沿综合评价因子d(k)，如图6所示；
由图5和图6可以看出，u(k)和d(k)的最大值所对应的离散采样时间索引ku＝1023，kd＝5133，由此估计出脉冲信号的到达时间τ0和脉宽τ分别为：
估计相对误差分别为
实施例2：
仿真信号参数分别设置为：信号幅度a＝3，初始相位
下面对仿真脉冲进行到达时间与脉宽估计：
在第(1)步中，设置滤波器带宽系数ξ＝0.1，由第(2)步提取的仿真脉冲信号包络如图7所示；
依据第(3)步，对y(n)进行归一化处理得到归一化的脉冲信号包络数据序列z(n)和差值序列σ(n)分别如图8和图9所示；
在第(4)步中，初始化wl＝15，wh＝45，η1＝0.45，η2＝0.8；
依据第(5)步，计算归一化包络z(n)的上升沿综合评价因子u(k)，如图10所示；
依据第(6)步，计算归一化包络z(n)的下降沿综合评价因子d(k)，如图11所示；
由图10和图11可以看出，u(k)和d(k)的最大值所对应的离散采样时间索引ku＝263，kd＝1276，由此估计出脉冲信号的到达时间τ0和脉宽τ分别为：
到达时间τ0和脉宽τ估计相对误差分别为：