时序分解 | Matlab基于WOA-MVMD鲸鱼算法优化多元变分模态分解

时序分解 | Matlab基于WOA-MVMD鲸鱼算法优化多元变分模态分解

效果一览

基本介绍

WOA-MVMD鲸鱼算法优化多元变分模态分解时间序列信号分解 可直接运行 分解效果好 适合作为创新点（Matlab完整源码和数据)，以包络熵为适应度函数。

1.利用鲸鱼优化算法优化参数k、alpha，分解效果好，包含边际谱、频率图、收敛曲线等图，满足您的需求，使用者较少，适合作为创新点。

2.采用西储大学数据集，运行主程序main即可。

3.包含WOA-MVMD迭代曲线图、MVMD分解图、IMF频域图、包含Hilbert(2D）边际谱图、包含Hilbert(3D）边际谱图。

程序设计

完整源码私信回复Matlab基于WOA-MVMD鲸鱼算法优化多元变分模态分解

nospace
%% 适应度函数
function [ff, min_entropy] = CostMVMD(c, X)
    alpha = c(1);   % 平滑参数
    K = round(c(2));% 模态数
    tau = 0;        % 拉格朗日乘子法步长
    DC = 0;         % 是否提取直流分量
    init = 1;       % 初始中心频率选择方法
    tol = 1e-7;     % 收敛阈值

    % 调用MVMD函数
    [u, ~, ~] = MVMD(X, alpha, tau, K, DC, init, tol);

    % 计算每个模态的包络熵
    fitness = zeros(1, K);
    for i = 1:K
        xx = abs(hilbert(u(i, :))); % 对IMF分量进行希尔伯特变换并求幅值
        xxx = xx / sum(xx); % 归一化
        ssum = 0;

        for ii = 1:size(xxx, 2)
            bb = xxx(1, ii) * log(xxx(1, ii)); % 计算包络熵
            ssum = ssum + bb;  % 求和
        end

        fitness(i) = -ssum;   % 加负号以使得最小化包络熵
    end

    ff = min(fitness); % 返回最小适应度值
    min_entropy = ff;  % 输出最终的最小包络熵值
end
%% MVMD函数
function [u, u_hat, omega] = MVMD(signal, alpha, tau, K, DC, init, tol)
[x, y] = size(signal);
if x > y
	C = y;% number of channels
    T = x;% length of the Signal
	signal = signal';
else
	C = x;% number of channels
    T = y;% length of the Signal
end
%---------- Preparations
% Sampling Frequency
fs = 1/T;
% Mirroring
f(:,1:T/2) = signal(:,T/2:-1:1);
f(:,T/2+1:3*T/2) = signal;
f(:,3*T/2+1:2*T) = signal(:,T:-1:T/2+1);
% Time Domain 0 to T (of mirrored signal)
T = size(f,2);
t = (1:T)/T;
% frequencies
freqs = t-0.5-1/T;
% Construct and center f_hat
f_hat = fftshift(fft(f,[],2),2);
f_hat_plus = f_hat;
f_hat_plus(:,1:T/2) = 0;
%------------ Initialization
% Maximum number of iterations 
N = 500;
% For future generalizations: individual alpha for each mode
Alpha = alpha*ones(1,K);
% matrix keeping track of every iterant 
u_hat_plus_00 = zeros(length(freqs), C, K);
u_hat_plus = zeros(length(freqs), C, K);
omega_plus = zeros(N, K);
% initialize omegas uniformly
switch init
	case 1
        omega_plus(1,:) = (0.5/K)*((1:K)-1);
    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(length(freqs), C, N); 
% other inits
uDiff = tol+eps; % update step
n = 1; % loop counter
sum_uk = zeros(length(freqs), C); % accumulator
%--------------- Algorithm of MVMD
while ( uDiff > tol &&  n < N ) % not converged and below iterations limit	
	% update modes
	for k = 1:K
        % update mode accumulator
        if k > 1
            sum_uk = u_hat_plus(:,:,k-1) + sum_uk - u_hat_plus_00(:,:,k);
        else
            sum_uk = u_hat_plus_00(:,:,K) + sum_uk - u_hat_plus_00(:,:,k);
        end
        % update spectrum of mode through Wiener filter of residuals
		for c = 1:C
			u_hat_plus(:,c,k) = (f_hat_plus(c,:).' - sum_uk(:,c) - lambda_hat(:,c,n)/2)./(1+Alpha(1,k)*(freqs.' - omega_plus(n,k)).^2);
		end
		% update first omega if not held at 0
        if DC || (k > 1)
            % center frequencies
            numerator = freqs(T/2+1:T)*(abs(u_hat_plus(T/2+1:T,:, k)).^2);
            denominator = sum(abs(u_hat_plus(T/2+1:T,:,k)).^2);
            temp1 = sum(numerator);
            temp2 = sum(denominator);
            omega_plus(n+1,k) = temp1/temp2;
        end
	end
	% Dual ascent
    lambda_hat(:,:,n+1) = lambda_hat(:,:,n) + tau*(sum(u_hat_plus,3) - f_hat_plus.');
	% loop counter
	n = n+1;
    u_hat_plus_m1 = u_hat_plus_00;
    u_hat_plus_00 = u_hat_plus;
	% converged yet?
    uDiff = u_hat_plus_00 - u_hat_plus_m1;
    uDiff = 1/T*(uDiff).*conj(uDiff);
	uDiff = eps+abs(sum(uDiff(:)));
end
%------ Post-processing and cleanup
% discard empty space if converged early
N = min(N,n);
omega = omega_plus(1:N,:);
% Signal reconstruction
u_hat = zeros(T, K, C);
for c = 1:C
	u_hat((T/2+1):T,:,c) = squeeze(u_hat_plus((T/2+1):T,c,:));
	u_hat((T/2+1):-1:2,:,c) = squeeze(conj(u_hat_plus((T/2+1):T,c,:)));
	u_hat(1,:,c) = conj(u_hat(end,:,c));
end
u = zeros(K,length(t),C);
for k = 1:K
	for c = 1:C
		u(k,:,c)=real(ifft(ifftshift(u_hat(:,k,c))));
	end
end
% remove mirror part
u = u(:,T/4+1:3*T/4,:);
% recompute spectrum
clear u_hat;
for k = 1:K
	for c = 1:C
		u_hat(:,k,c)=fftshift(fft(u(k,:,c)))';
	end
end
u_hat = permute(u_hat, [2 1 3]);
end

参考资料

[1] https://blog.csdn.net/kjm13182345320/article/details/129215161
[2] https://blog.csdn.net/kjm13182345320/article/details/128105718