非平稳信号的自适应局部迭代滤波（MATLAB）

仍以滚动轴承故障诊断为例，在滚动轴承的运行过程中，其振动信号包含了大量的系统运行状态信息。利用振动信号进行滚动轴承的故障诊断，实际上就是分析振动信号、提取信息的过程。由于非线性力的作用，滚动轴承的故障信号往往是非平稳、非线性的多分量信号。传统的频域分析是以傅里叶变换为基础的全局变换，只适用于用于平稳性信号的分析，使得其应用受到限制。时频分析是以信号局部变换为基础，适用于非平稳信号的分析，因此，时频分析方法是研究的热点问题之一。随着信号处理技术的发展，基于信号分解的时频分析方法已经得到了广泛的应用。这种方法的基本思路是先将信号进行分解，得到瞬时频率具有物理意义的若干个分量信号和一个余量信号。对各个分量分别进行时频变换，将得到的时频分布进行组合，得到整个信号的时频分布。这其中，最重要的一步便是信号分解方法。信号分解方法主要有两种思路：基于迭代的分解方法和基于优化的分解方法。

基于迭代的信号分解方法是应用最为广泛的一类方法，其中最经典的是EMD 方法。但是，EMD 方法也存在很多缺陷，例如欠包络、过包络、端点效应、模态混叠等。基于这样的原因，迭代滤波方法被提出，迭代滤波沿用了EMD 方法的分解框架，采用信号的全局极值尺度来构建滤波器，利用滤波函数与信号卷积来定义信号的均值曲线，通过迭代筛分来分解信号。对迭代滤波的筛分过程收敛性，提出了严格的数学证明。但是，迭代滤波方法是基于全局极值尺度来构建滤波器的，不适合非平稳信号的分析。因此自适应局部迭代滤波算法被提出，利用福克普朗克方程来构建滤波器，通过信号的局部极值尺度来获取滤波器参数，将IF方法推广到非平稳信号的分析中。与EMD等方法相比，自适应局部迭代滤波方法具有更强的数学基础，同时在分解精度、端点效应等方面，相比其他迭代方法具有一定的优势。

工学博士，担任《Mechanical System and Signal Processing》《中国电机工程学报》《控制与决策》等期刊审稿专家，擅长领域：现代信号处理，机器学习，深度学习，数字孪生，时间序列分析，设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。


function [IMF,mask_lengths] = ALIF(f,options,mask_lengths_given)
% 
%% We deal with the input

if nargin == 0,  help ALIFv5_1; return; end
if nargin == 1, options = Settings_ALIF; end
if nargin == 2, mask_lengths_given=[]; end

extensionType = 'p'; % used in the calculations of mins and maxs

N = length(f);
if size(f,1)>size(f,2)
    f = f.';
end
if size(f,1)>1
    disp('Wrong dataset, the signal must be a single row vector')
end
IMF =[];

if options.plots>0
    if options.saveplots==1
        nameFile=input('Please enter the name of the file as a string using '' and '' <<  '); %'v081_Ex2';
    end
end

%% Main code

fprintf(['\n\n         ****************** WARNING ******************\n\n '...
    'We assume periodicity in the signal and\n\n  its instantaneous periods\n\n' ...
    '         *********************************************\n'])
pause(0.5)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Adaptive Local Iterative Filtering           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic

load('prefixed_double_filter','MM');
if length(f)>length(MM)
    fprintf(['\n\n      ********   Warning  *********\n\n'...
        ' The filter MM should contain more points\n'...
        ' to properly decompose the given signal\n\n'...
        ' We will use interpolation to generate a\n'...
        ' filter with the proper number of points\n\n'...
        '      *****************************\n\n'])
end

%%%%%%%%%%%%%%% Identify extreme points %%%%%%%%%%%%%%
maxmins_f=Maxmins_v2(f,extensionType);


