一种基于峰值包络谱的傅里叶分解信号降噪方法

1.本发明涉及一种基于峰值包络谱的傅里叶分解信号降噪方法，属于数字信号降噪方法技术领域。


背景技术：

2.傅里叶分解(fdm)是一种非线性非平稳时间序列分析的信号分解方法。fdm已被应用到生物医学信号分析、机械故障诊断等领域。一般利用fdm将信号分解为一系列傅里叶本征带函数(fibfs)，即分量信号，然后通过对分量信号进行筛选并重构达到降噪的效果，或者提取筛选分量信号的特征进行进一步的分析。然而，fdm获得单个分量的条件过于容易满足，获得的分量信号未必满足瞬时频率具有物理意义的条件。对于所含频率组分较多的信号，fdm将分解产生大量无意义的fibfs，干扰对有用分量信号的分析。另外，由于fdm是根据频带所对应的信号的相位是否符合单调递增的条件来获得分量，计算非常耗时。
3.针对上述问题，fdm的一些改进方法已经被提出，如自适应经验傅里叶分解(aefd)、自适应功率谱傅里叶分解方法(apsfdm)等，这些方法较好的解决了fdm的过分解问题，但都还存在着一定的不足，如aefd边界集的设置方法还需近一步研究，apsfdm中功率谱的最优阶数选择还存在问题。


技术实现要素：

