【信号去噪】基于麻雀算法优化VMD信号去噪(目标函数为包络熵局部极小值)附Matlab代码

✅作者简介：热爱科研的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]. 重庆三峡学院。

⛳️ 代码获取关注我

❤️部分理论引用网络文献，若有侵权联系博主删除

❤️ 关注我领取海量matlab电子书和数学建模资料