function [u, u_hat, omega] = VMD(signal, alpha, tau, K, DC, init, tol)
% Variational Mode Decomposition
% 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):如果第一种模式被放置并保持在DC(0频率),则为真
% 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(收敛公差准则;通常大约e-6)
%
% Output:
% -------
% u - the collection of decomposed modes(分解模态的集合)
% u_hat - spectra of the modes(模谱)
% omega - estimated mode center-frequencies(中心频率的估计模态)
% abs是取绝对值函数
% ones是取矩阵函数,行列为1
% sort是升序排序函数
% squeeze是去除矩阵一维行列的函数
% conj函数计算复数的共轭值
% Period and sampling frequency of input signal(输入信号的周期和采样频率)
save_T = length(signal);
fs = 1/save_T;%分解算法中的采样频率和时间是标准化的,分解信号的采样时间为1s,然后就得到相应的采样频率。采点:length(signal),频率:1/ length(signal)。
% 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;%freqs 进行移位是由于进行傅里叶变换时,会有正负对称的频率,分析时一般只有正频率,所以看到的频谱图是没有负频率的。
% 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
变分模态分解_变分模态分解_vmd_
5星 · 超过95%的资源 79 浏览量
2021-09-29
14:59:01
上传
评论 4
收藏 4KB ZIP 举报
心若悬河
- 粉丝: 52
- 资源: 3958
最新资源
- 前端开发-什么是前端开发-关于前端开发的一些相关介绍
- Sora AI-关于文生视频的使用场景说明
- suno AI文生视频的相关教程和介绍使用
- 什么是后端开发-关于后端开发的一些小介绍分享
- Jurassic Pack Vol. II Dinosaurs 侏罗纪包卷恐龙二号Unity游戏模型资源unitypackage
- Jurassic Pack Vol. III Dinosaurs 侏罗纪包卷恐龙三号Unity游戏模型资源unitypackag
- Ultimate Seating Controller 终极座椅控制器Unity游戏开发插件资源unitypackage
- 什么是人工智能-关于人工智能的相关介绍说明
- Figma Converter for Unity适用Unity的Figma转换器Unity游戏开发插件unitypackage
- Creepy Animatronic Anims 令人毛骨悚然的电子动画Unity游戏动画插件资源unitypackage
资源上传下载、课程学习等过程中有任何疑问或建议,欢迎提出宝贵意见哦~我们会及时处理!
点击此处反馈
评论6