基于FFT累积量的循环平稳相关熵谱密度计算方法与流程

本发明涉及循环平稳信号技术领域，特别是涉及一种基于fft累积量的循环平稳相关熵谱密度计算方法。

背景技术：

循环平稳信号处理在工程应用中具有非常重要的现实意义，已在通信、机电设备故障诊断等领域得到了广泛的应用。目前，循环平稳信号处理主要采用基于信号二阶统计量的谱相关密度方法和循环平稳相关熵谱密度方法。谱相关密度估计常用的方法有循环相关图和循环周期图。为了提高谱相关密度的估计效率，学者们提出了多种估计算法。在诸多的谱相关密度的估计算法中，fft累积量法(fftaccumulationmethod，简称fam)，将谱相关平面分成多个循环谱分析区域，对不同区域的循环频率α利用fft进行细化，大大提高了计算效率。经过几十年的发展，基于谱相关密度方法的循环平稳信号处理技术，已日臻成熟。

而基于循环平稳相关熵谱的信号处理方法是最近几年才提出的，循环平稳相关熵谱密度的估计方法还没有成熟的高效算法，都采用巴西学者fontes提出的算法(aluiso.i.r.fontes,etal,cyclostationarycorrentropy:definitionandapplication,expertsystemswithapplication.2017，69(1)：110-117)，该算法采用循环周期图检测(cyclicperiodogramdetection，简称cpd)方法计算循环平稳相关熵谱密度，cpd算法采用多重嵌套循环语句计算循环平稳相关熵谱密度，存在计算效率低且计算方法复杂、占用计算机内存多的缺点，严重制约了循环平稳相关熵谱密度在信号处理领域的推广应用，因此，亟需探寻开发高效率的循环平稳相关熵谱密度估计方法。

技术实现要素：

本发明的目的是针对现有技术中存在的cpd算法计算效率低的缺点，而提供一种基于fft累积量的循环平稳相关熵谱密度计算方法。

为实现本发明的目的所采用的技术方案是：

一种基于fft累积量的循环平稳相关熵谱密度计算方法，包括以下步骤：

步骤1，给定信号x、核长σ、信号采样频率fs、循环频率分辨率δα、谱频率分辨率δf，分别计算信号x分段序列的采样点长度nw、相邻两个数据分段平移的采样点数r和将信号x进行平均分段的数量k；

步骤2，将信号x进行分段，得到各分段信号xxⁱ，xxⁱ的采样点数为nw；

步骤3，计算各分段信号xxⁱ的中心化相关熵exix；

步骤4，将所有exix作为列向量，创建中心化相关熵矩阵x；

步骤5，对矩阵x进行加窗处理，得到矩阵xw；

步骤6，按列计算矩阵xw的fft，得到矩阵xf1及其共轭矩阵cxf1；

步骤7，将计算矩阵xf1的每一行向量的元素与cxf1的每一行向量的元素相乘，得到矩阵xm；

步骤8，按行计算矩阵xm的fft并转置，得到循环平稳相关熵谱密度s。

在上述技术方案中，所述步骤1中，

在上述技术方案中，所述步骤2中，xxⁱ计算公式为

在上述技术方案中，所述步骤3中，计算的公式为：

其中κ(·)为高斯核函数，高斯核函数κ[xxⁱ(j),xxⁱ(j-τ)]的计算公式为：

其中σ为高斯核函数的核长，即步骤1中的核长σ。

在上述技术方案中，所述步骤4中，矩阵x的行数为nw、列数为k。

在上述技术方案中，所述步骤5中，xw＝diag[w(nw)]*x，其中w(nw)为长度为nw对称窗函数；diag(·)为对角矩阵算子，diag[w(nw)]是由长度为nw对称窗函数组成的对角矩阵；*为矩阵乘算子，矩阵xw的行数为nw、列数为k。

在上述技术方案中，所述对称窗函数为hamming窗函数或hanning窗函数。

在上述技术方案中，所述步骤6中，

其中：fft(·)为快速傅里叶变换算子；conj(·)为共轭算子；矩阵xf1和矩阵cxf1的行数为nw、列数为k。

在上述技术方案中，所述步骤7中，

xm((i-1)*nw+j,:)＝xf1(i,:).*cxf1(j,:)i,j＝1,2,3,…,nw

其中矩阵xm的行数为nw×nw、列数为k。

在上述技术方案中，所述步骤8中，

s＝[fft(xm,k,2)]^t

其中fft(·)为快速傅里叶变换算子；[·]^t为矩阵转置算子；矩阵s的行数为k、列数为nw×nw。

与现有技术相比，本发明的有益效果是：

本发明借鉴基于fft累积量的谱相关密度估计方法，提出了一种基于fft累积量的循环平稳相关熵谱密度计算方法，不仅提高了计算效率，而且提高了频谱分辨率。

附图说明

图1为本发明所述方法的流程图。

