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

本发明专利技术公开了一种基于FFT累积量的循环平稳相关熵谱密度计算方法，包括以下步骤：给定信号x、核长σ、信号采样频率f

【技术实现步骤摘要】
基于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进行分段，得到各分段信号xxi，xxi的采样点数为Nw；步骤3，计算各分段信号xxi的中心化相关熵Exix；步骤4，将所有Exix作为列向量，创建中心化相关熵矩阵X；步骤5，对矩阵X进行加窗处理，得到矩阵XW；步骤6，按列计算矩阵XW的FFT，得到矩阵XF1及其共轭矩阵CXF1；步骤7，将计算矩阵XF1的每一行向量的元素与CXF1的每一行向量的元素相乘，得到矩阵XM；步骤8，按行计算矩阵XM的FFT并转置，得到循环平稳相关熵谱密度S。在上述技术方案中，所述步骤1中，在上述技术方案中，所述步骤2中，xxi计算公式为在上述技术方案中，所述步骤3中，计算的公式为：其中κ(·)为高斯核函数，高斯核函数κ[xxi(j),xxi(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进行分段，得到各分段信号xxi，xxi的采样点数为Nw，xxi计算公式为步骤3，计算各分段信号xxi的中心化相关熵计算的公式为：其中κ(·)为高斯核函数，高斯核函数κ[xxi(j),xxi(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和矩本文档来自技高网...

【技术保护点】
1.一种基于FFT累积量的循环平稳相关熵谱密度计算方法，其特征在于，包括以下步骤：/n步骤1，给定信号x、核长σ、信号采样频率f

【技术特征摘要】
1.一种基于FFT累积量的循环平稳相关熵谱密度计算方法，其特征在于，包括以下步骤：
步骤1，给定信号x、核长σ、信号采样频率fs、循环频率分辨率Δα、谱频率分辨率Δf，分别计算信号x分段序列的采样点长度Nw、相邻两个数据分段平移的采样点数R和将信号x进行平均分段的数量K；
步骤2，将信号x进行分段，得到各分段信号xxi，xxi的采样点数为Nw；
步骤3，计算各分段信号xxi的中心化相关熵
步骤4，将所有作为列向量，创建中心化相关熵矩阵X；
步骤5，对矩阵X进行加窗处理，得到矩阵XW；
步骤6，按列计算矩阵XW的FFT，得到矩阵XF1及其共轭矩阵CXF1；
步骤7，将计算矩阵XF1的每一行向量的元素与CXF1的每一行向量的元素相乘，得到矩阵XM；
步骤8，按行计算矩阵XM的FFT并转置，得到循环平稳相关熵谱密度S。


2.如权利要求1所述的循环平稳相关熵谱密度计算方法，其特征在于，所述步骤1中，


3.如权利要求1所述的循环平稳相关熵谱密度计算方法，其特征在于，所述步骤2中，xxi计算公式为





4.如权利要求1所述的循环平稳相关熵谱密度计算方法，其特征在于，所述步骤3中，计算的公式为：



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





5.如权利要求1...

【专利技术属性】
技术研发人员：李辉，张宇平，
申请(专利权)人：天津职业技术师范大学中国职业培训指导教师进修中心，
类型：发明
国别省市：天津;12