资讯详情

多尺度排列熵

文章目录

  • 前言
  • 一、多尺度排列熵是什么?
  • 二、实验平台照片
  • 三、MATLAB代码
    • 3.1 熵的多尺度排列
    • 3.2 排列熵
  • 参考文献


前言

齿轮和齿轮箱作为机械设备常用的旋转机械设备,不仅可以传递更大的功率和载荷,而且具有更好的可靠性。但在高精度切割加工中,当齿轮在变速、变载等复杂条件下工作时,容易损坏、磨损、断齿,大大降低了加工精度。因此,齿轮的状态监测和故障诊断尤为重要。


一、多尺度排列熵是什么?

排列熵是由 BANDT C 一维时间序列复杂度的平均熵参数通常用于提取机械故障的特征。然而,齿轮和齿轮箱运行过程中的故障特征信息分布在多尺度中,只有单一尺度的熵分析才会遗漏其他尺度的故障特征信息。排列熵( Permutation Entropy,PE) 测量时间序列复杂度和检测动力学突变的平均熵函数,算法简单,计算速度快。故障引起的高频冲击成分越多,熵的排列就越大。因此,在许多论文中引入熵排列指标作为评价不同滤波器长度滤波效果的标准。 针对排列熵的不足,COSTA M 提出了多尺度排列熵的概念。多尺度排列熵是在排列熵的基础上对时间序列进行多尺度粗粒化,然后计算不同尺度下粗粒化序列的排列熵。具体计算过程如下: 在这里插入图片描述

二、实验平台照片

在实验平台下验证本文提到的方法的有效性( 图2) ,采集试验台减速箱齿轮不同状态下的实时振动信号。齿轮减速箱型号为 JZQ采样频率为200 8 200 Hz。选择齿轮正常、不平衡、松动、脱位、断齿、裂纹 6 实验验证不同状态的数据。

三、MATLAB代码

3.1 熵的多尺度排列

function [pe hist] = pec(y,m,t)  %  Calculate the permutation entropy  %  Input:   y: time series; %           m: order of permuation entropy %           t: delay time of permuation entropy,   % Output:  %           pe:    permuation entropy %           hist:  the histogram for the order distribution    ly = length(y); permlist = perms(1:m); c(1:length(permlist))=0;       for j=1:ly-t*(m-1)      [a,iv]=sort(y(j:t:j t*(m-1)));      for jj=1:length(permlist)
         if (abs(permlist(jj,:)-iv))==0
             c(jj) = c(jj) + 1 ;
         end
     end
 end

hist = c;
 
