VMD算法的MATLAB实现

%主程序%--------------- Preparationclear all;close all;clc;% Time Domain 0 to TT = 1000;fs = 1/T;t = (1:T)/T;freqs = 2*pi*(t-0.5-1/T)/(fs);% center frequencies of componentsf_1 = 2;f_2 = 24;f_3 = 288;% modesv_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 purposeswsub{1} = 2*pi*f_1;wsub{2} = 2*pi*f_2;wsub{3} = 2*pi*f_3;% composite signal, including noisef = v_1 + v_2 + v_3 + 0.1*randn(size(v_1));% some sample parameters for VMDalpha = 2000;        % moderate bandwidth constrainttau = 0;            % noise-tolerance (no strict fidelity enforcement)K = 4;              % 4 modesDC = 0;             % no DC part imposedinit = 1;           % initialize omegas uniformlytol = 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)+1subplot(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 codeimf = 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)+1subplot(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 signalsave_T = length(signal);fs = 1/save_T;% extend the signal by mirroringT = 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 discretizationfreqs = 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 modeAlpha = alpha*ones(1,K);% Construct and center f_hatf_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 memu_hat_plus = zeros(N, length(freqs), K);% Initialization of omega_komega_plus = zeros(N, K);switch initcase 1for i = 1:Komega_plus(1,i) = (0.5/K)*(i-1);endcase 2omega_plus(1,:) = sort(exp(log(fs) + (log(0.5)-log(fs))*rand(1,K)));otherwiseomega_plus(1,:) = 0;end% if DC mode imposed, set its omega to 0if DComega_plus(1,1) = 0;end% start with empty dual variableslambda_hat = zeros(N, length(freqs));% other initsuDiff = tol+eps; % update stepn = 1; % loop countersum_uk = 0; % accumulator% ----------- Main loop for iterative updateswhile ( uDiff > tol &&  n < N ) % not converged and below iterations limit% update first mode accumulatork = 1;sum_uk = u_hat_plus(n,:,K) + sum_uk - u_hat_plus(n,:,1);% update spectrum of first mode through Wiener filter of residualsu_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 0if ~DComega_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 modefor k=2:K% accumulatorsum_uk = u_hat_plus(n+1,:,k-1) + sum_uk - u_hat_plus(n,:,k);% mode spectrumu_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 frequenciesomega_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 ascentlambda_hat(n+1,:) = lambda_hat(n,:) + tau*(sum(u_hat_plus(n+1,:,:),3) - f_hat_plus);% loop countern = n+1;% converged yet?uDiff = eps;for i=1:KuDiff = 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)))';enduDiff = abs(uDiff);end%------ Postprocessing and cleanup% discard empty space if converged earlyN = min(N,n);omega = omega_plus(1:N,:);% Signal reconstructionu_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:Ku(k,:)=real(ifft(ifftshift(u_hat(:,k))));end% remove mirror partu = u(:,T/4+1:3*T/4);% recompute spectrumclear u_hat;for k = 1:Ku_hat(:,k)=fftshift(fft(u(k,:)))';endend

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

VMD算法的MATLAB实现相关推荐

VMD分解，matlab代码，包络线，包络谱，中心频率，峭度值，能量熵，近似熵，包络熵，频谱图，希尔伯特变换，包含所有程序MATLAB代码，-西储大学数据集为例
目录 1.选取数据 2.VMD函数-matlab代码 3.采用matlab脚本导入数据并做VMD分解 4. VMD分解图 5.计算中心频率 6.画包络线 7. 画包络谱 8. 计算峭度值 9.计算能量 ...
VMD分解及其matlab实现方法
VMD是一种新型的信号分解方法,它基于Hilbert Huang变换(HHT)理论,可以将一个信号分解成多个正交的模态,每个模态都有自己的中心频率和频率带宽.VMD的优点在于,能够克服传统的信号分解方 ...
fcm算法的MATLAB实现,FCM算法的matlab程序（初步）
FCM算法的matlab程序 1.采用iris数据库 iris_data.txt 5.1 3.5 1.4 0.2 4.9 3 1.4 0.2 4.7 3.2 1.3 0.2 4.6 3.1 1.5 0 ...
2018-4-8蚁群算法---包子阳《智能优化算法以及Matlab实现》第五章
资料来源: <智能优化算法以及matlab实现>包子阳余继周编著第五章-----蚁群算法是一种元启发式优化算法(自己理解:就是作为群体的单位个体也就是元,在里面充当着随机的选择搜 ...
matlab dfp法,DFP算法及Matlab程序.docx
DFP算法及Matlab程序作业二用DFP算法求解,取,.一.求解:求迭代点x1令,得的极小值点,所以得:于是,由DFP修正公式有下一个搜索方向为求迭代点x2令,得的极小值点于是得:,所以:,因H ...
matlab整定串级pid,PID算法在Matlab串级控制中的应用
PID算法在Matlab串级控制中的应用自114 1112002039 陈艳前言:这个专题是由王娟老师给我们授课,我感觉收获挺大的,尤其是matlab仿真软件的使用,为我以后的实验课打下良好的基础 ...
PSO-LSSVM算法及其MATLAB代码
挺完整的一篇博客,这里转载记录一下. 原文链接:PSO-LSSVM算法及其MATLAB代码一.PSO 1.概念粒子群优化算法(PSO:Particle swarm optimization)是一种 ...
matlab hist函数_算法工匠MATLAB专训营：Matlab绘图，小试牛刀
作者 | 蔡老师仿真秀专栏作者首发 | 仿真秀平台导读:正文之前,我在此详细说明一下,因为本文包含的程序太难得,网上肯定找不到这样的程序.随着讲课的越来越深入,我给出的程序会越来越实用,接近于实 ...
matlab音频基频的提取,(620512681) 自相关基频提取算法的MATLAB实现
第31卷总第80期西北民族大学学报(自然科学版) V01．31．No．4 1 0年1 2 0 2月 Journal of Nonhw铭t University for Nationalities(N ...

VMD算法的MATLAB实现

VMD算法的MATLAB实现相关推荐

最新文章

热门文章