VMD分解,matlab代码,包络线,包络谱,中心频率,峭度值,能量熵,近似熵,包络熵,频谱图,希尔伯特变换,包含所有程序MATLAB代码,-西储大学数据集为例
目录
1.选取数据
2.VMD函数-matlab代码
3.采用matlab脚本导入数据并做VMD分解
4. VMD分解图
5.计算中心频率
6.画包络线
7. 画包络谱
8. 计算峭度值
9.计算能量和能量熵
10.近似熵代码
11.包络熵
12.频谱图
13.结果展示
13.智能算法优化VMD参数
1.选取数据
选取1797转速下的内圈故障数据,也就是105.mat,数据集可以在官网下载。下载数据文件|凯斯工程学院 |凯斯西储大学 (case.edu)https://engineering.case.edu/bearingdatacenter/download-data-file
2.VMD函数-matlab代码
VMD函数的matlab代码实现,该代码作为函数实现,无需修改,直接使用即可。
function [u, u_hat, omega] = VMD(signal, alpha, tau, K, DC, init, tol)
% Variational Mode Decomposition% Input and Parameters:
% ---------------------
% signal - the time domain signal (1D) to be decomposed
% alpha - the balancing parameter of the data-fidelity constraint 惩罚因子
% tau - time-step of the dual ascent ( pick 0 for noise-slack )
% K - the number of modes to be recovered 模态分量
% DC - true if the first mode is put and kept at DC (0-freq)
% init - 0 = all omegas start at 0
% 1 = all omegas start uniformly distributed
% 2 = all omegas initialized randomly
% tol - tolerance of convergence criterion; typically around 1e-6
%
% Output:
% -------
% u - the collection of decomposed modes
% u_hat - spectra of the modes
% omega - estimated mode center-frequencies%---------- Preparations
% Period and sampling frequency of input signal
save_T = length(signal);
fs = 1/save_T;% extend the signal by mirroring
T = save_T;
f_mirror(1:T/2) = signal(T/2:-1:1);
f_mirror(T/2+1:3*T/2) = signal;
f_mirror(3*T/2+1:2*T) = signal(T:-1:T/2+1);
f = f_mirror;% Time Domain 0 to T (of mirrored signal)
T = length(f);
t = (1:T)/T;% Spectral Domain discretization
freqs = t-0.5-1/T;% Maximum number of iterations (if not converged yet, then it won't anyway)
N = 500;% For future generalizations: individual alpha for each mode
Alpha = alpha*ones(1,K);% Construct and center f_hat
f_hat = fftshift((fft(f)));
f_hat_plus = f_hat;
f_hat_plus(1:T/2) = 0;% matrix keeping track of every iterant // could be discarded for mem
u_hat_plus = zeros(N, length(freqs), K);% Initialization of omega_k
omega_plus = zeros(N, K);
switch initcase 1for i = 1:Komega_plus(1,i) = (0.5/K)*(i-1);endcase 2omega_plus(1,:) = sort(exp(log(fs) + (log(0.5)-log(fs))*rand(1,K)));otherwiseomega_plus(1,:) = 0;
end% if DC mode imposed, set its omega to 0
if DComega_plus(1,1) = 0;
end% start with empty dual variables
lambda_hat = zeros(N, length(freqs));% other inits
uDiff = tol+eps; % update step
n = 1; % loop counter
sum_uk = 0; % accumulator% ----------- Main loop for iterative updateswhile ( uDiff > tol && n < N ) % not converged and below iterations limit% update first mode accumulatork = 1;sum_uk = u_hat_plus(n,:,K) + sum_uk - u_hat_plus(n,:,1);% update spectrum of first mode through Wiener filter of residualsu_hat_plus(n+1,:,k) = (f_hat_plus - sum_uk - lambda_hat(n,:)/2)./(1+Alpha(1,k)*(freqs - omega_plus(n,k)).^2);% update first omega if not held at 0if ~DComega_plus(n+1,k) = (freqs(T/2+1:T)*(abs(u_hat_plus(n+1, T/2+1:T, k)).^2)')/sum(abs(u_hat_plus(n+1,T/2+1:T,k)).^2);end% update of any other modefor k=2:K% accumulatorsum_uk = u_hat_plus(n+1,:,k-1) + sum_uk - u_hat_plus(n,:,k);% mode spectrumu_hat_plus(n+1,:,k) = (f_hat_plus - sum_uk - lambda_hat(n,:)/2)./(1+Alpha(1,k)*(freqs - omega_plus(n,k)).^2);% center frequenciesomega_plus(n+1,k) = (freqs(T/2+1:T)*(abs(u_hat_plus(n+1, T/2+1:T, k)).^2)')/sum(abs(u_hat_plus(n+1,T/2+1:T,k)).^2);end% Dual ascentlambda_hat(n+1,:) = lambda_hat(n,:) + tau*(sum(u_hat_plus(n+1,:,:),3) - f_hat_plus);% loop countern = n+1;% converged yet?uDiff = eps;for i=1:KuDiff = uDiff + 1/T*(u_hat_plus(n,:,i)-u_hat_plus(n-1,:,i))*conj((u_hat_plus(n,:,i)-u_hat_plus(n-1,:,i)))';enduDiff = abs(uDiff);end%------ Postprocessing and cleanup% discard empty space if converged early
N = min(N,n);
omega = omega_plus(1:N,:);% Signal reconstruction
u_hat = zeros(T, K);
u_hat((T/2+1):T,:) = squeeze(u_hat_plus(N,(T/2+1):T,:));
u_hat((T/2+1):-1:2,:) = squeeze(conj(u_hat_plus(N,(T/2+1):T,:)));
u_hat(1,:) = conj(u_hat(end,:));u = zeros(K,length(t));for k = 1:Ku(k,:)=real(ifft(ifftshift(u_hat(:,k))));
end% remove mirror part
u = u(:,T/4+1:3*T/4);% recompute spectrum
clear u_hat;
for k = 1:Ku_hat(:,k)=fftshift(fft(u(k,:)))';
endend
3.采用matlab脚本导入数据并做VMD分解
该段代码将内圈故障数据导入,并进行了VMD分解。其中得到的u即为分解出来的IMF分量。
clc
clear
fs=12000;%采样频率
Ts=1/fs;%采样周期
L=1500;%采样点数
t=(0:L-1)*Ts;%时间序列
STA=1; %采样起始位置
%----------------导入内圈故障的数据-----------------------------------------
load 105.mat
X = X105_DE_time(1:L)'; %这里可以选取DE(驱动端加速度)、FE(风扇端加速度)、BA(基座加速度),直接更改变量名,挑选一种即可。%--------- some sample parameters forVMD:对于VMD样品参数进行设置---------------
alpha = 2500; % moderate bandwidth constraint:适度的带宽约束/惩罚因子
tau = 0; % noise-tolerance (no strict fidelity enforcement):噪声容限(没有严格的保真度执行)
K = 8; % modes:分解的模态数,可以自行设置,这里以8为例。
DC = 0; % no DC part imposed:无直流部分
init = 1; % initialize omegas uniformly :omegas的均匀初始化
tol = 1e-7;
%--------------- Run actual VMD code:数据进行vmd分解---------------------------
[u, u_hat, omega] = VMD(X, alpha, tau, K, DC, init, tol); %其中u为分解得到的IMF分量
4. VMD分解图
该段代码用于作图,图片可以参考下面的结果展示。
figure(1);
imfn=u;
n=size(imfn,1);
subplot(n+1,1,1);
plot(t,X); %故障信号
ylabel('原始信号','fontsize',12,'fontname','宋体');for n1=1:nsubplot(n+1,1,n1+1);plot(t,u(n1,:));%输出IMF分量,a(:,n)则表示矩阵a的第n列元素,u(n1,:)表示矩阵u的n1行元素ylabel(['IMF' int2str(n1)]);%int2str(i)是将数值i四舍五入后转变成字符,y轴命名
endxlabel('时间\itt/s','fontsize',12,'fontname','宋体');
5.计算中心频率
%----------------------计算中心频率确定分解个数K-----------------------------
average=mean(omega); %对omega求平均值,即为中心频率。
中心频率可以用来确定模态分量K的个数,average即为计算得出的中心频率。因为是要确定分解层数,将K设置不同的值,例如1-9,比较最后一个分量的频率。可以确定K值的依据为:一旦出现相似频率,此时的K值被确定为最佳K值。
6.画包络线
该段代码用来画包络线,图形参考下面的结果展示。
figure(2)
for i = 1:KHy(i,:)= abs(hilbert(u(i,:)));subplot(K,1,i);plot(t,u(i,:),'k',t,Hy(i,:),'r');xlabel('样点'); ylabel('幅值')grid; legend('信号','包络');
end
title('信号和包络');
set(gcf,'color','w');
7. 画包络谱
该段代码用来画包络谱,图形参考下面的结果展示。
画包络谱的时候,切记要将 “nfft” 设置为采样点数的一半。
% 画包络谱
figure('Name','包络谱','Color','white');
nfft=fix(L/2);
for i = 1:Kp=abs(fft(Hy(i,:))); %并fft,得到p,就是包络线的fft---包络谱 p = p/length(p)*2;p = p(1: fix(length(p)/2));subplot(K,1,i);plot((0:nfft-1)/nfft*fs/2,p) %绘制包络谱%xlim([0 600]) %展示包络谱低频段,这句代码可以自己根据情况选择是否注释if i ==1title('包络谱'); xlabel('频率'); ylabel('幅值')elsexlabel('频率'); ylabel('幅值')end
end
set(gcf,'color','w');
8. 计算峭度值
%% 计算峭度值
for i=1:Ka(i)=kurtosis(imfn(i,:));%峭度disp(['IMF',num2str(i),'的峭度值为:',num2str(a(i))])
end
figure
b = bar(a,0.3);
xlabel('模态函数'); ylabel('峭度值')
set(gca,'xtick',1:1:K);
set(gca,'xticklabel',{'IMF1','IMF2','IMF3','IMF4','IMF5','IMF6','IMF7','IMF8'});
xtips1 = b.XEndPoints;
ytips1 = b.YEndPoints; %获取Bar对象的XEndPoints和YEndPoints属性
labels1 = string(b.YData); %获取条形末端坐标
text(xtips1, ytips1, labels1, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom');
%指定垂直、水平对其,让值显示在条形末端居中
9.计算能量和能量熵
%% 能量熵
for i=1:KEimf(i) = sum(imfn(i,:).^2,2);
end
disp(['IMF分量的能量'])
disp(Eimf(1:K))
% 能量熵
E = sum(Eimf);
for j = 1:Kp(j) = Eimf(j)/E;HE(j)=-sum(p(j).*log(p(j)));
end
disp('EMD能量熵=%.4f');
disp(HE(1:K));
10.近似熵代码
近似熵函数,直接复制粘贴,不用改函数。
%近似熵函数的代码
function ApEn = kApproximateEntropy(data, dim, r)
% 计算近似熵(ApEn)
% data - 待分析数据,需要是一维数据
% dim - 模式维度
% r- 阈值大小,一船选择r=0.1-0.25,
data = data(:); %强制转化数据为列方向
N = length(data);%获取数据长度
result = zeros(1,2);%初始化参数
for j = 1:2
m = dim+j-1;
C = zeros(1,N-m+1);%C值初始化,
dataMat = zeros (m, N-m+1);%初始化
for i = 1:m
dataMat(i,:) = data(i:N-m+i); %
end
% counting similar patterns using distance calculation
for i = 1:N-m+1
tempMat = abs(dataMat - repmat (dataMat(:,i),1,N-m+1));
boolMat = any((tempMat>r),1);
C(i)=sum(~boolMat)/(N-m+1);
end
phi(j) = sum(log(C))/(N-m+1);
end
ApEn = phi(1)-phi(2);
end
近似熵调用:对每个IMF分量计算近似熵
for i1=1:Kimf0=imfn(i1,:);%x=imf0';ApEnx = kApproximateEntropy(x, 2, 0.15); %多尺度熵计算ApEn(1,i1)= ApEnx;disp(['IMF',num2str(i1),'的近似熵为:',num2str(ApEn(1,i1))])
end
11.包络熵
包络熵调用:对每个IMF分量计算包络熵
for i = 1:Kxx= abs(hilbert(u(i,:))); %最小包络熵计算公式!xxx = xx/sum(xx);ssum=0;for ii = 1:size(xxx,2)bb = xxx(1,ii)*log(xxx(1,ii));ssum=ssum+bb;endEnen(i,:) = -ssum;%每个IMF分量的包络熵disp(['IMF',num2str(i),'的包络熵为:',num2str(Enen(i,:))])
end
ff = min(Enen);%求取局部最小包络熵,一般用智能优化算法优化VMD,最为最小适应度函数使用
disp(['局部最小包络熵为:',num2str(ff)])
12.频谱图
figure('Name','频谱图','Color','white');
for i = 1:Kp=abs(fft(u(i,:))); %并fft,得到p,就是包络线的fft---包络谱 subplot(K,1,i);plot((0:L-1)*fs/L,p) %绘制包络谱xlim([0 fs/2]) %展示包络谱低频段,这句代码可以自己根据情况选择是否注释if i ==1title('频谱图'); xlabel('频率'); ylabel(['IMF' int2str(i)]);%int2str(i)是将数值i四舍五入后转变成字符,y轴命名elsexlabel('频率'); ylabel(['IMF' int2str(i)]);%int2str(i)是将数值i四舍五入后转变成字符,y轴命名end
end
set(gcf,'color','w');
13.结果展示
VMD分解图
包络线图
包络熵低频段
峭度值
包络熵计算
近似熵
13.智能算法优化VMD参数
智能算法优化VMD参数将在下一篇文章介绍!敬请关注谢谢!
(2条消息) 麻雀算法SSA,优化VMD,适应度函数为最小包络熵,包含MATLAB源代码,直接复制粘贴!_今天吃饺子的博客-CSDN博客
VMD分解,matlab代码,包络线,包络谱,中心频率,峭度值,能量熵,近似熵,包络熵,频谱图,希尔伯特变换,包含所有程序MATLAB代码,-西储大学数据集为例相关推荐
- 【动手教你学故障诊断:Python实现Tensorflow+CNN深度学习的轴承故障诊断(西储大学数据集)(含完整代码)】
项目名称 动手教你学故障诊断:Python实现基于Tensorflow+CNN深度学习的轴承故障诊断(西储大学数据集)(含完整代码) 项目介绍 该项目使用tensorflow和keras搭建深度学习C ...
- 西储大学轴承数据处理--附MATLAB代码
关于西储大学数据介绍的文章,已经有很多很多了.本篇文章旨在采用matlab语言,实现对西储大学数据处理,由原数据直接变成可供大家直接用来故障识别的数据. 首先附上一张西储大学数据集介绍的图片,从图中可 ...
- CWRU(凯斯西储大学轴承数据中心)数据集获取
CWRU(凯斯西储大学轴承数据中心)数据集获取网址 https://engineering.case.edu/bearingdatacenter/apparatus-and-procedures 数据 ...
- 小波变换之(凯斯西储大学CWRU数据转为时频图像)-matlab
凯西储大学(CWRU)数据集下载地址:CWUR完整版 1,挑选出所需数据,每种工况下的10类数据,也就是每个类别一个mat格式的文件 2,加载原始路径,预先定义好所需要的参数 采样频率12Khz ad ...
- 用matlab绘制P三曲线,知道曲线方程 怎么用matlab绘制三维图 一定要给出程序 , matlab怎样画三维曲线...
导航:网站首页 > 知道曲线方程 怎么用matlab绘制三维图 一定要给出程序 , matlab怎样画三维曲线 知道曲线方程 怎么用matlab绘制三维图 一定要给出程序 , matlab怎样画 ...
- emd,eemd,vmd,频谱图,分解图对比matlab代码
作为自己学习的一个记录吧. 对这个信号进行实验,其中公式是截图截的,懒得改了,f1就是s1. 对这个s信号进行分解.下面开始代码操作: 原始信号生成:运行该段代码,生成一个s.mat数据,并作图. c ...
- 基于EMD-HHT包络谱与堆栈降噪自编码SDAE轴承故障诊断
目录 0.引言 1.方法原理 1.1 EMD-HHT包络谱 1.2 PCA 1.3 SDAE 2.流程 3.具体实现 3.1 数据准备 3.2 EMD-HHT包络谱求取 3.3 PCA降维 3.4 故 ...
- 码元速率的盲估计-包络谱法
码元速率的盲估计 实现自论文 主要实现的是包络谱的检测方法. 算法的大致原理请看论文,本文这里是帮助实现了论文代码.分享出来. 源代码 clc clear close all % 这是一个对码元速率的 ...
- 中心频率法确认VMD的K值的python实现
变分模态分解VMD 1.1原理 VMD算法是利用VMD是一种新的信号分解方法,就是求解变分问题.首先构造变分问题,将信号f分解为K个IMF分量,且每个模态分量都有各自的中心频率,同时求解模态的估计带宽 ...
最新文章
- Javascript到PHP加密通讯的简单实现
- OVS DPDK vhost-user详解(十一)
- java der pem_JAVA解析各种编码密钥对(DER、PEM、openssh公钥)
- 联想gen系列服务器,Hpe Microserver Gen10 Plus开箱
- 如何利用ffmpeg提供的API函数进行多媒体文件的解封装
- 给 WordPress 博客添加随机文章的3个方法
- 《Iphone SDK3开发快速上手》
- Mybatis JPA-集成方案+代码解析
- 【IIOT】欧姆龙PLC数采之NX/NJ系列
- 使用UCSC基因组浏览器可视化测序深度分布数据
- 内置平头哥玄铁的WiFi和蓝牙芯片
- sqli-labs闯关指南 1—10
- 2016计算机cpu,2016年12月电脑CPU天梯图一览
- 手动抛出异常回滚事务,且返回数据给前端
- BCK生态链第一条价值共识侧链CK正式上线
- python程序设计与基础教程第六章上机实验_《Python程序设计与算法基础教程》教学大纲.doc...
- IDM怎么下载电影电视剧?快来GET
- 姚班和北大计算机哪个更牛,北大图灵班与清华姚班到底选择哪一个?
- 用java写米与英尺,JAVA 英制变换(英尺与英寸)
- 计算机工具栏使用说明,2345好压主菜单:工具菜单使用说明
热门文章
- Timeline Service V2.0 Reader 启动不了
- Python的特点有哪些?
- 乱码问题之文件,文本文件以及编码
- Artifactory Jfrog与Nexus
- 医疗器械公司计算机系统和功能说明,医疗器械计算机信息管理系统基本情况介绍和功能说明...
- WPF开发:WindowsFormsHost上浮动控件方法
- Cordova iOS 项目中微信/qq/Apple ID插件的安装以及登录的实现
- listbox控件的一些操作
- 什么是VOL版?什么是VLK版
- There are multiple modules with names that only differ in casing. This can lead to unexpected behavi