✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。
🍎个人主页：Matlab科研工作室
🍊个人信条：格物致知。
更多Matlab仿真内容点击👇
⛄ 内容介绍
然后计算每个麻雀位置下的平均包网络密度;更新最小平均包网络密度,即获得当前群体中最佳个人体的位置更新前麻雀个体的空间位置;输出最佳麻雀个体的位置方向,即得VMD的解析参数组合;根据所得到的分析参数组合对信号进行VMD分析,将分出IMF分量相加到重组信号,即得去后的信号.本发明解决了现在有技术存在的原始VMD算法的解析参数需要根据经验进行人员确定，导致无法得到最佳解析结果，从而影响信号的去有效的问题。以西储大学数据集为例，选用105.mat中的X105_BA_time.mat数据。首先进行VMD分解，采用麻雀优化算法（SSA）对VMD的两个关键参数（惩罚因子α和模态分解数K）进行优化，以最小包络熵为适应度值。其他智能优化算法同样适用，关键要学会最小包络熵代码的编写，实验过程中，会实时显示每次寻优后的最小包络熵值和VMD对应的两个最佳参数。本次寻优共100次（自己可以随意更改寻优次数）。
函数 [u, u_hat, omega] = VMD(signal, alpha, tau, K, DC, init, tol)
% 变分模式分解
% 作者：Konstantin Dragomiretskiy 和 Dominique Zosso
% zosso@math.ucla.edu --- http://www.math.ucla.edu/~zosso
% 初始版本 2013-12-12 (c) 2013
%
% 输入和参数：
% ----------------------
% signal - 要分解的时域信号（1D）
% alpha - 数据保真度约束的平衡参数
% tau - 双上升的时间步长（选择 0 表示噪声松弛）
% K - 要恢复的模式数
% DC - 如果第一个模式被放置并保持在 DC (0-freq)
% init - 0 = 所有 omegas 从 0 开始
% 1 = 所有 omegas 开始均匀分布
% 2 = 所有 omegas 随机初始化
% tol - 收敛标准的公差；通常在 1e-6 左右
%
％ 输出：
% ------
% u - 分解模式的集合
% u_hat - 模式的光谱
% omega - 估计模式中心频率
%
% 使用此代码时，请务必引用我们的论文：
% ----------------------------------------------
% K. Dragomiretskiy, D. Zosso, 变分模式分解, IEEE Trans.
信号处理百分比（印刷中）
% 请在此处查看更新参考：
% http://dx.doi.org/10.1109/TSP.2013.2288675
%------------准备工作
% 输入信号的周期和采样频率
save_T = 长度（信号）；
fs = 1/保存_T；
% 通过镜像扩展信号
T = 保存_T；
f_mirror(1:T/2) = 信号(T/2:-1:1);
f_mirror(T/2+1:3*T/2) = 信号；
f_mirror(3*T/2+1:2*T) = 信号(T:-1:T/2+1);
f = f_mirror;
% 时域 0 到 T（镜像信号）
T = 长度（f）；
t = (1:T)/T;
% 频谱域离散化
频率 = t-0.5-1/T；
% 最大迭代次数（如果还没有收敛，那么无论如何都不会）
N = 500;
% 对于未来的概括：每个模式的单独 alpha
Alpha = alpha*ones(1,K);
% 构建并居中 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 // 可以为 mem 丢弃
u_hat_plus = zeros(N, length(freqs), K);
% omega_k 的初始化
omega_plus = zeros(N, K);
切换初始化
情况1
因为我 = 1:K
omega_plus(1,i) = (0.5/K)*(i-1);
结尾
案例2
omega_plus(1,:) = sort(exp(log(fs) + (log(0.5)-log(fs))*rand(1,K)));
否则
omega_plus(1,:) = 0;
结尾
% 如果施加 DC 模式，将其 omega 设置为 0
如果直流
omega_plus(1,1) = 0;
结尾
% 从空对偶变量开始
lambda_hat = zeros(N, 长度(频率));
% 其他初始化
uDiff = tol+eps；% 更新步骤
n = 1; %循环计数器
sum_uk = 0; % 累加器
% ---------- 迭代更新的主循环
while ( uDiff > tol && n < N ) % 未收敛且低于迭代限制
% 更新第一模式累加器
k = 1;
sum_uk = u_hat_plus(n,:,K) + sum_uk - u_hat_plus(n,:,1);
% 通过残差的维纳滤波器更新第一模式的频谱
u_hat_plus(n+1,:,k) = (f_hat_plus - sum_uk - lambda_hat(n,:)/2)./(1+Alpha(1,k)*(freqs - omega_plus(n,k)).^2 );
% 如果不保持在 0，则更新第一个 omega
如果~DC
omega_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);
结尾
% 任何其他模式的更新
对于 k=2:K
% 累加器
sum_uk = u_hat_plus(n+1,:,k-1) + sum_uk - u_hat_plus(n,:,k);
%模式谱
u_hat_plus(n+1,:,k) = (f_hat_plus - sum_uk - lambda_hat(n,:)/2)./(1+Alpha(1,k)*(freqs - omega_plus(n,k)).^2 );
%中心频率
omega_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);
结尾
% 双重上升
lambda_hat(n+1,:) = lambda_hat(n,:) + tau*(总和(u_hat_plus(n+1,:,:),3) - f_hat_plus);
%循环计数器
n = n + 1;
%收敛了吗？
uDiff = eps；
对于 i=1:K
uDiff = 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,:,我）））';
结尾
uDiff = abs(uDiff);
结尾
%------ 后处理和清理
% 如果早期收敛则丢弃空白空间
N = 最小值（N，n）；
omega = omega_plus(1:N,:);
% 信号重建
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));
对于 k = 1:K
u(k,:)=real(ifft(ifftshift(u_hat(:,k))));
结尾
% 删除镜像部分
u = u(:,T/4+1:3*T/4);
% 重新计算频谱
清除你的帽子；
对于 k = 1:K
u_hat(:,k)=fftshift(fft(u(k,:)))';
结尾
结尾
⛄ 运行结果
⛄ 参考文采
[1]孟博, 余中舟, 侯战, 等. 基于VMD-SSA-LSSVM的齿轮箱故障检测模型:, CN114964783A[P]. 2022.
[2]逆全波。基于改进VMD的滚动轴承故障无线监测系统[D]. 重庆三峡学院。
