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

本发明涉及信号处理
技术领域：
，特别是涉及一种基于相关图的循环平稳相关熵谱密度计算方法，适用于特定应用的数字计算或信号处理的设备或方法。
背景技术：
：循环平稳相关熵谱密度在2015年提出后，已在无线通信信号处理等领域得到广泛应用。在已发表的循环平稳相关熵谱密度学术论文中，采用循环周期图检测(cyclicperiodogramdetection，简称cpd)算法计算循环平稳相关熵谱密度，在算法上采用多重嵌套循环语句计算循环平稳相关熵谱密度，存在计算效率低且计算方法复杂、占用计算机内存多的缺点，严重制约了循环平稳相关熵谱密度在信号处理领域的推广应用。针对现有技术的不足，亟需探寻更加简单、高效的方法计算循环平稳相关熵谱密度。技术实现要素：本发明的目的是针对现有技术中存在的技术缺陷，而提供一种基于相关图的循环平稳相关熵谱密度计算方法，通过综合利用相关图(correlogram)算法和toeplitz矩阵与hankel矩阵高效计算能力，提高了计算效率。为实现本发明的目的所采用的技术方案是：一种基于相关图的循环平稳相关熵谱密度计算方法，包括以下步骤：步骤1，给定信号x、核长σ、窗函数的采样点数nw、分段序列重叠的采样点数noverlap和计算fft时的采样点数nfft；步骤2，将信号x进行分段，得到各分段信号xxi，xwi的采样点数为nw；步骤3，计算各分段信号xxi的中心化相关熵步骤4，对加窗处理，得到步骤5，计算的toeplitz矩阵及其hankel矩阵步骤6，计算矩阵与矩阵的点积，得到各分段信号时变相关熵矩阵步骤7，计算所有分段信号的均值，得到时变相关熵矩阵vm；步骤8，对矩阵vm各行计算nfft点fft，得到循环相关熵函数矩阵ccef；步骤9，对矩阵ccef各列计算nfft点fft，得到循环平稳相关熵谱密度矩阵s。在上述技术方案中，所述步骤2中，按welch方法将信号x进行分段，得到各分段信号xxi。在上述技术方案中，所述步骤2中，各分段信号xxi的采样点数为nw，xxi计算公式为：r为连续两个分段平移的采样点数，r＝nw-noverlap；k为信号x分段的总数量，表示不大于(lx-r)/r的最大整数；lx为信号x的总采样点数，即信号x的长度。在上述技术方案中，所述步骤3中，其中κ[·]为高斯核函数，高斯核函数κ[xxi(j)-xxi(j-m)]的计算公式为：在上述技术方案中，所述步骤4中，其中w(·)表示长度为nw对称窗函数，例如hamming窗函数、hanning窗函数等。在上述技术方案中，所述步骤5中，矩阵和矩阵的计算公式为：其中：toeplitz(·)为toeplitz矩阵算子，表示创建以为矩阵元素且行数和列数都等于nw的方阵，矩阵的第一行和第一列皆为而且矩阵中任一条平行于主对角线的直线上的元素相同；hankel(·)为hankel矩阵算子，表示创建以为矩阵元素且行数和列数都等于nw的方阵，矩阵的第一行和第一列皆为而且矩阵中任一条平行于副对角线的直线上的元素相同。在上述技术方案中，所述步骤6中，的计算公式为：其中.*为点积算子。在上述技术方案中，所述步骤7中，vm的计算公式为：在上述技术方案中，所述步骤8中，ccef的计算公式为：ccef＝fft(vm,nfft,2)，其中：fft(·)为快速傅里叶变换算子。在上述技术方案中，所述步骤9中，s的计算公式为：s＝fft(ccef,nfft,1)。与现有技术相比，本发明的有益效果是：1.本发明利用toeplitz矩阵和hankel矩阵快速计算时变相关熵矩阵，避免了传统方法计算相关熵矩阵采用多重循环语句、计算速度慢的弊端。2.本发明提供的基于相关图的循环平稳相关熵谱密度计算方法，计算速度快且占用计算机内存少。附图说明图1为本发明所述方法的流程图。图2为实施例2仿真调幅信号x(t)的时域波形。图3为实施例2仿真调幅信号x(t)的傅里叶变换(fft)。图4为实施例2.1仿真调幅信号x(t)，当σ＝0.3时的循环平稳相关熵谱密度的平面轮廓图。图5为实施例2.1仿真调幅信号x(t)，当σ＝0.3时的循环平稳相关熵谱密度的三维立体图。图6为实施例2.2仿真调幅信号x(t)，当σ＝1时的循环平稳相关熵谱密度的平面轮廓图。图7为实施例2.2仿真调幅信号x(t)，当σ＝1时的循环平稳相关熵谱密度的三维立体图。图8为实施例2.3仿真调幅信号x(t)，当σ＝5时的循环平稳相关熵谱密度的平面轮廓图。图9为实施例2.3仿真调幅信号x(t)，当σ＝5时的循环平稳相关熵谱密度的三维立体图。图10为对比例1仿真调幅信号x(t)，根据cpd算法计算的循环平稳相关熵谱密度平面轮廓图(σ＝0.3)。图11为对比例1仿真调幅信号x(t)，根据cpd算法计算的循环平稳相关熵谱密度三维立体图(σ＝0.3)；具体实施方式以下结合具体实施例对本发明作进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。实施例1步骤1，给定信号x、核长σ、窗函数的采样点数nw、分段序列重叠的采样点数noverlap和计算fft时的采样点数nfft。步骤2，将信号x进行分段，得到各分段信号xxi，xwi的采样点数为nw；按welch方法将信号x进行分段，得到各分段信号xxi，各分段信号xxi的采样点数为nw，xxi计算公式为：r为连续两个分段平移的采样点数，r＝nw-noverlap；k为信号x分段的总数量，表示不大于(lx-r)/r的最大整数；lx为信号x的总采样点数，即信号x的长度。步骤3，计算各分段信号xxi的中心化相关熵其中κ[·]为高斯核函数，高斯核函数κ[xxi(j)-xxi(j-m)]的计算公式为：σ为高斯核函数的核长，即为步骤1中给定的核长σ。步骤4，对加窗处理，得到其中w(·)表示长度为nw对称窗函数，例如hamming窗函数、hanning窗函数等。步骤5，计算的toeplitz矩阵及其hankel矩阵矩阵和矩阵的计算公式为：其中：toeplitz(·)为toeplitz矩阵算子，表示创建以为矩阵元素且行数和列数都等于nw的方阵，矩阵的第一行和第一列皆为而且矩阵中任一条平行于主对角线的直线上的元素相同；hankel(·)为hankel矩阵算子，表示创建以为矩阵元素且行数和列数都等于nw的方阵，矩阵的第一行和第一列皆为而且矩阵中任一条平行于副对角线的直线上的元素相同。步骤6，计算矩阵与矩阵的点积，得到各分段信号时变相关熵矩阵的计算公式为：其中.*为点积算子。步骤7，计算所有分段信号的均值，得到时变相关熵矩阵vm；vm的计算公式为：步骤8，对矩阵vm各行计算nfft点fft，得到循环相关熵函数矩阵ccef；ccef的计算公式为：ccef＝fft(vm,nfft,2)，其中：fft(·)为快速傅里叶变换算子。步骤9，对矩阵ccef各列计算nfft点fft，得到循环平稳相关熵谱密度矩阵s；s的计算公式为：s＝fft(ccef,nfft,1)。实施例2本实施例中，仿真调幅信号为x(t)＝a(1+cos(2πf0t))cos(2πfct)，采样时间t＝0.5s，采样频率为fs＝6000hz，信号x(t)长度n＝3000，调幅频率f0＝120hz，载波频率fc＝1000hz，幅值a＝1。窗函数的采样点数nw＝512、计算fft时的采样点数nfft＝1024和分段序列重叠的采样点数noverlap＝3nw/4。本实施例仿真调幅信号x(t)的时域波形，如图2所示；本实施例仿真调幅信号x(t)的频域波形，如图3所示。为验证核长σ对循环平稳相关熵谱密度的影响，本实施例给出了几种不同σ值时，循环平稳相关熵谱密度的平面轮廓图及其三维立体图。本实施例中，利用实施例1的方法对信号x(t)进行计算处理,不同实施例的σ值如表1所示：表1不同实施例的σ取值实施例σ2.10.32.212.35图4和图5为实施例2.1仿真调幅信号x(t)，当σ＝0.3时的循环平稳相关熵谱密度的平面轮廓图和三维立体图。图6和图7为实施例2.2仿真调幅信号x(t)，当σ＝1时的循环平稳相关熵谱密度的平面轮廓图和三维立体图。图8和图9为实施例2.3仿真调幅信号x(t)，当σ＝5时的循环平稳相关熵谱密度的平面轮廓图和三维立体图。对比图4、图6和图8可以看出：随着σ的增大，循环相关熵谱密度图中包含调幅频率f0的高次谐波成分在逐渐减少。当σ的较大(如σ≥5)时，循环相关熵谱密度图退化为传统谱相关密度图。对于实施例2仿真调幅信号x(t)和给定参数，采用计算机参数为：inteli7-2600，双核cpu3.4ghz，内存3gb，windowxp操作系统，使用的软件为matlabrelease12.1，程序运行时间为0.859s。对比例1为对比本发明提出的方法与fontes提出的循环平稳相关熵谱密度估计cpd方法(aluiso.i.r.fontes,etal,cyclostationarycorrentropy:definitionandapplication,expertsystemswithapplication.2017，69(1)：110-117)的计算效率，对于实施例2仿真调幅信号x(t)和给定参数，图10和图11分别给出了采用cpd方法计算的循环平稳相关熵谱密度的平面轮廓图及其三维立体图，程序运行时间为1707.672s。因此。采用本发明提出的方法，大大提高了循环平稳相关熵谱密度计算速度。对比图4、图5与图10、图11可以看出：本发明方法计算的循环平稳相关熵谱密度，在频谱图中是相互分离的谱峰，频谱分辨率高、无频谱泄漏现象；而cpd方法的计算结果谱峰之间是互连的，表明存在比较严重的频谱泄漏现象，频谱分辨率较低。因此。采用本发明提出的方法，不仅提高了频谱分辨率，而且还避免了频谱泄漏现象。以上所述仅是本发明的优选实施方式，应当指出的是，对于本
技术领域：
的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。当前第1页12