Fast_Kurtogram.m

function c = Fast_Kurtogram（xnlevelFsopt1opt2）
% Fast_Kurtogram（xnlevelFs）
% Computes the fast kurtogram of signal x up to level ‘nlevel‘
% Maximum number of decomposition levels is log2（length（x）） but it is 
% recommended to stay by a factor 1/8 below this.
% Fs = sampling frequency of signal x （default is Fs = 1）
% opt1 = 1: the kurtogram is computed via a fast decimated filterbank tree
% opt1 = 2: the kurtogram is computed via the short-time Fourier transform
% （if there is any difference in the kurtogram between the two measures this is
% due to the presence of impulsive additive noise）
% opt2 = 1: classical kurtosis based on 4th order statistics 
% opt2 = 2: robust kurtosis based on 2nd order statistics of the envelope
% （option 1 is faster and has more flexibility than option 2 in the design of the
% analysis filter: a short filter in option 1 gives virtually the same results as option 2）
%
% -------------------
% J. Antoni : 02/2005
% -------------------

N = length（x）;
N2 = log2（N） - 7;%最大允许分解层
if nlevel > N2%判断输入分解层是否大于最大允许分解层
   error（‘Please enter a smaller number of decomposition levels‘）;
end

if nargin < 5%判断输入变量
   opt2 = input（‘Choose the kurtosis measure （classic = 1 ; robust = 2）: ‘）;
   if nargin < 4
      opt1  = input（‘Choose the algorithm （filterbank = 1 ; stft-based = 2）: ‘）;
      if nargin < 3
         Fs = 1;
      end
   end
 end

% Fast computation of the kurtogram
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if opt1 == 1
   % 1） Filterbank-based kurtogram
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Analytic generating filters
   N = 16;			fc = .4;					% a short filter is just good enough!
   h = fir1（Nfc）.*exp（2i*pi*（0:N）*.125）;%二分段低通滤波器
   n = 2:N+1;
   g = h（1+mod（1-nN））.*（-1）.^（1-n）;%二分段高通滤波器
   % 
   N = fix（3/2*N）;
   h1 = fir1（N2/3*fc）.*exp（2i*pi*（0:N）*.25/3）;%三分段第一段滤波器
   h2 = h1.*exp（2i*pi*（0:N）/6）;%三分段第二段滤波器
   h3 = h1.*exp（2i*pi*（0:N）/3）;%三分段第三段滤波器  
   % 
   if opt2 == 1
      Kwav = K_wpQ（xhgh1h2h3nlevel‘kurt2‘）;				% kurtosis of the complex envelope
   else
      Kwav = K_wpQ（xhgh1h2h3nlevel‘kurt1‘）;				% variance of the envelope magnitude
   end
   Kwav = Kwav.*（Kwav>0）;												% keep positive values only!
   
else 
   % 2） STFT-based kurtogram
   %%%%%%%%%%%%%%%%%%%%%%%%%
   Nfft = 2.^[3:nlevel+2];				% level 1 of wav_kurt roughly corresponds to a 4-sample hanning window with stft_kurt
   %											  or a 8-sample flattop	第一层分8段	
   temp = [3*Nfft（1）/2 3*Nfft（1:end-2）;Nfft（2:end）];%第一行为三分层段数，第二行为二分层段数
   Nfft = [Nfft（1） temp（:）‘];%用矩阵计算而不是循环赋值的方法，所有分层段数
   if opt2 == 1
      Kstft = Kf_fft（xNfft1‘kurt2‘）;							% kurtosis of the complex envelope
      Kx = kurt（x‘kurt2‘）;%未分层前的信号峭度
   else
      Kstft = Kf_fft（xNfft1‘kurt1‘）;							% variance of the envelope magnitude
      Kx = kurt（x‘kurt1‘）;
   end
   Kstft = [Kx*ones（1