• shirley03
    了解作者
  • matlab
    开发工具
  • 2.6KB
    文件大小
  • rar
    文件格式
  • 0
    收藏次数
  • 5 积分
    下载积分
  • 1
    下载次数
  • 2022-01-04 16:56
    上传日期
变分模态分解程序,直接调用,适用于多分量信号分解
VMD.rar
  • VMD
  • lizi.m
    852B
  • VMD.m
    3.4KB
  • fenjielizi.m
    695B
内容介绍
function [u, u_hat, omega] = VMD(signal, alpha, tau, K, DC, init, tol) % VMD % Input and Parameter % signal - 要分解的信号 % alpha - 数据保真约束的平衡参数 % tau - 双上升的时间步长(选择0作为噪声碎片) % K - 要得到的模式数目 % DC - 如果第一种模式被放置并保持在DC(0-freq),则为true % init - 0= 所有的omega从0开始 % 1 = 所有的omega开始均匀分布 % 2 =所有的omega都是随机初始化的 % tol - 收敛准则的限度;通常在1e-6左右 % Output: % u - 分解模态的集合 % u_hat -模态的频谱 % omega - 估计的模态中心频率 %输入信号的周期和采样频率 save_T = length(signal); fs = 1/save_T; % 镜像延拓 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; % 镜像信号的时域0-T T = length(f); t = (1:T)/T; % 谱域离散化 freqs = t-0.5-1/T; % 最大迭代次数 N = 1000; % 进一步推广:每个模式都有单独的alpha 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; % 矩阵跟踪每次迭代// u_hat_plus = zeros(N, length(freqs), K); %初始化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 % 如果有DC模态,设它的omega为0 if DC omega_plus(1,1) = 0; end % 从空对偶变量开始 lambda_hat = zeros(N, length(freqs)); % other inits uDiff = tol+eps; % update step n = 1; % loop counter sum_uk = 0; % accumulator % 迭代更新的主循环 while ( uDiff > tol && n < N ) % 没有收敛且低于迭代限制 %更新第一模态的累加数 k = 1; sum_uk = u_hat_plus(n,:,K) + sum_uk - u_hat_plus(n,:,1); % 通过残差维纳滤波器更新第一模态频谱 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 % 更新其他模态 for k=2:K % 累加数 sum_uk = u_hat_plus(n+1,:,k-1) + sum_uk - u_hat_plus(n,:,k); % 模态频谱 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); % 中心频率 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 % 双重提升 lambda_hat(n+1,:) = lambda_hat(n,:) + tau*(sum(u_hat_plus(n+1,:,:),3) - f_hat_plus); % 循环计数器 n = n+1; % 收敛否? 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 % 后处理和清理 % 如果提前收敛,则丢弃多余部分 N = min(N,n); omega = omega_plus(1:N,:); % 信号恢复 u_hat = zeros(T, K); u_hat((T/2+1):T,:) = squeeze(u_hat_plus(N,(T/2+1):T,:)); %squeeze:除去size为1的维度 u_hat((T/2+1):-1:2,:) = squeeze(conj(u_hat_plus(N,(T/2+1):T,:))); u_hat(1,:) = conj(u_hat(end,:)); %conj:求复共轭 u = zeros(K,length(t)); for k = 1:K u(k,:)=real(ifft(ifftshift(u_hat(:,k)))); end % 除去镜像部分 u = u(:,T/4+1:3*T/4); % 重新计算频谱 clear u_hat; for k = 1:K u_hat(:,k)=fftshift(fft(u(k,:)))'; end end
评论
    相关推荐