Vmd的分解代码

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