请选择 进入手机版 | 继续访问电脑版

Hello Mat

 找回密码
 立即注册
查看: 4001|回复: 0

VMD(Variational Mode Decomposition)变模式分解

[复制链接]

1064

主题

1239

帖子

13

金钱

管理员

Rank: 9Rank: 9Rank: 9

积分
21393
发表于 2017-10-10 21:34:58 | 显示全部楼层 |阅读模式
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. %---------- Preparations
  25. % Period and sampling frequency of input signal
  26. save_T = length(signal);
  27. fs = 1/save_T;

  28. % extend the signal by mirroring
  29. T = save_T;
  30. f_mirror(1:T/2) = signal(T/2:-1:1);
  31. f_mirror(T/2+1:3*T/2) = signal;
  32. f_mirror(3*T/2+1:2*T) = signal(T:-1:T/2+1);
  33. f = f_mirror;

  34. % Time Domain 0 to T (of mirrored signal)
  35. T = length(f);
  36. t = (1:T)/T;

  37. % Spectral Domain discretization
  38. freqs = t-0.5-1/T;

  39. % Maximum number of iterations (if not converged yet, then it won't anyway)
  40. N = 500;

  41. % For future generalizations: individual alpha for each mode
  42. Alpha = alpha*ones(1,K);

  43. % Construct and center f_hat
  44. f_hat = fftshift((fft(f)));
  45. f_hat_plus = f_hat;
  46. f_hat_plus(1:T/2) = 0;

  47. % matrix keeping track of every iterant // could be discarded for mem
  48. u_hat_plus = zeros(N, length(freqs), K);

  49. % Initialization of omega_k
  50. omega_plus = zeros(N, K);
  51. switch init
  52.     case 1
  53.         for i = 1:K
  54.             omega_plus(1,i) = (0.5/K)*(i-1);
  55.         end
  56.     case 2
  57.         omega_plus(1,:) = sort(exp(log(fs) + (log(0.5)-log(fs))*rand(1,K)));
  58.     otherwise
  59.         omega_plus(1,:) = 0;
  60. end

  61. % if DC mode imposed, set its omega to 0
  62. if DC
  63.     omega_plus(1,1) = 0;
  64. end

  65. % start with empty dual variables
  66. lambda_hat = zeros(N, length(freqs));

  67. % other inits
  68. uDiff = tol+eps; % update step
  69. n = 1;      % loop counter
  70. sum_uk = 0; % accumulator

  71. % ----------- Main loop for iterative updates
  72. while ( uDiff > tol &&  n < N ) % not converged and below iterations limit
  73.    
  74.     % update first mode accumulator
  75.     k = 1;
  76.     sum_uk = u_hat_plus(n,:,K) + sum_uk - u_hat_plus(n,:,1);
  77.    
  78.     % update spectrum of first mode through Wiener filter of residuals
  79.     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);
  80.    
  81.     % update first omega if not held at 0
  82.     if ~DC
  83.         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);
  84.     end
  85.    
  86.     % update of any other mode
  87.     for k=2:K
  88.         
  89.         % accumulator
  90.         sum_uk = u_hat_plus(n+1,:,k-1) + sum_uk - u_hat_plus(n,:,k);
  91.         
  92.         % mode spectrum
  93.         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);
  94.         
  95.         % center frequencies
  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.         
  98.     end
  99.    
  100.     % Dual ascent
  101.     lambda_hat(n+1,:) = lambda_hat(n,:) + tau*(sum(u_hat_plus(n+1,:,:),3) - f_hat_plus);
  102.    
  103.     % loop counter
  104.     n = n+1;
  105.    
  106.     % converged yet?
  107.     uDiff = eps;
  108.     for i=1:K
  109.         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)))';
  110.     end
  111.     uDiff = abs(uDiff);
  112.    
  113. end

  114. %------ Postprocessing and cleanup

  115. % discard empty space if converged early
  116. N = min(N,n);
  117. omega = omega_plus(1:N,:);

  118. % Signal reconstruction
  119. u_hat = zeros(T, K);
  120. u_hat((T/2+1):T,:) = squeeze(u_hat_plus(N,(T/2+1):T,:));
  121. u_hat((T/2+1):-1:2,:) = squeeze(conj(u_hat_plus(N,(T/2+1):T,:)));
  122. u_hat(1,:) = conj(u_hat(end,:));

  123. u = zeros(K,length(t));

  124. for k = 1:K
  125.     u(k,:)=real(ifft(ifftshift(u_hat(:,k))));
  126. end

  127. % remove mirror part
  128. u = u(:,T/4+1:3*T/4);

  129. % recompute spectrum
  130. clear u_hat;
  131. for k = 1:K
  132.     u_hat(:,k)=fftshift(fft(u(k,:)))';
  133. end
复制代码

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







算法QQ  3283892722
群智能算法链接http://halcom.cn/forum.php?mod=forumdisplay&fid=73
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Python|Opencv|MATLAB|Halcom.cn

GMT+8, 2022-7-3 14:14 , Processed in 0.184113 second(s), 24 queries .

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表