图2为实施例仿真调幅信号x(t)的时域波形。

图3为实施例仿真调幅信号x(t)的傅里叶变换(fft)。

图4为实施例仿真调幅信号x(t)，根据本发明计算的循环平稳相关熵谱密度平面轮廓图(σ＝0.2)。

图5为实施例仿真调幅信号x(t)，根据本发明计算的循环平稳相关熵谱密度三维立体图(σ＝0.2)。

图6为实施例仿真调幅信号x(t)，根据cpd算法计算的循环平稳相关熵谱密度平面轮廓图(σ＝0.2)。

图7为实施例仿真调幅信号x(t)，根据cpd算法计算的循环平稳相关熵谱密度三维立体图(σ＝0.2)。

具体实施方式

以下结合具体实施例对本发明作进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

实施例1

一种基于fft累积量的循环平稳相关熵谱密度计算方法，包括如下步骤：

步骤1，给定信号x、核长σ、信号采样频率fs、循环频率分辨率δα、谱频率分辨率δf；

根据已知参数，分别计算信号x分段序列的采样点长度nw、相邻两个数据分段平移的采样点数r和将信号x进行平均分段的数量k。

步骤2，将信号x进行分段，得到各分段信号xxⁱ，xxⁱ的采样点数为nw，xxⁱ计算公式为

步骤3，计算各分段信号xxⁱ的中心化相关熵计算的公式为：

其中κ(·)为高斯核函数，高斯核函数κ[xxⁱ(j),xxⁱ(j-τ)]的计算公式为：

其中σ为高斯核函数的核长，即步骤1中的核长σ。

步骤4，将所有作为列向量，创建中心化相关熵矩阵x，即

矩阵x的行数为nw、列数为k。

步骤5，对矩阵x进行加窗处理，得到矩阵xw；

xw＝diag[w(nw)]*x

其中w(nw)为长度为nw对称窗函数，例如hamming窗函数、hanning窗函数等；diag(·)为对角矩阵算子，diag[w(nw)]是由长度为nw对称窗函数组成的对角矩阵；*为矩阵乘算子。矩阵xw的行数为nw、列数为k。

步骤6，按列计算矩阵xw的fft，得到矩阵xf1及其共轭矩阵cxf1；

其中：fft(·)为快速傅里叶变换算子；conj(·)为共轭算子；矩阵xf1和矩阵cxf1的行数为nw、列数为k。

步骤7，将计算矩阵xf1的每一行向量的元素与cxf1的每一行向量的元素相乘，得到矩阵xm；

xm((i-1)*nw+j,:)＝xf1(i,:).*cxf1(j,:)i,j＝1,2,3,…,nw

其中矩阵xm的行数为nw×nw、列数为k

步骤8，按行计算矩阵xm的fft并转置，得到循环平稳相关熵谱密度s。

s＝[fft(xm,k,2)]^t

其中fft(·)为快速傅里叶变换算子；[·]^t为矩阵转置算子；矩阵s的行数为k、列数为nw×nw。

实施例2

对实施例1的方法进行仿真验证，采用的实施例仿真调幅信号为x(t)＝a(1+cos(2πf0t))cos(2πfct)，采样时间t＝2s，采样频率为fs＝1024hz，信号x(t)长度n＝2048，调幅频率f0＝30hz，载波频率fc＝250hz，幅值a＝1。循环频率分辨率δα＝2hz，谱频率分辨率δf＝16hz。实施例仿真调幅信号x(t)的时域波形，如图2所示；实施例仿真调幅信号x(t)的频域波形，如图3所示。

图4和图5分别为实施例2仿真调幅信号x(t)，当高斯核函数的核长σ＝0.2时，根据本发明计算的循环平稳相关熵谱密度的平面轮廓图及其三维立体图。采用计算机参数为：inteli7-2600，双核cpu3.4ghz，内存3gb，windowxp操作系统，使用的软件为matlabrelease12.1，计算时间为0.641s。

对比例1

为对比本发明提出的循环平稳相关熵谱密度估计方法与cpd方法(aluiso.i.r.fontes,etal,cyclostationarycorrentropy:definitionandapplication,expertsystemswithapplication.2017，69(1)：110-117)的计算效率，当高斯核函数的核长σ＝0.2时，图6和图7分别给出了实施例采用cpd方法计算的循环平稳相关熵谱密度的平面轮廓图及其三维立体图，计算时间为334.422s。

对比图4与图6可以看出：本发明方法计算的循环平稳相关熵谱密度，在频谱图中是相互分离的谱峰，频谱分辨率高、无频谱泄漏现象；而cpd方法的计算结果谱峰之间是互连的，表明存在比较严重的频谱泄漏现象，频谱分辨率较低。因此。采用本发明提出的方法，不仅提高了频谱分辨率，而且还避免了频谱泄漏现象。

以上所述仅是本发明的优选实施方式，应当指出的是，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。