✅作者简介：热爱科研的Matlab仿真开发者，擅长数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。
🍎 往期回顾关注个人主页：Matlab科研工作室
👇 关注我领取海量matlab电子书和数学建模资料
🍊个人信条：格物致知,完整Matlab代码获取及仿真咨询内容私信
🔥 内容介绍
在数字信号处理的广阔领域中，信号去噪技术犹如一把神奇的 “降噪神器”，发挥着举足轻重的作用。无论是在音频处理中，让我们听到纯净清晰的音乐和语音；还是在图像处理里，呈现出细腻逼真的画面；亦或是在通信系统中，确保信息准确无误地传输，信号去噪都扮演着不可或缺的角色。
传统的信号去噪方法，如滤波器、小波变换等，在过去的几十年里为信号处理领域做出了重要贡献。滤波器通过对信号的频率成分进行筛选，能够有效地去除特定频率范围内的噪声；小波变换则利用小波函数的多分辨率特性，对信号进行分解和重构，从而实现去噪的目的。然而，这些传统方法也存在着一些局限性。它们往往需要人工设计参数，对信号的先验知识要求较高。一旦信号的特性发生变化，或者噪声的类型较为复杂，传统方法的去噪效果就会大打折扣，甚至可能会对信号本身造成损伤。
为了克服传统去噪方法的不足，科研人员不断探索新的技术和方法。近年来，基于智能优化算法的信号去噪技术逐渐成为研究的热点。其中，基于雪橇犬算法（SDO）优化变分模态分解（VMD）的数字信号去噪方法，以其独特的优势和出色的性能，受到了广泛的关注。
认识变分模态分解（VMD）
在深入探讨基于雪橇犬算法优化变分模态分解的数字信号去噪方法之前，我们先来认识一下变分模态分解（VMD）。
（一）VMD 原理剖析
Image
（二）VMD 关键参数
在 VMD 中，有两个关键参数对分解效果起着至关重要的作用，它们分别是 K 值（IMF 数量）和 alpha 值（惩罚因子）。
K 值决定了信号被分解成的 IMF 数量。如果 K 值设置过小，会导致模态混叠，即不同频率成分被压缩到同一模态中，难以有效分离特征；而当 K 值过大时，又会引入虚假分量并增加计算负担，导致计算冗余。例如，在对机械振动信号进行分析时，如果 K 值设置不当，可能会使原本能够清晰反映故障特征的信号变得模糊不清，无法准确判断设备的运行状态。
alpha 值则控制着惩罚项的权重，决定了对重构误差的惩罚程度，从而影响到各个模态分量之间的带宽约束。当 alpha 较小时，意味着对重构误差的惩罚较轻，这使得 VMD 更加灵活地调整各模态分量的能量分布，可能导致更广泛的频谱覆盖范围以及更高的分辨率，但也可能带来过拟合的风险，即某些噪声成分被误认为是有效信号的一部分；相反，较大的 alpha 值会增加对重构误差的惩罚力度，促使 VMD 尽力保持原始信号的整体结构不变形，这样可以减少不必要的细节波动，提高稳定性并降低过拟合的可能性，但另一方面也可能抑制一些细微特征的表现，造成一定程度上的欠拟合现象。
雪橇犬算法（SDO）探秘
（一）SDO 灵感来源
雪橇犬算法（SDO）是一种新型群智能优化算法，它的灵感源于北极雪橇犬群在极地环境中协同拉雪橇的独特行为模式。在北极地区，雪橇犬是当地居民重要的运输伙伴，它们在冰天雪地中执行各种任务，展现出了卓越的耐力与协作能力。
研究人员通过深入观察雪橇犬群的行为，发现它们在拉雪橇过程中存在着明确的分工。位于队伍最前端的头犬，犹如智慧的领航者，凭借着出色的方向辨识能力和丰富的经验，引领着整个队伍穿越复杂的地形，为团队指明前进的方向；梯队犬则紧密协作，各自发挥力量，齐心协力拉动雪橇，它们相互配合，保持着队伍的稳定行进；而退役犬虽然不再直接参与拉雪橇的工作，但它们的经验却通过某种方式传承给了年轻的雪橇犬，为整个群体的发展贡献着自己的智慧。
基于这些有趣的观察，研究团队巧妙地将雪橇犬的行为特征转化为数学模型，构建了雪橇犬算法。该算法模拟了雪橇犬群的多种行为，包括头犬引导的全局探索、能耗自适应的局部开发以及梯队轮换的多样性保持策略等。通过这些行为的模拟，SDO 算法能够在解空间中进行高效的搜索，寻找最优解。
（二）SDO 核心机制
种群初始化：在 SDO 算法中，首先需要随机生成 N 个候选解，这些候选解就如同雪橇犬个体，它们共同组成了初始种群。每个候选解都对应着问题的一个可能解，而种群规模 N 的设置至关重要。如果种群规模过小，算法可能无法充分探索解空间，容易陷入局部最优解；相反，若种群规模过大，虽然能够增加搜索的多样性，但也会导致计算量大幅增加，降低算法的运行效率。因此，需要根据具体问题的规模和复杂程度，自适应地设置种群规模，以确保种群多样性与搜索效率之间的平衡。
适应度评估：为了衡量每个候选解的优劣，SDO 算法以优化目标为基础构建了适应度函数。在不同的应用场景中，优化目标可能各不相同。例如，在流水车间调度问题中，优化目标可能是最小化最大完工时间（Makespan）；在信号去噪问题中，优化目标则可能是提高信号的信噪比或降低噪声的影响。通过适应度函数，算法可以计算出每个候选解的适应度值，适应度值越高（或越低，取决于优化目标是最大化还是最小化），表明该候选解越接近最优解。
核心搜索阶段：
Image
SDO 优化 VMD 的神奇之旅
（一）优化思路
在信号去噪的征程中，如何让 VMD 的 K 值和 alpha 值达到最优，成为了关键的挑战。而雪橇犬算法（SDO）的出现，为解决这一难题带来了新的曙光。SDO 以其强大的全局搜索能力而著称，它能够在广阔的解空间中高效地寻找最优解。就像一群训练有素的雪橇犬，在茫茫的雪地里精准地找到前进的道路。利用 SDO 的这一优势，我们可以将其应用于 VMD 的参数优化中，通过不断地搜索和迭代，寻找出使 VMD 分解效果最佳的 K 值和 alpha 值，从而提高信号的分解和去噪效果。
（二）实现步骤
定义目标函数：在 SDO 优化 VMD 的过程中，目标函数的定义至关重要。这里，我们选择包络熵作为适应度函数。包络熵是一种用于衡量信号复杂度和不确定性的指标，它能够反映信号的稀疏性。具体来说，包络熵越小，意味着信号的复杂度越低，稀疏性越好，也就表示信号中包含的有效信息越多，噪声越少。例如，在一个理想的无噪声信号中，包络熵的值会非常小，因为信号的变化是规律且可预测的；而当信号中存在大量噪声时，包络熵会显著增大，因为噪声的存在使得信号变得更加复杂和无序。通过将包络熵作为适应度函数，SDO 可以在迭代过程中，根据包络熵的大小来评估不同 K 值和 alpha 值组合下的 VMD 分解效果，从而引导算法朝着使包络熵最小的方向搜索，即寻找最优的 K 值和 alpha 值，以实现最佳的信号去噪效果。
参数范围设定：合理设定 K 值和 alpha 值的范围是 SDO 优化 VMD 的重要前提。对于 K 值，其取值范围通常需要根据信号的特性和经验来确定。一般来说，K 值的下限可以设置为 2，这是因为当 K 值为 1 时，VMD 实际上并没有对信号进行分解，而当 K 值为 2 时，至少可以将信号分解为两个模态，开始体现 VMD 的分解作用。K 值的上限则需要考虑信号的复杂程度和计算资源的限制。如果 K 值过大，不仅会增加计算量，还可能导致过分解，产生虚假模态。例如，对于一些简单的信号，K 值设置为 5 - 10 可能就足够了；而对于复杂的多模态信号，K 值可能需要设置在 10 - 20 之间。
对于 alpha 值，其取值范围也需要谨慎选择。alpha 值的下限可以设置为一个较小的值，如 100，这是因为当 alpha 值过小时，VMD 对重构误差的惩罚较轻，可能会导致模态混叠，无法有效地分离信号的不同频率成分。alpha 值的上限则可以设置为一个较大的值，如 10000，当 alpha 值过大时，虽然能够避免模态混叠，但可能会造成局部信息丢失，使分解结果过于平滑，无法准确反映信号的细节特征。在实际应用中，可以通过多次试验和分析，结合信号的具体特点，来确定最合适的 K 值和 alpha 值范围。
- 迭代优化过程：在迭代优化过程中，SDO 首先初始化种群，生成一组随机的 K 值和 alpha 值组合，这些组合就像是雪橇犬在雪地里的不同出发点。然后，对于每个组合，算法会将其代入 VMD 中进行信号分解，并计算分解后各模态分量的包络熵，以此作为适应度值。适应度值越小，说明该组合下的 VMD 分解效果越好，对应的 K 值和 alpha 值也就越优。
接下来，SDO 根据适应度值对种群中的个体进行分类，选出适应度最优的个体作为引导犬，适应度中等的个体作为搜索犬，剩余的个体作为协作犬。引导犬根据自身位置和全局最优解的信息，按照一定的公式更新自己的位置，从而带动整个种群向全局最优解的方向进化；搜索犬在引导犬导航的区域周边进行局部探索，通过引入一定的随机性，避免种群陷入局部最优解；协作犬则跟随引导犬与搜索犬的搜索方向，维持种群的协作性与多样性。在每一次迭代中，算法都会更新种群中每个个体的位置，即调整 K 值和 alpha 值，并重新计算适应度值。不断重复这个过程，直到达到最大迭代次数或适应度值不再明显变化为止。此时，算法输出的全局最优解就是经过优化的 K 值和 alpha 值，将其应用于 VMD 中，即可实现对信号的高效去噪。
⛳️ 运行结果
Image
Image
Image
Image
📣 部分代码
function [u, u_hat, omega] = VMD(signal, alpha, tau, K, DC, init, tol)
%
%
% 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
%
% When using this code, please do cite our paper:
% -----------------------------------------------
% K. Dragomiretskiy, D. Zosso, Variational Mode Decomposition, IEEE Trans.
% on Signal Processing (in press)
% please check here for update reference:
% http://dx.doi.org/10.1109/TSP.2013.2288675
%---------- 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(3T/2+1:2T) = 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 init
case 1
for i = 1:K
omega_plus(1,i) = (0.5/K)*(i-1);
end
case 2
omega_plus(1,:) = sort(exp(log(fs) + (log(0.5)-log(fs))*rand(1,K)));
otherwise
omega_plus(1,:) = 0;
end
% if DC mode imposed, set its omega to 0
if DC
omega_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 updates
while ( uDiff > tol && n < N ) % not converged and below iterations limit
% update first mode accumulator
k = 1;
sum_uk = u_hat_plus(n,:,K) + sum_uk - u_hat_plus(n,:,1);
% update spectrum of first mode through Wiener filter of residuals
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);
% update first omega if not held at 0
if ~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);
end
% update of any other mode
for k=2:K
% accumulator
sum_uk = u_hat_plus(n+1,:,k-1) + sum_uk - u_hat_plus(n,:,k);
% mode spectrum
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);
% center frequencies
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);
end
% Dual ascent
lambda_hat(n+1,:) = lambda_hat(n,:) + tau*(sum(u_hat_plus(n+1,:,:),3) - f_hat_plus);
% loop counter
n = n+1;
% converged yet?
uDiff = eps;
for 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,:,i)))';
end
uDiff = 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:K
u(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:K
u_hat(:,k)=fftshift(fft(u(k,:)))';
end
end
🔗 参考文献
🎈 部分理论引用网络文献，若有侵权联系博主删除
🏆团队擅长辅导定制多种科研领域MATLAB仿真，助力科研梦：
