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(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