c=c(find(c~=0));
p = c/sum(c);
pe = -sum(p .* log(p));
function outdata = PE( indata, delay, order, windowSize )
% @brief PE efficiently [1] computes values of permutation entropy [2] in 
% maximally overlapping sliding windows
%
% INPUT 
%   - indata - considered time series 
%   - delay - delay between points in ordinal patterns (delay = 1 means successive points)
%   - order - order of the ordinal patterns (order+1 - number of points in ordinal patterns)
%   - windowSize - size of sliding window
% OUTPUT 
%   - outdata - values of normalised permutation entropy as defined in [2]
%
% REFERENCES
% [1] Unakafova, V.A., Keller, K., 2013. Efficiently measuring complexity 
% on the basis of real-world data. Entropy, 15(10), 4392-4415.
% [2] Bandt C., Pompe B., 2002. Permutation entropy: a natural complexity 
% measure for time series. Physical review letters, APS
%
% @author Valentina Unakafova
% @date 10.11.2017
% @email UnakafovaValentina(at)gmail.com

 load( ['table' num2str( order ) '.mat'] ); % the precomputed table
 patternsTable = eval( ['table' num2str( order )] );    
 nPoints    = numel( indata );              % length of the time series
 opOrder1   = order + 1; 
 orderDelay = order*delay; 
 nPatterns  = factorial( opOrder1 );   % amount of ordinal patterns              
 patternsDistribution = zeros( 1, nPatterns ); % distribution of ordinal patterns
 outdata = zeros( 1, nPoints );        % values of permutation entropy    
 inversionNumbers = zeros( 1, order ); % inversion numbers of ordinal pattern (i1,i2,...,id)
 prevOP = zeros( 1, delay );           % previous ordinal patterns for 1:opDelay
 ordinalPatternsInWindow = zeros( 1, windowSize ); % ordinal patterns in the window
 ancNum = nPatterns./factorial( 2:opOrder1 );  % ancillary numbers       
 peTable( 1:windowSize ) = -( 1:windowSize ).*log( 1:windowSize  ); % table of precomputed values  
 peTable( 2:windowSize ) = ( peTable( 2:windowSize ) - peTable( 1:windowSize - 1 ) )./windowSize;       
 for iTau = 1:delay 
     cnt = iTau; 
     inversionNumbers( 1 ) = ( indata( orderDelay + iTau - delay ) >= indata( orderDelay + iTau ) );
     for j = 2:order
         inversionNumbers( j ) = sum( indata( ( order - j )*delay + iTau ) >= ...
           indata( ( opOrder1 - j )*delay + iTau:delay:orderDelay + iTau ) );
     end        
     ordinalPatternsInWindow( cnt ) = sum( inversionNumbers.*ancNum ); % first ordinal patterns
     patternsDistribution( ordinalPatternsInWindow( cnt )+ 1 ) = ...
       patternsDistribution( ordinalPatternsInWindow( cnt ) + 1 ) + 1;  
     for j = orderDelay + delay + iTau:delay:windowSize + orderDelay   % loop for the first window
         cnt = cnt + delay;                               
         posL = 1; % the position of the next point
         for i = j - orderDelay:delay:j - delay
             if( indata( i ) >= indata( j ) ) 
                 posL = posL + 1; 
             end
         end  
         ordinalPatternsInWindow( cnt ) = ...
           patternsTable( ordinalPatternsInWindow( cnt - delay )*opOrder1 + posL );
         patternsDistribution( ordinalPatternsInWindow( cnt ) + 1 ) = ...
           patternsDistribution( ordinalPatternsInWindow( cnt ) + 1 ) + 1;            
     end  
     prevOP( iTau ) = ordinalPatternsInWindow( cnt );
 end    
 ordDistNorm = patternsDistribution/windowSize;
 tempSum = 0;
 for iPattern = 1:nPatterns
   if ( ordDistNorm( iPattern ) ~= 0 )
     tempSum = tempSum - ordDistNorm( iPattern )*log( ordDistNorm( iPattern ) );
   end
 end
 outdata( windowSize + delay*order ) = tempSum;

 iTau = mod( windowSize, delay ) + 1;   % current shift 1:delay
 patternPosition = 1;                   % position of the current pattern in the window
 for t = windowSize + delay*order + 1:nPoints % loop over all points
     posL = 1;                          % the position of the next point
     for j = t-orderDelay:delay:t-delay
         if( indata( j ) >= indata( t ) ) 
             posL = posL + 1; 
         end
     end                          
     nNew = patternsTable( prevOP( iTau )*opOrder1 + posL ); % incoming ordinal pattern         
     nOut = ordinalPatternsInWindow( patternPosition );      % outcoming ordinal pattern 
     prevOP( iTau ) = nNew;
     ordinalPatternsInWindow( patternPosition ) = nNew; 
     nNew = nNew + 1;
     nOut = nOut + 1;       
     if ( nNew ~= nOut ) % update the distribution of ordinal patterns    
         patternsDistribution( nNew ) = patternsDistribution( nNew ) + 1; % incoming ordinal pattern
         patternsDistribution( nOut ) = patternsDistribution( nOut ) - 1; % outcoming ordinal pattern
         outdata( t ) = outdata( t - 1 ) + ( peTable( patternsDistribution( nNew ) ) - ...
                                        peTable( patternsDistribution( nOut ) + 1  ) );
     else
         outdata( t ) = outdata( t - 1 );
     end      
     iTau = iTau + 1; 
     patternPosition = patternPosition + 1;
     if ( iTau > delay ) 
       iTau = 1; 
     end
     if ( patternPosition > windowSize ) 
       patternPosition = 1; 
     end
 end 
 outdata = outdata( windowSize + delay*order:end )/log( factorial( order + 1 ) );

3.2 排列熵