while  length(maxmins_f) > (options.ALIF.ExtPoints) && size(IMF,1) < (options.ALIF.NIMFs) && toc < (options.maxTime)
    %% Outer loop
    
    if options.verbose>0
        fprintf('\n   ======================== IMF # %1.0d ========================\n',size(IMF,1)+1)
    end
    
    h = f;
    
    %%%%%%%%%%  Computing the mask length  %%%%%%%%%%%%%%%%
    
    if not(isempty(mask_lengths_given)) && size(mask_lengths_given,1)>size(IMF,1)
        if options.verbose>0
            fprintf('\n    ---------- Mask length provided by the user ----------\n')
        end
        
        iT_f=options.ALIF.xi*mask_lengths_given(size(IMF,1)+1,:);
    else
        if options.verbose>0
            fprintf('\n    --------------- Mask length computation ---------------\n')
        end
        
        T_f=[diff(maxmins_f) (maxmins_f(1)+N-maxmins_f(end))];
        temp_T_f=[T_f T_f T_f T_f T_f T_f T_f T_f T_f T_f T_f];
        temp_maxmins_f=[maxmins_f maxmins_f+N maxmins_f+2*N maxmins_f+3*N ...
            maxmins_f+4*N maxmins_f+5*N maxmins_f+6*N ...
            maxmins_f+7*N maxmins_f+8*N maxmins_f+9*N maxmins_f+10*N];
        temp_iT_f= interp1(temp_maxmins_f,temp_T_f,1:11*N,'cubic');
        iT_f = temp_iT_f(5*N+1:6*N);
        
        
        nTry=1;
        
        iT_f0=iT_f;
        
        if options.verbose>0
            fprintf('\n       ~~~~~~~~~~~~~~ Mask length using IF ~~~~~~~~~~~~~~\n')
        end
        
        OK=0;
        while OK==0
            opts=Settings_IF('IF.ExtPoints',3,'IF.NIMFs',nTry,'verbose',options.verbose,'IF.alpha',1,'IF.extensionType','p');
            IMF_iT_f = IF_v2(iT_f0,opts);  % We use IF algo for periodic signals to compute the mask length
            if 0>=min(IMF_iT_f(end,:)) && (size(IMF_iT_f,1)-1)==nTry
                nTry=nTry+1;
            elseif 0>=min(IMF_iT_f(end,:)) && not((size(IMF_iT_f,1)-1)==nTry)
                disp('Negative mask length')
                return
            else
                OK=1;
            end
        end
        
        if options.verbose>0
            fprintf('\n       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n')
        end
        
        iT_f = IMF_iT_f(end,:);
        nn = 2*options.ALIF.xi;
        iT_f = nn*iT_f;
        
        if options.plots>0
            figMask=figure;
            title(['Mask length IMF_{' num2str(size(IMF,1)+1) '}'])
            plot(iT_f,'b')
            hold on
            plot(maxmins_f,nn*T_f,'kx')
            plot(nn*sum(IMF_iT_f,1),'r')
            legend('Mask length','periods of the signal','Instantaneous periods')
            set(figMask,'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
            if options.saveplots==1
                saveas(figMask,[nameFile '_MaskLength_IMF_' num2str(size(IMF,1)+1)], 'fig')
                saveas(figMask,[nameFile '_MaskLength_IMF_' num2str(size(IMF,1)+1)], 'epsc')
                saveas(figMask,[nameFile '_MaskLength_IMF_' num2str(size(IMF,1)+1)], 'png')
            end
            pause(0.01)
        end
    end
    
    mask_lengths(size(IMF,1)+1,:)=iT_f;
    
    if options.verbose>0
        fprintf('\n    ------------------------------------------------------\n')
    end
    
    inStepN=0;
    
    W=zeros(N,N);
    for i=1:N
        wn = get_mask_v1(MM, iT_f(i));
        wn=wn/norm(wn,1);
        len = (length(wn)-1)/2;
        wn = [reshape(wn,1,2*len+1) zeros(1,N-2*len-1)];
        W(i,:)=circshift(wn,[0 (i-len-1)]);
    end
    
    %%
    
    if options.verbose>0
        fprintf('\n IMF # %1.0d\n',size(IMF,1)+1)
        fprintf('\n  step #            SD\n\n')
    end
    
    SD=Inf;
    
    if options.plots>0
        figM=figure;
    end
    while SD > options.ALIF.delta && inStepN<1000
        %% Inner loop
        
        inStepN=inStepN+1;
        
        if options.verbose>0
            
        end
        
        %%%%%%%% Compute the average %%%%%%%%%
        
        ave = W*h';
        
        %%%%%%%%%%%%%%%%%% generating h_n %%%%%%%%%%%%%%%%%%
        
        SD = norm(ave)^2/norm(h)^2;
        h = h - ave';
        
        if options.verbose>0
            fprintf('    %2.0d      %1.14f\n',inStepN,SD)
        end
        
        if options.plots>0 && rem(inStepN,1)==0
            
            textLeg{1}='f-h';
            textLeg{2}='h';
            titLeg=sprintf('IMF # %2.0d  Inner step =   %2.0d',size(IMF,1)+1,inStepN);
            
            title(sprintf('IMF # %2.0d  Inner step =   %2.0d',size(IMF,1)+1,inStepN))
            plot_imf_v8([h;f-h],1:length(h),titLeg,textLeg,figM);
            set(figM,'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
            
            if options.saveplots==1
                saveas(figM,[nameFile '_IMF' num2str(size(IMF,1)+1) '_step' num2str(inStepN)], 'fig')
                saveas(figM,[nameFile '_IMF' num2str(size(IMF,1)+1) '_step' num2str(inStepN)], 'epsc')
                saveas(figM,[nameFile '_IMF' num2str(size(IMF,1)+1) '_step' num2str(inStepN)], 'png')
            end
            pause(0.01)
        end
    end
    
    %%%%%%%%%%%%%%%%%%%%%% Adding the new IMF %%%%%%%%%%%%%%%%%%%%%
    IMF=[IMF; h];
    
    %%%%%%%%%%%%%%%%%%%%% Updating the signal %%%%%%%%%%%%%%%%%%%%%%
    f = f-h;
    
    %%%%%%%%%%%%%%% Identify extreme points %%%%%%%%%%%%%%
    maxmins_f=Maxmins_v2(f,extensionType);
    
    if options.saveEnd == 1
        save('Decomp_ALIF_v5.mat')
    end
    
end

%%%%%%%%%%%%%%%%%%%%%%%%% Adding the trend %%%%%%%%%%%%%%%%%%%%%%%
IMF =[IMF; f];

if options.saveEnd == 1
    save('Decomp_ALIF_v5.mat')
end

end


%% Auxiliar functions

function a=get_mask_v1(y,k)
% get the mask with length 2*k+1
% k could be an integer or not an integer
% y is the area under the curve for each bar

n=length(y);
m=(n-1)/2;

if k<=m % The prefixed filter contains enough points
    
    if mod(k,1)==0     % if the mask_length is an integer
        
        a=zeros(1,2*k+1);
        
        for i=1:2*k+1
            s=(i-1)*(2*m+1)/(2*k+1)+1;
            t=i*(2*m+1)/(2*k+1);
            
            %s1=s-floor(s);
            s2=ceil(s)-s;
            
            t1=t-floor(t);
            %t2=ceil(t)-t;
            
            if floor(t)<1
                disp('Ops')
            end
            a(i)=sum(y(ceil(s):floor(t)))+s2*y(ceil(s))+t1*y(floor(t));
        end
        
    else   % if the mask length is not an integer
        new_k=floor(k);
        extra = k-new_k;
        c=(2*m+1)/(2*new_k+1+2*extra);
        
        a=zeros(1,2*new_k+3);
        
        t=extra*c+1;
        t1=t-floor(t);
        %t2=ceil(t)-t;
        if k<0
            disp('Ops')
        end
        a(1)=sum(y(1:floor(t)))+t1*y(floor(t));
        
        for i=2:2*new_k+2
            s=extra*c+(i-2)*c+1;
            t=extra*c+(i-1)*c;
            %s1=s-floor(s);
            s2=ceil(s)-s;
            
            t1=t-floor(t);
            
            a(i)=sum(y(ceil(s):floor(t)))+s2*y(ceil(s))+t1*y(floor(t));
        end
        t2=ceil(t)-t;
        
        a(2*new_k+3)=sum(y(ceil(t):n))+t2*y(ceil(t));
    end
else % We need a filter with more points than MM, we use interpolation
    dx=0.01;
    % we assume that MM has a dx = 0.01, if m = 6200 it correspond to a
    % filter of length 62*2 in the physical space
    f=y/dx; % function we need to interpolate
    dy=m*dx/k;
    b=interp1(0:m,f(m+1:2*m+1),0:m/k:m);
    if size(b,1)>size(b,2)
        b=b.';
    end
    if size(b,1)>1
        fprintf('\n\nError!')
        disp('The provided mask is not a vector!!')
        a=[];
        return
    end
    a=[fliplr(b(2:end)) b]*dy;
    if abs(norm(a,1)-1)>10^-14
        %         fprintf('\n\nError\n\n')
        %         fprintf('Area under the mask = %2.20f\n',norm(a,1))
        %         fprintf('it should be equal to 1\nWe rescale it using its norm 1\n\n')
        a=a/norm(a,1);
    end
end

end


function varargout = Maxmins(f,extensionType)


% Identify the maxima and minima of a signal f

if nargin == 1, extensionType = 'c'; end
N = length(f);
maxmins=zeros(1,N);
df = diff(f);


h = 1;
cIn=0;
if strcmp(extensionType,'p') && df(1) == 0 && df(end) == 0
    while df(h)==0
        cIn=cIn+1;
        h=h+1;
    end
end

c = 0;
cmaxmins=0;
for i=h:N-2
    if   df(i)*df(i+1) <= 0
        if df(i+1) == 0
            if c == 0
                posc = i;
            end
            c = c + 1;
        else
            if c > 0
                cmaxmins=cmaxmins+1;
                maxmins(cmaxmins)=posc+floor((c-1)/2)+1;
                c = 0;
            else
                cmaxmins=cmaxmins+1;
                maxmins(cmaxmins)=i+1;
            end
        end
    end
end
if c > 0
    cmaxmins=cmaxmins+1;
    maxmins(cmaxmins)=mod(posc+floor((c+cIn-1)/2)+1,N);
    if maxmins(cmaxmins)==0
        maxmins(cmaxmins)=N;
    end
end

maxmins=maxmins(1:cmaxmins);

if strcmp(extensionType,'p') % we deal with a periodic signal
    if isempty(maxmins)
        maxmins = 1;
    else
        if maxmins(1)~=1 && maxmins(end)~=N
            if (f(maxmins(end)) > f(maxmins(end)+1) && f(maxmins(1)) > f(maxmins(1)-1)) || (f(maxmins(end)) < f(maxmins(end)+1) && f(maxmins(1)) < f(maxmins(1)-1))
                maxmins=[1 maxmins];
            end
        end
    end
elseif strcmp(extensionType,'c')
    if isempty(maxmins)
        maxmins = [1, N];
    else
        if maxmins(1) ~= f(1) && maxmins(end) ~= f(end)
            maxmins = [1, maxmins, N];
        elseif f(maxmins(1)) ~= f(1)
            maxmins = [1, maxmins];
        elseif  f(maxmins(end)) ~= f(end)
            maxmins = [maxmins, N];
        end
    end
elseif strcmp(extensionType,'r')
    if isempty(maxmins)
        maxmins = [1, N];
    else
        if maxmins(1) ~= f(1) && maxmins(end) ~= f(end)
            maxmins = [1, maxmins, N];
        elseif f(maxmins(1)) ~= f(1)
            maxmins = [1, maxmins];
        elseif  f(maxmins(end)) ~= f(end)
            maxmins = [maxmins, N];
        end
    end
end

if nargout<=1
    varargout{1}=maxmins;
elseif nargout==2
    if f(maxmins(1))>f(maxmins(2))
        % Maxes
        varargout{1}=maxmins(1:2:end);
        % Mins
        varargout{2}=maxmins(2:2:end);
    elseif f(maxmins(2))>f(maxmins(1))
        % Maxes
        varargout{1}=maxmins(2:2:end);
        % Mins
        varargout{2}=maxmins(1:2:end);
    else
        disp('Problem in identifying Extrema')
        varargout{1}=[];
    end
end

end

知乎学术咨询：
https://www.zhihu.com/consult/people/792359672131756032?isMe=1