仍以滚动轴承故障诊断为例，在滚动轴承的运行过程中，其振动信号包含了大量的系统运行状态信息。利用振动信号进行滚动轴承的故障诊断，实际上就是分析振动信号、提取信息的过程。由于非线性力的作用，滚动轴承的故障信号往往是非平稳、非线性的多分量信号。传统的频域分析是以傅里叶变换为基础的全局变换，只适用于用于平稳性信号的分析，使得其应用受到限制。时频分析是以信号局部变换为基础，适用于非平稳信号的分析，因此，时频分析方法是研究的热点问题之一。随着信号处理技术的发展，基于信号分解的时频分析方法已经得到了广泛的应用。这种方法的基本思路是先将信号进行分解，得到瞬时频率具有物理意义的若干个分量信号和一个余量信号。对各个分量分别进行时频变换，将得到的时频分布进行组合，得到整个信号的时频分布。这其中，最重要的一步便是信号分解方法。信号分解方法主要有两种思路：基于迭代的分解方法和基于优化的分解方法。
基于迭代的信号分解方法是应用最为广泛的一类方法，其中最经典的是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 inputif nargin == 0, help ALIFv5_1; return; end
if nargin == 1, options = Settings_ALIF; end
if nargin == 2, mask_lengths_given=[]; endextensionType = 'p'; % used in the calculations of mins and maxsN = length(f);
if size(f,1)>size(f,2)f = f.';
end
if size(f,1)>1disp('Wrong dataset, the signal must be a single row vector')
end
IMF =[];if options.plots>0if options.saveplots==1nameFile=input('Please enter the name of the file as a string using '' and '' << '); %'v081_Ex2';end
end%% Main codefprintf(['\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 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ticload('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 loopif options.verbose>0fprintf('\n ======================== IMF # %1.0d ========================\n',size(IMF,1)+1)endh = f;%%%%%%%%%% Computing the mask length %%%%%%%%%%%%%%%%if not(isempty(mask_lengths_given)) && size(mask_lengths_given,1)>size(IMF,1)if options.verbose>0fprintf('\n ---------- Mask length provided by the user ----------\n')endiT_f=options.ALIF.xi*mask_lengths_given(size(IMF,1)+1,:);elseif options.verbose>0fprintf('\n --------------- Mask length computation ---------------\n')endT_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>0fprintf('\n ~~~~~~~~~~~~~~ Mask length using IF ~~~~~~~~~~~~~~\n')endOK=0;while OK==0opts=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 lengthif 0>=min(IMF_iT_f(end,:)) && (size(IMF_iT_f,1)-1)==nTrynTry=nTry+1;elseif 0>=min(IMF_iT_f(end,:)) && not((size(IMF_iT_f,1)-1)==nTry)disp('Negative mask length')returnelseOK=1;endendif options.verbose>0fprintf('\n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n')endiT_f = IMF_iT_f(end,:);nn = 2*options.ALIF.xi;iT_f = nn*iT_f;if options.plots>0figMask=figure;title(['Mask length IMF_{' num2str(size(IMF,1)+1) '}'])plot(iT_f,'b')hold onplot(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==1saveas(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')endpause(0.01)endendmask_lengths(size(IMF,1)+1,:)=iT_f;if options.verbose>0fprintf('\n ------------------------------------------------------\n')endinStepN=0;W=zeros(N,N);for i=1:Nwn = 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>0fprintf('\n IMF # %1.0d\n',size(IMF,1)+1)fprintf('\n step # SD\n\n')endSD=Inf;if options.plots>0figM=figure;endwhile SD > options.ALIF.delta && inStepN<1000%% Inner loopinStepN=inStepN+1;if options.verbose>0end%%%%%%%% Compute the average %%%%%%%%%ave = W*h';%%%%%%%%%%%%%%%%%% generating h_n %%%%%%%%%%%%%%%%%%SD = norm(ave)^2/norm(h)^2;h = h - ave';if options.verbose>0fprintf(' %2.0d %1.14f\n',inStepN,SD)endif options.plots>0 && rem(inStepN,1)==0textLeg{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==1saveas(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')endpause(0.01)endend%%%%%%%%%%%%%%%%%%%%%% 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 == 1save('Decomp_ALIF_v5.mat')endend%%%%%%%%%%%%%%%%%%%%%%%%% Adding the trend %%%%%%%%%%%%%%%%%%%%%%%
IMF =[IMF; f];if options.saveEnd == 1save('Decomp_ALIF_v5.mat')
endend%% Auxiliar functionsfunction 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 barn=length(y);
m=(n-1)/2;if k<=m % The prefixed filter contains enough pointsif mod(k,1)==0 % if the mask_length is an integera=zeros(1,2*k+1);for i=1:2*k+1s=(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)<1disp('Ops')enda(i)=sum(y(ceil(s):floor(t)))+s2*y(ceil(s))+t1*y(floor(t));endelse % if the mask length is not an integernew_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<0disp('Ops')enda(1)=sum(y(1:floor(t)))+t1*y(floor(t));for i=2:2*new_k+2s=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));endt2=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 interpolationdx=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 spacef=y/dx; % function we need to interpolatedy=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.';endif size(b,1)>1fprintf('\n\nError!')disp('The provided mask is not a vector!!')a=[];returnenda=[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
endendfunction varargout = Maxmins(f,extensionType)% Identify the maxima and minima of a signal fif 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) == 0while df(h)==0cIn=cIn+1;h=h+1;end
endc = 0;
cmaxmins=0;
for i=h:N-2if df(i)*df(i+1) <= 0if df(i+1) == 0if c == 0posc = i;endc = c + 1;elseif c > 0cmaxmins=cmaxmins+1;maxmins(cmaxmins)=posc+floor((c-1)/2)+1;c = 0;elsecmaxmins=cmaxmins+1;maxmins(cmaxmins)=i+1;endendend
end
if c > 0cmaxmins=cmaxmins+1;maxmins(cmaxmins)=mod(posc+floor((c+cIn-1)/2)+1,N);if maxmins(cmaxmins)==0maxmins(cmaxmins)=N;end
endmaxmins=maxmins(1:cmaxmins);if strcmp(extensionType,'p') % we deal with a periodic signalif isempty(maxmins)maxmins = 1;elseif maxmins(1)~=1 && maxmins(end)~=Nif (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];endendend
elseif strcmp(extensionType,'c')if isempty(maxmins)maxmins = [1, N];elseif 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];endend
elseif strcmp(extensionType,'r')if isempty(maxmins)maxmins = [1, N];elseif 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];endend
endif nargout<=1varargout{1}=maxmins;
elseif nargout==2if f(maxmins(1))>f(maxmins(2))% Maxesvarargout{1}=maxmins(1:2:end);% Minsvarargout{2}=maxmins(2:2:end);elseif f(maxmins(2))>f(maxmins(1))% Maxesvarargout{1}=maxmins(2:2:end);% Minsvarargout{2}=maxmins(1:2:end);elsedisp('Problem in identifying Extrema')varargout{1}=[];end
endend知乎学术咨询：
https://www.zhihu.com/consult/people/792359672131756032?isMe=1