%% 主函数调用模糊熵函数求时间序列的模糊熵值
clc;clear;
close all;
tic;
% 产生仿真信号
fs=100;   %数据采样率Hz
t=1:1/fs:4096*1/fs; %对数据进行采样
n = length(t);  %数据的采样数?
f1 =0.25; %信号的频率

f2=0.005;
x=2*sin(2*pi*f1*t+cos(2*pi*f2*t)); %产?原始信号
nt=0.2*randn(1,n); %?斯?噪声?成
y=x+nt; %含噪信号
% EEMD分解
Nstd=0.2;
NE=20;
X=eemd(y,Nstd,NE,100); % EEMD分解函数在本?的资源?可供下载
%%
% 相空间重构:eDim为嵌?维数
% 当X具有多列和多?时,每列将被视为独?的时间序列,该算法对X的每?列假设相同的嵌?维度和时间延迟,并以标量返回ESTDIM和ESTLAG。
[~,~,eDim] = phaseSpaceReconstruction(X');
% 求时间序列X的模糊熵,模糊熵的输?时间序列为?向量
X=X'; % 将信号y的各个分量转置
[m,n]=size(X);
r0=0.15; % r为相似容限度
Out_FuzEn=zeros(1,m);
for i=1:m
r=r0*std(X(i,:));
FuzEn(i) = FuzzyEntropy(X(i,:),eDim,r,2,1);
end
toc;
%% 模糊熵函数
function FuzEn = FuzzyEntropy(data,dim,r,n,tau)
%
% This function calculates fuzzy entropy (FuzEn) of a univariate signal data
%
% Inputs:
%
% data: univariate signal - a vector of size 1 x N (the number of sample points)
% dim: embedding dimension
% r: threshold (it is usually equal to 0.15 of the standard deviation of a signal - because we normalize signals to have a standard deviation of 1, here, r is usually equal to 0.15)
% n: fuzzy power (it is usually equal to 2)
% tau: time lag (it is usually equal to 1)
% 模糊熵算法的提出者:Chen Weiting,Wang Zhizhong,XieHongbo,et al. Characterization of surfaceEMG signal based on fuzzy entropy. IEEE Transactions on Neural Systems and Rehabilitation Engineering. 2007,15(2):266-272.
%
if nargin == 4, tau = 1; end
if nargin == 3, n = 2; tau=1; end
if tau > 1, data = downsample(data, tau); end
N = length(data);
result = zeros(1,2);
for m = dim:dim+1
count = zeros(N-m+1,1);
dataMat = zeros(N-m+1,m);
% 设置数据矩阵,构造成m维的?量
for i = 1:N-m+1
dataMat(i,:) = data(1,i:i+m-1);
end
% 利?距离计算相似模式数
for j = 1:N-m+1
% 计算切?雪夫距离,不包括?匹配情况
dataMat=dataMat-mean(dataMat,2);
tempmat=repmat(dataMat(j,:),N-m+1,1);
dist = max(abs(dataMat - tempmat),[],2);
D=exp(-(dist.^n)/r);
count(j) = (sum(D)-1)/(N-m);
end
result(m-dim+1) = sum(count)/(N-m+1);
end

% 计算得到的模糊熵值
FuzEn = log(result(1)/result(2));
end

%%
function [modos its]=eemd(x,Nstd,NR,MaxIter)
%--------------------------------------------------------------------------
%WARNING: this code needs to include in the same
%directoy the file emd.m developed by Rilling and Flandrin.
%This file is available at %http://perso.ens-lyon.fr/patrick.flandrin/emd.html
%We use the default stopping criterion.
%We use the last modification: 3.2007
% -------------------------------------------------------------------------
%  OUTPUT
%   modos: contain the obtained modes in a matrix with the rows being the modes
%   its: contain the iterations needed for each mode for each realization
%
%  INPUT
%  x: signal to decompose
%  Nstd: noise standard deviation
%  NR: number of realizations
%  MaxIter: maximum number of sifting iterations allowed.
% -------------------------------------------------------------------------
%   Syntax
%
%   modos=eemd(x,Nstd,NR,MaxIter)
%  [modos its]=eemd(x,Nstd,NR,MaxIter)
% ------------------------------------
        标签: keller压力变送器paa

锐单商城拥有海量元器件数据手册IC替代型号,打造 电子元器件IC百科大全!

锐单商城 - 一站式电子元器件采购平台