非平稳信号的自适应局部迭代滤波(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 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

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.xdnf.cn/news/1473559.html

如若内容造成侵权/违法违规/事实不符,请联系一条长河网进行投诉反馈,一经查实,立即删除!

相关文章

记录一次麒麟V10 安装sysbench各种报错(关于MySQL)处理过程

sysbench手工下载&#xff1a; https://github.com/akopytov/sysbench 下载.zip文件&#xff0c;上传到服务器上 解压、安装&#xff1a; unzip sysbench-master.zipcd sysbench-master/sh autogen.sh./configure 报错&#xff1a;没有mysql驱动 configure: error: mysql_c…

UEC++ 虚幻5第三人称射击游戏(二)

UEC++ 虚幻5第三人称射击游戏(二) 派生榴弹类武器 新建一个继承自Weapon的子类作为派生榴弹类武器 将Weapon类中的Fire函数添加virtual关键字变为虚函数让榴弹类继承重写 在ProjectileWeapon中重写Fire函数,新建生成投射物的模版变量 Fire函数重写逻辑 代码//生成的投射物U…

关注推送---Feed流,推模式实现的个人分析及其思考。

本篇文章记录我们实际开发过程中&#xff0c;关注推送场景的个人思考&#xff0c;以及解析。 文章目录 前言一、关注推送是什么&#xff1f;是什么是Feed流&#xff1f;二、解决关注推送问题的技术方案1.理论模型的选取2.数据类型的选取 三、理论模型的选取三、数据类型的选取总…

前端八股文 对$nextTick的理解

$nexttick是什么? 获取更新后的dom内容 为什么会有$nexttick ? vue的异步更新策略 (这也是vue的优化之一 要不然一修改数据就更新dom 会造成大量的dom更新 浪费性能) 这是因为 message &#xff08;data&#xff09;数据在发现变化的时候&#xff0c;vue 并不会立刻去更…

echarts阶段仪表图

echarts阶段仪表图 – 效率图 1、先上效果展示 2、完整源码奉上 Vue2 echarts 5 <template><div ref"gaugeChart" style"width: 100%; height: 100%"></div> </template><script> import * as echarts from "echar…

智能化客户服务:提升效率与体验的新模式

在数字化浪潮的推动下&#xff0c;客户服务领域正经历着一场深刻的变革。智能化客户服务的兴起&#xff0c;不仅重塑了企业与客户之间的互动方式&#xff0c;更在提升服务效率与增强客户体验方面展现出了巨大潜力。本文将深入探讨智能化客户服务的新模式&#xff0c;分析其如何…

Arc for Windows 无法使用?一篇文章教会你!

&#x1f44b; 大家好&#xff0c;我是 Beast Cheng &#x1f4eb; 联系我&#xff1a;458290771qq.com &#x1f331; 接合作、推广…… 什么是Arc浏览器&#xff1f; Arc浏览器是The Browser Conpany使用Swift语言开发的一款浏览器&#xff0c;Arc浏览器由其漂亮的侧边栏闻名…

上万组风电,光伏,用户负荷数据分享

上万组风电&#xff0c;光伏&#xff0c;用户负荷数据分享 可用于风光负荷预测等研究 获取链接&#x1f517; https://pan.baidu.com/s/1izpymx6R3Y8JsFdx42rL0A 提取码&#xff1a;381i 获取链接&#x1f517; https://pan.baidu.com/s/1izpymx6R3Y8JsFdx42rL0A 提取…

vue2响应式原理+模拟实现v-model

效果 简述原理 配置对象传入vue实例 模板解析&#xff0c;遍历出所有文本节点&#xff0c;利用正则替换插值表达式为真实数据 data数据代理给vue实例&#xff0c;以后通过this.xxx访问 给每个dom节点增加观察者实例&#xff0c;由观察者群组管理&#xff0c;内部每一个键值…

django高校教务系统-计算机毕业设计源码81661

目 录 摘要 1 绪论 1.1 研究背景 1.2目的及意义 1.3论文结构与章节安排 2 高校教务系统设计分析 2.1 可行性分析 2.1.1 技术可行性分析 2.1.2 经济可行性分析 2.1.3 法律可行性分析 2.2 系统功能分析 2.2.1 功能性分析 2.2.2 非功能性分析 2.3 系统用例分析 2.4…

【云WAF为您的Web防御保驾护航】

在这个数字时代&#xff0c;网络就像是一张没有尽头的大网&#xff0c;将整个世界都联系在了一起。但是&#xff0c;在这个网络的背后&#xff0c;却潜藏着数不清的安全隐患。恶意攻击、数据泄漏、网站瘫痪……各种隐患就像是隐藏在暗处的毒蛇&#xff0c;时刻都会对没有任何防…

金蝶云苍穹-插件开发(一)加载数据

前言 此系列博客是进行金蝶云苍穹开发时的插件开发的教程&#xff0c;一是在明年要是还要参加软件杯金蝶A6赛题的话&#xff0c;可以看此系列教程的博客来进行复习&#xff0c;同时如果要是我实验室的学弟学妹要参加的话&#xff0c;我这个系列的博客可以给他们提供学习参考&a…

破解宇宙终极奥秘,战胜昊天无上束缚

在幽邃的暗夜下&#xff0c;细品着夫子与昊天跨越千年的智勇交锋&#xff0c;我的思绪不禁飘向了更加深远的宇宙边际&#xff0c;回响起那些关于人类如何挑战天命、战胜上天的过往。 宇宙奥秘 在浩瀚无垠的宇宙深渊中&#xff0c;隐藏着一段超越凡尘的规则。昊天&#xff0c;…

2024/7/6 英语每日一段

More than half of late-teens are specifically calling for more youth work that offers “fun”, with older teenagers particularly hankering for more jollity, according to a study carried out by the National Youth Agency. One in 10 said they have zero option…

【Go】excelize库实现excel导入导出封装(四),导出时自定义某一列或多列的单元格样式

大家好&#xff0c;这里是符华~ 查看前三篇&#xff1a; 【Go】excelize库实现excel导入导出封装&#xff08;一&#xff09;&#xff0c;自定义导出样式、隔行背景色、自适应行高、动态导出指定列、动态更改表头 【Go】excelize库实现excel导入导出封装&#xff08;二&…

AI Earth应用—— 在线使用sentinel数据VV和VH波段进行水体提取分析(昆明抚仙湖、滇池为例)

AI Earth 本文的主要目的就是对水体进行提取,这里,具体的操作步骤很简单基本上是通过,首页的数据检索,选择需要研究的区域,然后选择工具箱种的水体提取分析即可,剩下的就交给阿里云去处理,结果如下: 这是我所选取的一景影像: 详情 卫星: Sentinel-1 级别: 1 …

新产品或敏捷项目过程 SOP,附带流程图及流程规范

一、项目启动 项目背景和目标明确 市场调研结果分析&#xff0c;确定新产品的需求和市场机会。制定明确的项目目标&#xff0c;包括产品特性、上市时间、预期收益等。 组建项目团队 确定项目经理、产品经理、开发人员、测试人员、市场人员等角色。明确各成员的职责和权限。 项目…

Java语言程序设计基础篇(第10版)编程练习题13.18(使用 Rational 类)

第十三章第十八题(使用 Rational 类) 题目要求&#xff1a; 编写程序&#xff0c;使用 Rational 类计算下面的求和数列: 你将会发现输出是不正确的 &#xff0c;因为整数溢出(太大了)。为了解决这个问题 &#xff0c;参见编程练习題13.15。代码参考&#xff1a; package cha…

蜂窝物联农业气象站,守护丰收每一步

现代农业的革新者——农业自动气象站&#xff0c;正以其多功能的传感器、高效的数据采集传输系统、智能的数据云平台以及可靠的供电供网系统&#xff0c;成为农业生产中的得力助手。这些传感器能够实时监测温度、湿度、风速、风向、气压、土壤温度、土壤湿度、土壤PH值、土壤盐…

CorelDRAW2024新版本来咯!你的设计神助手

&#x1f389; 设计界的朋友们&#xff0c;注意啦&#xff01;你们的新宠——CorelDRAW 2024 来咯&#xff01; &#x1f31f; 一、设计神器再进化 亲爱的设计小伙伴们&#xff0c;有没有感觉每天与那些不配合的软件战斗&#xff0c;像是在打怪升级&#xff1f;&#x1f409; …