VMD算法的MATLAB实现

%主程序

%--------------- Preparation

clear all;

close all;

clc;

% Time Domain 0 to T

T = 1000;

fs = 1/T;

t = (1:T)/T;

freqs = 2*pi*(t-0.5-1/T)/(fs);

% center frequencies of components

f_1 = 2;

f_2 = 24;

f_3 = 288;

% modes

v_1 = (cos(2*pi*f_1*t));

v_2 = 1/4*(cos(2*pi*f_2*t));

v_3 = 1/16*(cos(2*pi*f_3*t));

% for visualization purposes

wsub{1} = 2*pi*f_1;

wsub{2} = 2*pi*f_2;

wsub{3} = 2*pi*f_3;

% composite signal, including noise

f = v_1 + v_2 + v_3 + 0.1*randn(size(v_1));

% some sample parameters for VMD

alpha = 2000;        % moderate bandwidth constraint

tau = 0;            % noise-tolerance (no strict fidelity enforcement)

K = 4;              % 4 modes

DC = 0;             % no DC part imposed

init = 1;           % initialize omegas uniformly

tol = 1e-7;


%--------------- Run actual VMD code

[u, u_hat, omega] = VMD(f, alpha, tau, K, DC, init, tol);

subplot(size(u,1)+1,2,1);

plot(t,f,'k');grid on;

title('VMD分解');

subplot(size(u,1)+1,2,2);

plot(freqs,abs(fft(f)),'k');grid on;

title('对应频谱');

for i = 2:size(u,1)+1

    subplot(size(u,1)+1,2,i*2-1);

    plot(t,u(i-1,:),'k');grid on;

    subplot(size(u,1)+1,2,i*2);

    plot(freqs,abs(fft(u(i-1,:))),'k');grid on;

end


%---------------run EMD code

imf = emd(f);

figure;

subplot(size(imf,1)+1,2,1);

plot(t,f,'k');grid on;

title('EMD分解');

subplot(size(imf,1)+1,2,2);

plot(freqs,abs(fft(f)),'k');grid on;

title('对应频谱');

for i = 2:size(imf,1)+1

    subplot(size(imf,1)+1,2,i*2-1);

    plot(t,imf(i-1,:),'k');grid on;

    subplot(size(imf,1)+1,2,i*2);

    plot(freqs,abs(fft(imf(i-1,:))),'k');grid on;

end




%%%%函数部分

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

扫码关注我的公众号，一起学习，一起成长！