本发明属于信号处理技术领域,尤其涉及一种脉冲信号频域自适应滤波包络提取方法。
背景技术:
信号包络是表征脉冲信号时域特征的重要参数,准确提取脉冲信号包络可为进一步精确估计脉冲信号的到达时间和脉冲宽度提供基础,并可利用两路脉冲信号包络的时延差估计信号的来波方向,在声纳、雷达以及电子战等领域具有重要的理论和应用价值,尤其在雷达和声纳信号处理中所占的地位更为突出。
常用的信号包络提取方法主要有四类,分别为:(1)希尔伯特(hilbert)变换法;(2)小波变换法;(3)改进的hilbert变换法,如基于经验模式分解的hilbert变换(或称为基于希尔伯特-黄变换)的包络提取法,基于盲源提取的hilbert变换法;(4)基于数学形态学处理的包络提取法。其中hilbert变换法首先将原始信号转变为复解析信号,然后取复解析信号的模值作为信号的包络,该方法以fft为基础,实现简单,运算量小,高信噪比下包络提取效果好,但该方法是在整个频率区间上进行的,抑制噪声的能力差,且这种方法对窄带载波信号的包络提取是很有用的,但对于时变的宽带信号,hilbert变换法不可避免的存在两个缺点:①将信号中的高频成分同时提取出来,从而产生不光滑的毛刺;②不具备自适应分析能力。
小波变换法首先对信号进行小波变换得到解析信号,小波变换相当于利用一系列的带通滤波器对信号进行多通带滤波,因此得的解析信号是经过带通滤波的,由于得到的解析信号的实部和虚部正交,再取解析信号的模值即可获得信号的包络。与hilbert变换法相比,该方法提取出的包络更为理想,但该方法的难点在于小波变换尺度的选择,尺度选取的合适可以使带通滤波器覆盖频带与信号频率相匹配,包络提取效果好;尺度选取不合适使得带通滤波器覆盖频带与信号频率失配,包络提取效果急剧变差,而尺度的选取很大程度上取决于使用者的经验,通用性差。改进的hilbert变换法首先是利用某一特定的预处理方法对信号进行滤波去噪处理,然后进一步利用hilbert变换提取包络,如前述的经验模式分解和盲源提取等,这类方法可提高hilbert变换法噪声抑制能力,但当前常见的改进方法,实现复杂,且对不同的信号通用性差。基于数学形态学的处理方法,首先利用形态学滤波对信号进行预处理,然后利用形态学闭运算提取信号包络,最后应用形态学开运算来消除噪声包络,这种方法的关键是选取合适的结构元素,而对多类型的脉冲信号,很难找到一个适合各种不同信号的结构元素。
技术实现要素:
发明目的:针对现有技术存在的问题与不足,本发明提供一种脉冲信号频域自适应滤波包络提取方法,以满足雷达和声纳信号处理的需求。并且,该方法与现有hilbert变换包络提取法的运算量相当,但相较hilbert变换法,具有更好的噪声抑制能力和自适应能力,原理清晰,实现简单,工程适用性强。
技术方案:为了实现上述发明目的,本发明采用以下技术方案:一种脉冲信号频域自适应滤波包络提取方法,包括如下步骤:
(1)获取待处理的脉冲信号采样数据序列:从传感器接收n个采样点的实时采集数据作为待处理的数据序列x(n),n=0,1,…,n-1,或从存储器中提取包含整个脉冲信号的n个采样点数据作为待处理的数据序列x(n),n=0,1,…,n-1,所述的n为包含整个脉冲信号的采样点数,取值为2的整数次幂。
(2)根据所述数据序列x(n)计算脉冲信号功率谱p(k)。
(3)迭代平滑与滤波器带宽相关参数初始化:设置功率谱迭代平滑的最大迭代次数门限,精度控制门限和迭代平滑窗长分别为i1,ε1和m1,初始包络迭代平滑的最大迭代次数门限,精度控制门限和迭代平滑窗长分别为i2,ε2和m2,窄带信号带宽系数ω1和宽带信号带宽系数ω2,判断窄带脉冲和宽带脉冲主瓣特征参数门限η。
(4)对p(k)进行迭代平滑处理得到平滑功率谱s(k)。
(5)提取原始功率谱和平滑后功率谱主瓣特征参数γp和γs。
(6)生成自适应滤波参数h(k)。
(7)依据所述滤波参数h(k)对脉冲信号进行自适应滤波。
(8)提取脉冲信号包络y(n)。
优选的,在步骤(2)中,对所述数据序列x(n)做快速傅里叶变换,计算得到数据序列的离散傅里叶变换x(l)和功率谱p(k),具体包括如下步骤:
步骤2-1:计算x(n)的离散傅里叶变换为:
其中i为x(l)的离散频率索引,j表示虚数单位,即
步骤2-2:依据x(l)计算x(n)的功率谱:
其中,k为功率谱p(k)的离散频率索引,||代表取模值运算,|x(k)|为信号幅度谱。
优选的,步骤(2)中,步骤2-1计算x(n)的离散傅里叶变换x(l)可通过快速傅里叶变换实现;
优选的,步骤(3)中,设置功率谱迭代平滑的最大迭代次数门限,精度控制门限和迭代平滑窗长分别为i1,ε1和m1,初始包络迭代平滑的最大迭代次数门限,精度控制门限和迭代平滑窗长分别为i2,ε2和m2,窄带信号带宽系数ω1和宽带信号带宽系数ω2,判断窄带脉冲和宽带脉冲主瓣特征参数门限η,其中i1和i2取值为大于2的正整数,ε1和ε2取值为小于1的正数,功率谱迭代平滑窗长为满足3≤m1≤n/2-2的奇数,包络迭代平滑窗长为满足5≤m2≤49的奇数,ω1为满足0≤ω1≤0.05的任意数,ω2为满足0.05<ω2≤0.1的任意数,η为满足1.5≤η≤2.5的任意数。优选的i1取值为10,ε1取值为0.01,m1取值为5,i2取值为5,ε2取值为0.02,m2取值为11,ω1的取值为0.05,ω2的取值为0.1,η取值为2.0。
优选的,步骤(4)中,对功率谱p(k)进行迭代平滑处理得到平滑功率谱s(k),具体包括如下步骤:
步骤4-1:初始化功率谱迭代平滑次数i=0,功率谱迭代平滑结果s0(k)为:
s0(k)=p(k),k=0,1,2…,n/2-1
步骤4-2:令迭代次数i=i+1,并对功率谱的上次迭代平滑结果si-1(k)进行平滑处理得到当次平滑结果si(k):
步骤4-3:判断是否达到功率谱最大迭代平滑次数,即判断i≤i1是否成立,若成立进入步骤4-4,否则跳到步骤4-6;
步骤4-4:分别计算功率谱上次平滑结果si-1(k)和当次平滑结果si(k)与原始功率谱p(k)的残差平方和ji-1和ji,分别为
步骤4-5:判断|ji-1-ji|≤ε1ji是否成立,若成立进入步骤4-6,否则返回步骤4-2;其中ε1为功率谱迭代平滑精度控制门限;
步骤4-6:令s(k)=si(k),k=0,1,2…,n/2-1,得到平滑后功率谱s(k)。
优选的,步骤(5)中,提取原始功率谱和平滑后功率谱主瓣特征参数γp和γs,具体包括如下步骤:
步骤5-1:搜索平滑后功率谱s(k)最大值所对应的离散频率索引km
其中
步骤5-2:将s(k)减去最大值的一半得半功率差功率谱z(k):
z(k)=s(k)-s(km)/2,k=0,1,2…,n/2-1
步骤5-3:搜索z(k)峰值左侧第一过零点离散频率索引kl,搜索过程包括如下步骤:
步骤5-3-1:初始化z(k)峰值左侧第一过零点频率起始搜索索引kl=km-1;
步骤5-3-2:判断kl=0或z(kl)<0是否成立,若成立跳到步骤5-3-4,否则进入步骤5-3-3;
步骤5-3-3:令kl=kl-1,并返回步骤5-3-2;
步骤5-3-4:令kl=kl,得z(k)峰值左侧第一过零点离散频率索引kl;
步骤5-4:搜索z(k)峰值右侧第一过零点离散频率索引kr,搜索过程包括如下步骤:
步骤5-4-1:初始化z(k)峰值右侧第一过零点频率起始搜索索引kr=km+1;
步骤5-4-2:判断kr=n/2-1或z(kr)<0是否成立,若成立跳到步骤5-4-4,否则进入步骤5-4-3;
步骤5-4-3:令kr=kr+1,并返回步骤5-4-2;
步骤5-4-4:令kr=kr,得z(k)峰值右侧第一过零点离散频率索引kr。
步骤5-5:计算原始功率谱p(k)主瓣特征参数γp
步骤5-6:计算平滑后功率谱s(k)主瓣特征参数γs
优选的,步骤(6)中,生成目适应滤波参数h(k),具体包括如下步骤:
步骤6-1:初始化自适应滤波器参数
h(k)=0,k=0,1,2…,n-1
步骤6-2:分别计算窄带脉冲信号滤波器过渡带起始和终止离散频率索引tn1和tn2,以及通带起始和终止离散频率索引pn1和pn2:
tn1=max(round[(1-0.5ω1)km],0),tn2=min(round[(1+0.5ω1)km],n/2-1)
pn1=max(km-1,tn1),pn2=min(km+1,tn2)
其中round[]代表四舍五入取整运算,max(,)和min(,)分别表示取大值和取小值运算。
步骤6-3:分别计算宽带脉冲信号滤波器过渡带起始和终止离散频率索引tw1和tw2,以及通带起始和终止离散频率索引pw1和pw2:
tw1=max(round[(1-0.5ω2)km],0),tw2=min(round[(1+0.5ω2)km],n/2-1)
pw1=max(kl-1,tw1),pw2=max(kr+1,tw2)
步骤6-4:依据信号类型获得自适应滤波器过渡带起始和终止离散频率索引t1和t2,以及通带起始和终止离散频率索引p1和p2:
即判断max(γp,γs)/min(γp,γs)>η是否成立,若成立判为窄带脉冲信号并令:
t1=tn1,t2=tn2,p1=pn1,p2=pn2
否则判为宽带脉冲信号并令:
t1=tw1,t2=tw2,p1=pw1,p2=pw2
步骤6-5:对自适应滤波器系数h(k)的过渡带部分进行修正,即令k满足t1≤k≤t2的部分为:
步骤6-6:对自适应滤波器系数h(k)的通带部分进行修正,即令k满足p1≤k≤p2的部分为:h(k)=1。
优选的,步骤(7)中,对脉冲信号离散傅里叶变换结果x(k)进行自适应滤波得滤波后的离散傅里叶变换y(k):
y(k)=x(k)h(k),k=0,1,2…,n-1
优选的,步骤(8)中,提取脉冲信号包络y(n),n=0,1…,n-1,具体包括如下步骤:
步骤8-1:对滤波后的离散傅里叶变换y(k)做离散傅里叶逆变换得脉冲信号滤波后的时域复数信号yc(n):
步骤8-2:取复数yc(n)的模值得到脉冲信号初始包络数据序列y0(n):
y0(n)=|yc(n)|,n=0,1…,n-1
步骤8-3:初始化包络数据迭代平滑次数q=0;
步骤8-4:令迭代次数q=q+1,并对包络数据序列的上次迭代平滑结果yq-1(n)进行平滑处理得到当次平滑结果yq(n):
步骤8-5:判断是否达到包络最大迭代平滑次数,即判断q≤i2是否成立,若成立进入步骤8-6,否则跳到步骤8-8;
步骤8-6:分别计算包络数据序列上次平滑结果yq-1(n)和当次平滑结果yq(n)与初始包络数据序列y0(n)的残差平方和ψq-1和ψq,分别为:
步骤8-7:判断|ψq-1-ψq|≤ε2ψq是否成立,若成立进入步骤8-8,否则返回步骤8-4;
步骤8-8:令y(n)=yq(n),n=0,1,2…,n-1,得到信号包络数据序列:
y(n),n=0,1,2…,n-1。
优选的,步骤(8)中,步骤8-1对滤波后的离散傅里叶变换y(k)做离散傅里叶逆变换得脉冲信号滤波后的时域复数信号yc(n)可通过快速傅里叶变换实现。
有益效果:本发明与现有的方法相比,有以下几点有益效果:
(1)本发明的包络提取方法通过对脉冲信号主瓣特征的实时提取,在非合作条件下,对脉冲信号的类型进行自动判别,并实时估计信号带宽参数,生成与脉冲信号类型与频率相匹配的自适应滤波器参数,对脉冲信号进行自适应滤波,可获得近似匹配滤波的信号处理增益,较常规hilbert变换法,仅以增加十分小的运算量为代价,大幅提高了包络提取的噪声抑制能力和对多类型脉冲信号的适应能力。
(2)本发明的包络提取方法在对滤波后的信号离散傅里叶变换结果进行离散傅里叶逆变换,并取模值获得信号初始包络数据序列后,进一步的对初始包络做了迭代平滑处理,减少了初始包络存在的毛刺,有效地改善了常规hilbert变换法由于高频分量的存在使得包络中毛刺较多的问题。
附图说明
图1为本发明方法的流程示意图;
图2为实施例1仿真脉冲信号原始功率谱;
图3为实施例1仿真脉冲信号平滑后功率谱;
图4为实施例1自适应滤波器系数;
图5为实施例1提取出的仿真信号包络;
图6为实施例2仿真脉冲信号原始功率谱;
图7为实施例2仿真脉冲信号平滑后功率谱;
图8为实施例2自适应滤波器系数;
图9为实施例2提取出的仿真信号包络。
具体实施方式
下面结合附图和实施例对本发明做进一步的说明:
如图1所示,一种脉冲信号频域自适应滤波包络提取方法,包括如下步骤:
(1)获取待处理的脉冲信号采样数据序列:从传感器接收n个采样点的实时采集数据作为待处理的数据序列x(n),n=0,1,…,n-1,或从存储器中提取包含整个脉冲信号的n个采样点数据作为待处理的数据序列x(n),n=0,1,…,n-1,所述的n为包含整个脉冲信号的采样点数,取值为2的整数次幂。
(2)根据所述数据序列x(n)计算脉冲信号功率谱p(k):对所述数据序列x(n)做快速傅里叶变换,计算得到数据序列的离散傅里叶变换x(l)和功率谱p(k),具体包括如下步骤:
步骤2-1:计算x(n)的离散傅里叶变换为:
其中i为x(l)的离散频率索引,j表示虚数单位,即
步骤2-2:依据x(l)计算x(n)的功率谱:
其中k为功率谱p(k)的离散频率索引,||代表取模值运算,|x(k)|为信号幅度谱。
在第(2)步中,步骤2-1计算x(n)的离散傅里叶变换x(l)可通过快速傅里叶变换实现,提高算法的计算效率。
(3)迭代平滑与滤波器带宽相关参数初始化:设置功率谱迭代平滑的最大迭代次数门限,精度控制门限和迭代平滑窗长分别为i1,ε1和m1,初始包络迭代平滑的最大迭代次数门限,精度控制门限和迭代平滑窗长分别为i2,ε2和m2,窄带信号带宽系数ω1和宽带信号带宽系数ω2,判断窄带脉冲和宽带脉冲主瓣特征参数门限η。其中i1和i2取值为大于2的正整数,ε1和ε2取值为小于1的正数,功率谱迭代平滑窗长为满足3≤m1≤n/2-2的奇数,包络迭代平滑窗长为满足5≤m2≤49的奇数,η为满足1.5≤η≤2.5的任意数,ω1为满足0≤ω1≤0.05的任意数,ω2为满足0.05<ω2≤0.1的任意数。优选的i1取值为10,ε1取值为0.01,m1取值为5,i2取值为5,ε2取值为0.02,m2取值为11,ω1的取值为0.05,ω2的取值为0.1,η取值为2.0。
(4)对p(k)进行迭代平滑处理得到平滑功率谱s(k):对功率谱p(k)进行迭代平滑处理得到平滑功率谱s(k),具体包括如下步骤:
步骤4-1:初始化功率谱迭代平滑次数i=0,功率谱迭代平滑结果s0(k)为:
s0(k)=p(k),k=0,1,2…,n/2-1
步骤4-2:令迭代次数i=i+1,并对功率谱的上次迭代平滑结果si-1(k)进行平滑处理得到当次平滑结果si(k):
步骤4-3:判断是否达到功率谱最大迭代平滑次数,即判断i≤i1是否成立,若成立进入步骤4-4,否则跳到步骤4-6;
步骤4-4:分别计算功率谱上次平滑结果si-1(k)和当次平滑结果si(k)与原始功率谱p(k)的残差平方和ji-1和ji,分别为:
步骤4-5:判断|ji-1-ji|≤ε1ji是否成立,若成立进入步骤4-6,否则返回步骤4-2;其中ε1为功率谱迭代平滑精度控制门限;
步骤4-6:令s(k)=si(k),k=0,1,2…,n/2-1,得到平滑后功率谱s(k)。
(5)提取原始功率谱和平滑后功率谱主瓣特征参数γp和γs,具体包括如下步骤:
步骤5-1:搜索平滑后功率谱s(k)最大值所对应的离散频率索引km
其中
步骤5-2:将s(k)减去最大值的一半得半功率差功率谱z(k):
z(k)=s(k)-s(km)/2,k=0,1,2…,n/2-1
步骤5-3:搜索z(k)峰值左侧第一过零点离散频率索引kl,搜索过程包括如下步骤:
步骤5-3-1:初始化z(k)峰值左侧第一过零点频率起始搜索索引kl=km-1;
步骤5-3-2:判断kl=0或z(kl)<0是否成立,若成立跳到步骤5-3-4,否则进入步骤5-3-3;
步骤5-3-3:令kl=kl-1,并返回步骤5-3-2;
步骤5-3-4:令kl=kl,得z(k)峰值左侧第一过零点离散频率索引kl。
步骤5-4:搜索z(k)峰值右侧第一过零点离散频率索引kr,搜索过程包括如下步骤:
步骤5-4-1:初始化z(k)峰值右侧第一过零点频率起始搜索索引kr=km+1;
步骤5-4-2:判断kr=n/2-1或z(kr)<0是否成立,若成立跳到步骤5-4-4,否则进入步骤5-4-3;
步骤5-4-3:令kr=kr+1,并返回步骤5-4-2;
步骤5-4-4:令kr=kr,得z(k)峰值右侧第一过零点离散频率索引kr。
步骤5-5:计算原始功率谱p(k)主瓣特征参数γp:
步骤5-6:计算平滑后功率谱s(k)主瓣特征参数γs:
(6)生成目适应滤波参数h(k),具体包括如下步骤:
步骤6-1:初始化自适应滤波器参数
h(k)=0,k=0,1,2…,n-1
步骤6-2:分别计算窄带脉冲信号滤波器过渡带起始和终止离散频率索引tn1和tn2,以及通带起始和终止离散频率索引pn1和pn2:
tn1=max(round[(1-0.5ω1)km],0),tn2=min(round[(1+0.5ω1)km],n/2-1)
pn1=max(km-1,tn1),pn2=min(km+1,tn2)
其中round[]代表四舍五入取整运算,max(,)和min(,)分别表示取大值和取小值运算。
步骤6-3:分别计算宽带脉冲信号滤波器过渡带起始和终止离散频率索引tw1和tw2,以及通带起始和终止离散频率索引pw1和pw2:
tw1=max(round[(1-0.5ω2)km],0),tw2=min(round[(1+0.5ω2)km],n/2-1)
pw1=max(kl-1,tw1),pw2=max(kr+1,tw2)
步骤6-4:依据信号类型获得自适应滤波器过渡带起始和终止离散频率索引t1和t2,以及通带起始和终止离散频率索引p1和p2:
即判断max(γp,γs)/min(γp,γs)>η是否成立,若成立判为窄带脉冲信号并令:
t1=tn1,t2=tn2,p1=pn1,p2=pn2
否则判为宽带脉冲信号并令:
t1=tw1,t2=tw2,p1=pw1,p2=pw2
步骤6-5:对自适应滤波器系数h(k)的过渡带部分进行修正,即令k满足t1≤k≤t2的部分为:
步骤6-6:对自适应滤波器系数h(k)的通带部分进行修正,即令k满足p1≤k≤p2的部分为:h(k)=1;
(7)依据所述滤波器参数h(k)对脉冲信号进行自适应滤波,即对脉冲信号离散傅里叶变换结果x(k)进行自适应滤波得滤波后的离散傅里叶变换y(k):
y(k)=x(k)h(k),k=0,1,2…,n-1
(8)提取脉冲信号包络y(n),具体包括如下步骤:
步骤8-1:对滤波后的离散傅里叶变换y(k)做离散傅里叶逆变换得脉冲信号滤波后的时域复数信号yc(n):
步骤8-2:取复数yc(n)的模值得到脉冲信号初始包络数据序列y0(n):
y0(n)=|yc(n)|,n=0,1…,n-1
步骤8-3:初始化包络数据迭代平滑次数q=0;
步骤8-4:令迭代次数q=q+1,并对包络数据序列的上次迭代平滑结果yq-1(n)进行平滑处理得到当次平滑结果yq(n)
步骤8-5:判断是否达到包络最大迭代平滑次数即判断q≤i2是否成立,若成立进入步骤8-6,否则跳到步骤8-8;
步骤8-6:分别计算包络数据序列上次平滑结果yq-1(n)和当次平滑结果yq(n)与初始包络数据序列y0(n)的残差平方和ψq-1和ψq,分别为:
步骤8-7:判断|ψq-1-ψq|≤ε2ψq是否成立,若成立进入步骤8-8,否则返回步骤8-4;
步骤8-8:令y(n)=yq(n),n=0,1,2…,n-1,得到信号包络数据序列:
y(n),n=0,1,2…,n-1。
其中,步骤8-1可通过快速傅里叶变换实现,提高算法的计算效率。
本发明的实施例中,仿真矩形包络脉冲信号模型为:
其中a为信号幅度,
以采样频率fs对上述脉冲信号进行离散采样可得到脉冲信号采样数据序列:
其中,n0=int(τ0fs),m=int(τfs),n=int(tfs)。
实施例1:
仿真信号参数分别设置为:信号幅度a=1,初始相位
依据第(2)步,计算所述数据序列x(n)的功率谱p(k),如图2所示。
在第(3)步中,设置功率谱迭代平滑最大迭代次数门限i1=10,精度控制门限ε1=0.01,迭代平滑窗长m1=5,包络迭代平滑最大迭代次数门限i2=5,精度控制门限ε2=0.02,迭代平滑窗长m2=11,窄带信号带宽系数ω1=0.05,宽带信号带宽系数ω2=0.1,判断窄带脉冲和宽带脉冲主瓣特征参数门限η=2。
依据第(4)步,对p(k)进行迭代平滑处理得到平滑功率谱s(k),如图3所示,对比图2和图3可以看出:经过迭代平滑处理后,窄带脉冲信号的功率谱主瓣明显展宽,整个频带谱形更为平滑。
依据第(5)步,提取原始功率谱和平滑后功率谱主瓣特征参数γp和γs分别为
γp=45.6772,γs=1.4108
依据第(6)步,生成自适应滤波器参数h(k)如图4所示。
依据第(7)和第(8)步,提取出脉冲信号包络如图5所示,由图5可以看出,提取出的窄带脉冲信号包络前后沿特征明显,可为进一步提取窄带脉冲信号到达时间和结束时间提供基础。
实施例2:
仿真信号参数分别设置为:信号幅度a=2,初始相位
依据第(2)步,计算所述数据序列x(n)的功率谱p(k),如图6所示。
在第(3)步中,设置功率谱迭代平滑最大迭代次数门限i1=10,精度控制门限ε1=0.01,迭代平滑窗长m1=7,包络迭代平滑最大迭代次数门限i2=10,精度控制门限ε2=0.02,迭代平滑窗长m2=13,窄带信号带宽系数ω1=0.05,宽带信号带宽系数ω2=0.1,判断窄带脉冲和宽带脉冲主瓣特征参数门限η=2。
依据第(4)步,对p(k)进行迭代平滑处理得到平滑功率谱s(k),如图7所示,对比图6和图7可以看出:经过迭代平滑处理后,宽带脉冲信号的功率谱主瓣变化不明显,整个频带谱形更为平滑。
依据第(5)步,提取原始功率谱和平滑后功率谱主瓣特征参数γp和γs分别为:
γp=1.1081,γs=1.1974
依据第(6)步,生成自适应滤波器参数h(k)如图8所示。
依据第(7)和第(8)步,提取出脉冲信号包络如图9所示,由图9可以看出,提取出的宽带脉冲信号包络前后沿特征明显,可为进一步提取宽带脉冲信号到达时间和结束时间提供基础。