4.本发明为了解决现有的傅里叶分解方法确定单个分量的条件过于容易满足，导致对信号过度分解而影响信号降噪效果，以及傅里叶分解方法计算耗时的问题，提供一种基于峰值包络谱的傅里叶分解信号降噪方法。
5.本发明的技术方案：
6.一种基于峰值包络谱的傅里叶分解信号降噪方法，步骤如下：
7.s1，对离散时间序列长度为n的信号x(k)进行快速傅里叶变换得到x(k)的复系数x(n)，接着由复系数x(n)获得信号的频谱a(n)；
8.s2，对信号的频谱a(n)进行上峰值包络处理；
9.s3，使用改进的locmaxmin频谱分割方法检测频谱分割边界；
10.s4，利用检测得到的分割边界将频谱a(n)分割为多个频带，将各个频带进行傅里叶逆变换得到降噪处理后的时域信号。
11.进一步限定，s1中离散时间序列长度为n的信号x(k)的复系数x(n)和频谱a(n)分别由下式(1)和(2)获得；
[0012][0013][0014]
式中，n是信号的长度，k是时域离散值序号，k＝0,1,2,
···
,n-1，n是频域离散
值的序号。
[0015]
进一步限定，s2的操作过程为：
[0016]
首先，减少频谱中的局部极大值的数目：设置最小峰值距离l来忽略频谱中可能存在大局部峰值附近的较小的峰值，最小峰值距离l的公式如下：
[0017]
l＝round(α
×
β
×
n)
ꢀꢀꢀ
(3)
[0018]
式中，α和β为用于调节最小峰值距离的系数，n是信号的长度；
[0019]
然后，设置α和β的初始值，在频谱幅值序列中找到至少由l个幅值序列点数分隔的局部极大值，获得的局部极大值的数目为k；
[0020]
最后，利用三次样条曲线对频谱幅值序列的两个端点以及k个局部最大值进行插值，使总的插值点数与频谱幅值序列点数相同，得到峰值包络谱。
[0021]
更进一步限定，α的初始值为0.01。
[0022]
更进一步限定，β的初始值为1.5。
[0023]
进一步限定，若k《10，以0.1的步长减小β，接着以新的l在频谱幅值序列中再次寻找符合要求的局部极大值，直至k不小于10。
[0024]
进一步限定，s3具体步骤为：
[0025]
首先，将峰值包络谱的前10个最大的局部最大值对应的频率按升序索引，将其表示为[ω1,ω2,
···
,ω
10
]，另外，定义ω0＝0；
[0026]
然后，获得初始分割边界集{ns}，其中s∈[1,10]，ns为频谱上峰值包络谱中两个相邻局部最大值中的最小值对应的频率的位次，如下式(4)所示：
[0027]ns
＝argminexs(n)
ꢀꢀꢀ
(4)
[0028]
式中，exs(n)表示峰值包络谱在频率区间[ω
s-1
,ωs]内的幅值序列，argmin(
·
)表示取最小幅值的频率在[0,fs/2)的位次，fs为信号的采样频率；
[0029]
再然后，将分量个数设置为m，m《10，将n1,n2,
···
,n
10
随机排列得到{n's}，取{n's}的前m-1个并按照升序排列得到[b1,b2,
···
,b
m-1
]；定义b0＝0，bm＝(n-1)/2，则分割边界集合为：{b
l
},l∈[0,m]，其中，{b
l
}中分割边界按照频率的位次升序排列。利用{b
l
}对信号频谱进行分割得到m个频带，计算m个频带的能量并表示为[e1,e2,
···em
],计算它们的方差,用v表示；
[0030]
最后，将{ns}进行1000次随机排列，与之对应获得1000个分割边界集合{b
l
}，相应的得到频带能量方差的集合{vm},m∈[1,1000]，最小的v所对应的分割边界集{b
l
}为最终的频谱分割边界集。
[0031]
进一步限定，s4具体步骤为：
[0032]
由最终的分割边界集{b
l
}分割频谱得到的每个频带区间为λ
l
＝δf
×
[b
l-1
,b
l
),l∈[1,m]，其中δf＝fs/n为频率分辨率，将各个频段进行傅里叶逆变换得到降噪处理后的时域信号x
l
(k)，如下式(5)所示：
[0033][0034]
式中，n是信号的长度，k是时域离散值序号，n是频域离散值的序号，l＝1,2,
···
,m。
[0035]
本发明具有以下有益效果：
[0036]
本发明对信号的频谱进行峰值包络得到频谱的峰值包络谱，使用基于频带能量方差最小化的locmaxmin频谱分割方法在峰值包络谱上自适应确定频谱分割边界，再将分割边界之间的频带重构为多个单分量信号，提出了一种基于峰值包络谱的傅里叶分解信号降噪方法，解决了傅里叶分解的过分解以及计算耗时的问题。本发明提出的峰值包络谱的傅里叶分解信号降噪方法不会产生过窄或过宽的分量，能够有效将故障信息丰富的频带与噪声频带分来，能取得良好的降噪效果。
附图说明
[0037]
图1为基于峰值包络谱的傅里叶分解信号降噪方法的流程图；
[0038]
图2为实施例1中轴承外圈故障仿真信号各组成分量及叠加信号的时域图；
[0039]
图3为实施例1中利用pesfdm分解仿真信号得到的分量信号的时域图；
[0040]
图4为实施例1中利用pesfdm分解仿真信号得到的分量信号的包络谱图；
[0041]
图5为对比例1中利用fdm分解仿真信号得到的分量信号的时域图；
[0042]
图6为对比例1中利用fdm分解仿真信号得到的分量信号的包络谱图；
[0043]
图7为对比例2中利用emd分解仿真信号得到的分量信号的时域图；
[0044]
图8为对比例2中利用emd分解仿真信号得到的分量信号的包络谱图；
[0045]
图9为对比例2中emd和pesfdm各自信噪比最高的分量的信噪比。
具体实施方式
[0046]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
[0047]
下述实施例中所使用的实验方法如无特殊说明均为常规方法。所用材料、试剂、方法和仪器，未经特殊说明，均为本领域常规材料、试剂、方法和仪器，本领域技术人员均可通过商业渠道获得。
[0048]
具体实施方式一：
[0049]
如图1所示，本实施例的具体步骤如下：
[0050]
步骤1，利用式(1)对离散时间序列长度为n的待分解信号x(k)进行快速傅里叶变换(fft),得到其复系数x(n)，其中fft的变换长度为信号本身的长度，
[0051][0052]
然后，利用式(2)得到信号的频谱，
[0053][0054]
其中，k是时域离散值序号，k＝0,1,2,
···
,n-1，n是频域离散值的序号。
[0055]
步骤2，对信号频谱进行上峰值包络处理，获得峰值包络谱，具体过程为：
[0056]
首先，通过式(3)设置最小峰值距离系数l：
[0057]
l＝round(α
×
β
×
n)
ꢀꢀꢀ
(3)
[0058]
将α设置为0.01，β初始值设置为1.5，在α＝0.01，β＝1.5的初始条件下，在频谱上
找到至少由l个频谱幅值序列点数分开的局部极大值，得到得局部极大值数目为k,若k小于10，以0.1为步长减小β，再以新的l在频谱上寻找彼此间隔至少为l的局部极大值，直至k不小于10为止。
[0059]
然后，利用三次样条曲线对频谱幅值序列的两个端点以及k个局部最大值进行插值，使总的插值点数与频谱幅值序列点数相同，得到峰值包络谱。
[0060]
步骤3，在峰值包络谱上检测频谱分割边界，具体步骤为：
[0061]
首先，将步骤2中峰值包络谱的前10个最大的局部最大值对应的频率按升序索引，将其表示为[ω1,ω2,
···
,ω
10
]，另外，定义ω0＝0；
[0062]
然后利用式(4)获得初始分割边界集{ns}；
[0063]ns
＝argminexs(n)
ꢀꢀꢀ
(4)
[0064]
式中，s∈[1,10]，ns为频谱上峰值包络谱中两个相邻局部最大值中的最小值对应的频率的在[0,fs/2)的位次，exs(n)表示峰值包络谱在频率区间[ω
s-1
,ωs]内的幅值序列，argmin(
·
)表示取最小幅值的频率在[0,fs/2)的位次，fs为信号的采样频率。
[0065]
再然后，设置分量个数为m，将n1,n2,
···
,n
10
随机排列得到{n's}，取{n's}的前m-1个并按照升序排列得到[b1,b2,
···
,b
m-1
]。此外，定义b0＝0，bm＝(n-1)/2，则分割边界集合为{b
l
},l∈[0,m]，{b
l
}中分割边界按照频率的位次升序排列；
[0066]
再然后，利用{b
l
}对信号频谱进行分割得到m个频带，计算m个频带的能量并表示为[e1,e2,
···em
],计算它们的方差,用v表示；
[0067]
最后，将{ns}进行1000次随机排列，与之对应获得1000个分割边界集合{b
l
}，相应的得到频带能量方差的集合{vm},m∈[1,1000]，选择最小的v所对应的分割边界集{b
l
}作为最终的频谱分割边界。
[0068]
步骤4，利用得到的分割边界集{b
l
}对频谱进行分割得到m个频带，每个频带区间可表示为λ
l
＝δf
×
[b
l-1
,b
l
),l∈[1,m]，其中δf＝fs/n为频率分辨率。利用式(5)将各个频带进行傅里叶逆变换得到降噪处理后的时域信号。
[0069][0070]
其中，x
l
(k)即为获得的时域分量信号。
[0071]
实施例1：
[0072]
为验证本发明的有效性，构建了信号长度为1s，采样频率为12khz的轴承外圈故障仿真信号，并利用本发明提出的pesfdm对其进行分解。
[0073]
轴承外圈故障仿真信号由四部分组成，分别为轴承外圈缺陷周期性冲击信号b(t)、随机意外冲击信号d(t)、模拟受到其它机械部件的离散谐波干扰信号h(t)、模拟轴承工作背景噪声的高斯白噪声信号n(t)。b(t)、d(t)、h(t)的具体形式分别如公式(6-8)所示，轴承外圈故障仿真信号模型如公式(9)。
[0074][0075]
式(6)中，j是外圈故障脉冲产生的次数，被设置为100，故障特征频率fm为100hz，共振频率f1和阻尼系数β1分别被设置为1700hz、1000。故障脉冲幅值aj是在[0.8,1]的均匀分布中随机选择的。
[0076][0077]
式(7)中，随机意外冲击的个数m1被设置为5，幅值dj是在均值为1.5、标准差为2的高斯分布里随机选取的，随机冲击产生的时刻ti由(0,1)区间内的随机数生成。另外，阻尼系数β2为800，共振频率f2为2000hz。
[0078][0079]
式(8)中，p1，p2，p3分别表示一次、二次和三次谐波的幅值，分别被设置为0.025，0.025，0.05。
[0080]
s(t)＝b(t) d(t) h(t) n(t)
ꢀꢀꢀ
(9)
[0081]
式(9)中，高斯白噪声n(t)的均值为0，标准差为0.3。
[0082]
轴承外圈故障仿真信号的各分量信号以及叠加信号的时域波形如图2所示，叠加信号中轴承外圈缺陷周期性冲击信号严重受到其它分量干扰，观察不到明显的周期性冲击。
[0083]
利用pesfdm对该信号进行分解，其中分量个数设置为6，α＝0.01，β自动调整为1，分解获得的6个分量信号的时域波形图以及包络谱图分别如图3、图4所示。
[0084]
其中，分量信号的包络谱是由式(10-12)计算；
[0085]
q(k)＝x(k) i
·
x
′
(k)k＝0，1，2，
…
，n-1
ꢀꢀꢀ
(10)
[0086]
式中，对信号x(k)进行hilber变换得到x'(k)，进而得到x(k)的解析信号q(k)；
[0087][0088]
式中，|q(k)|是x(k)的包络；
[0089][0090]
式中，对|q(k)|进行fft得到包络信号的复数q(n)，则|q(n)|即为x(k)的包络谱。
[0091]
在图3中，第二个分量周期性冲击明显且干扰噪声少。在图4中，第二个分量可以清楚的观察到轴承外圈故障特征频率的一倍频到九倍频的波峰。对于该仿真信号，本发明提出的方法起到了良好的降噪效果。另外，pesfdm计算速度快，对于该信号，用时不足1s。
[0092]
对比例1：
[0093]
为了证明实施例1提供的基于峰值包络谱的傅里叶分解信号降噪方法相比于fdm降噪效果的有效性和优越性，采用fdm对相同仿真信号进行分解。
[0094]
对于该轴承外圈故障仿真信号，fdm过分解问题严重，从高频到低频和低频到高频的方式都获得52个分量信号。由于分量个数过多，图5、图6分别只展示了从高频到低频方式获得的前12个分量的时域波形图和包络谱。
[0095]
由于fdm存在过分解问题，产生了很多瞬时频率没有意义的分量，将故障信息丰富的频带分解为几个分量，这些分量的包络谱观察不到故障特征频率的波峰或者只能观察到极少的故障特征频率的波峰。另外，对于该信号，fdm用时166.49s，相比于本发明提出的pesfdm计算非常耗时。
[0096]
通过对pesfdm和fdm对轴承外圈故障仿真信号分解得到的分量的时域波形以及包
络谱的对比，pesfdm将故障信息丰富的频带与噪声频带分割开来，取得了更好的信号降噪效果。另外，本发明提出的基于pesfdm的信号降噪方法不存在fdm的过分解和计算耗时的问题。
[0097]
对比例2：
[0098]
为了证明基于峰值包络谱的傅里叶分解信号降噪方法降噪效果的有效性和优越性，采用在数字信号降噪方法技术领域被广泛应用的经验模态分解(emd)对相同仿真信号进行分解。
[0099]
对于该仿真信号，emd获得了10个分量，分量的时域波形图和包络谱分别如图7、图8所示。图7中，emd的前5个分量包含仿真信号的主要成分，后5个分量都是虚假分量。另外，emd的前三个分量中包含的噪声多，故障冲击不明显。图8中，第一个分量可以观察到故障特征频率的前五倍频的波峰，在第二个分量中，可以看到故障特征频率的一倍频到六倍频的波峰，但这些波峰的幅值均低于pesfdm的第二个分量对应波峰的幅值。
[0100]
然后，利用式(13)、(14)计算包络信号的信噪比，来对降噪效果进行定量评价。
[0101][0102][0103]
其中，p
signal
是信号功率，p
noise
是噪声功率，p[
·
]表示用于计算信噪比的分量的包络谱的幅值，f为轴承理论故障特征频率，δf是频率分辨率，fds为包络谱中故障特征频率的一倍频所在的幅值序列的位次，h表示故障特征频率第h个倍频。
[0104]
由于emd后面五个分量不包含表示故障特征频率的波峰，因此只对pesfdm和emd各自前5个分量的信噪比进行计算并挑选各自信噪比最高的分量进行对比。emd和pesfdm各自第二个分量的信噪比最高，这两个分量的信噪比如图9所示。在图9中，除了pesfdm的故障特征频率的一倍频的信噪比略低于emd的一倍频的信噪比外，其余故障特征频率的信噪比，pesfdm均高于emd。因此，本发明提出的基于pesfdm的信号降噪方法取得了比emd更好的降噪效果。
[0105]
以上所述仅为本发明的较佳实施例而已，鉴于本发明所属领域的技术人员可以对上述实施方式进行适当的变更和修改，因此，本发明并不局限于上面所述的具体实施方式，对本发明的一些修改和变更也应当落入本发明的权利要求的保护范围之内。