VMD分解的matlab程序

基于matlab的VMD分解程序，直接可以调用，效果比EMD分解和小波包分解的效果更好

function [u u_hat omega] = VMD（signal alpha tau K DC init tol）
% Variational Mode Decomposition
% Authors: Konstantin Dragomiretskiy and Dominique Zosso
% zosso@math.ucla.edu --- http://www.math.ucla.edu/~zosso
% Initial release 2013-12-12 （c） 2013
%
% 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（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（1K）;

% 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（1i） = （0.5/K）*（i-1）;
        end
    case 2
        omega_plus（1:） = sort（exp（log（fs） + （log（0.5）-log（fs））*rand（1K）））;
    otherwise
        omega_plus（1:） = 0;
end

% if DC mode imposed set its omega to 0
if DC
    omega_plus（11） = 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 - la

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       4645  2018-03-23 10:18  VMD.m

----------- ---------  ---------- -----  ----

                 4645                    1