VMD（Variational Mode Decomposition）变模式分解

Halcom

VMD（Variational Mode Decomposition）变模式分解
软件：MATLAB

1. function [u, u_hat, omega] = VMD(signal, alpha, tau, K, DC, init, tol)
   
2. % Variational Mode Decomposition
   
3. % Input and Parameters:
   
4. % ---------------------
   
5. % signal  - the time domain signal (1D) to be decomposed
   
6. % alpha   - the balancing parameter of the data-fidelity constraint
   
7. % tau     - time-step of the dual ascent ( pick 0 for noise-slack )
   
8. % K       - the number of modes to be recovered
   
9. % DC      - true if the first mode is put and kept at DC (0-freq)
   
10. % init    - 0 = all omegas start at 0
    
11. %                    1 = all omegas start uniformly distributed
    
12. %                    2 = all omegas initialized randomly
    
13. % tol     - tolerance of convergence criterion; typically around 1e-6
    
14. %
    
15. % Output:
    
16. % -------
    
17. % u       - the collection of decomposed modes
    
18. % u_hat   - spectra of the modes
    
19. % omega   - estimated mode center-frequencies
    
20. %
    
21. % When using this code, please do cite our paper:
    
22. % K. Dragomiretskiy, D. Zosso, Variational Mode Decomposition, IEEE Trans. on Signal Processing (in press)
    
23. % please check here for update reference: http://dx.doi.org/10.1109/TSP.2013.2288675
    
24. 
    
25. %---------- Preparations
    
26. % Period and sampling frequency of input signal
    
27. save_T = length(signal);
    
28. fs = 1/save_T;
    
29. 
    
30. % extend the signal by mirroring
    
31. T = save_T;
    
32. f_mirror(1:T/2) = signal(T/2:-1:1);
    
33. f_mirror(T/2+1:3*T/2) = signal;
    
34. f_mirror(3*T/2+1:2*T) = signal(T:-1:T/2+1);
    
35. f = f_mirror;
    
36. 
    
37. % Time Domain 0 to T (of mirrored signal)
    
38. T = length(f);
    
39. t = (1:T)/T;
    
40. 
    
41. % Spectral Domain discretization
    
42. freqs = t-0.5-1/T;
    
43. 
    
44. % Maximum number of iterations (if not converged yet, then it won't anyway)
    
45. N = 500;
    
46. 
    
47. % For future generalizations: individual alpha for each mode
    
48. Alpha = alpha*ones(1,K);
    
49. 
    
50. % Construct and center f_hat
    
51. f_hat = fftshift((fft(f)));
    
52. f_hat_plus = f_hat;
    
53. f_hat_plus(1:T/2) = 0;
    
54. 
    
55. % matrix keeping track of every iterant // could be discarded for mem
    
56. u_hat_plus = zeros(N, length(freqs), K);
    
57. 
    
58. % Initialization of omega_k
    
59. omega_plus = zeros(N, K);
    
60. switch init
    
61.     case 1
    
62.         for i = 1:K
    
63.             omega_plus(1,i) = (0.5/K)*(i-1);
    
64.         end
    
65.     case 2
    
66.         omega_plus(1,:) = sort(exp(log(fs) + (log(0.5)-log(fs))*rand(1,K)));
    
67.     otherwise
    
68.         omega_plus(1,:) = 0;
    
69. end
    
70. 
    
71. % if DC mode imposed, set its omega to 0
    
72. if DC
    
73.     omega_plus(1,1) = 0;
    
74. end
    
75. 
    
76. % start with empty dual variables
    
77. lambda_hat = zeros(N, length(freqs));
    
78. 
    
79. % other inits
    
80. uDiff = tol+eps; % update step
    
81. n = 1;      % loop counter
    
82. sum_uk = 0; % accumulator
    
83. 
    
84. % ----------- Main loop for iterative updates
    
85. while ( uDiff > tol &&  n < N ) % not converged and below iterations limit
    
86.    
    
87.     % update first mode accumulator
    
88.     k = 1;
    
89.     sum_uk = u_hat_plus(n,:,K) + sum_uk - u_hat_plus(n,:,1);
    
90.    
    
91.     % update spectrum of first mode through Wiener filter of residuals
    
92.     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);
    
93.    
    
94.     % update first omega if not held at 0
    
95.     if ~DC
    
96.         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);
    
97.     end
    
98.    
    
99.     % update of any other mode
    
100.     for k=2:K
     
101.         
     
102.         % accumulator
     
103.         sum_uk = u_hat_plus(n+1,:,k-1) + sum_uk - u_hat_plus(n,:,k);
     
104.         
     
105.         % mode spectrum
     
106.         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);
     
107.         
     
108.         % center frequencies
     
109.         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);
     
110.         
     
111.     end
     
112.    
     
113.     % Dual ascent
     
114.     lambda_hat(n+1,:) = lambda_hat(n,:) + tau*(sum(u_hat_plus(n+1,:,:),3) - f_hat_plus);
     
115.    
     
116.     % loop counter
     
117.     n = n+1;
     
118.    
     
119.     % converged yet?
     
120.     uDiff = eps;
     
121.     for i=1:K
     
122.         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)))';
     
123.     end
     
124.     uDiff = abs(uDiff);
     
125.    
     
126. end
     
127. 
     
128. %------ Postprocessing and cleanup
     
129. 
     
130. % discard empty space if converged early
     
131. N = min(N,n);
     
132. omega = omega_plus(1:N,:);
     
133. 
     
134. % Signal reconstruction
     
135. u_hat = zeros(T, K);
     
136. u_hat((T/2+1):T,:) = squeeze(u_hat_plus(N,(T/2+1):T,:));
     
137. u_hat((T/2+1):-1:2,:) = squeeze(conj(u_hat_plus(N,(T/2+1):T,:)));
     
138. u_hat(1,:) = conj(u_hat(end,:));
     
139. 
     
140. u = zeros(K,length(t));
     
141. 
     
142. for k = 1:K
     
143.     u(k,:)=real(ifft(ifftshift(u_hat(:,k))));
     
144. end
     
145. 
     
146. % remove mirror part
     
147. u = u(:,T/4+1:3*T/4);
     
148. 
     
149. % recompute spectrum
     
150. clear u_hat;
     
151. for k = 1:K
     
152.     u_hat(:,k)=fftshift(fft(u(k,:)))';
     
153. end

复制代码

参考链接：
【1】K. Dragomiretskiy, D. Zosso, Variational Mode Decomposition, IEEE Trans. on Signal Processing.