• zhjhappy
    了解作者
  • matlab
    开发工具
  • 2KB
    文件大小
  • zip
    文件格式
  • 0
    收藏次数
  • 1 积分
    下载积分
  • 6
    下载次数
  • 2021-03-28 22:33
    上传日期
盲源分离方面的算法,使用的是VMD算法对单通道信号进行盲分离。
w.zip
  • w
  • VMD2.m
    7KB
内容介绍
%主程序 %--------------- Preparation clear all; close all; clc; % Time Domain 0 to T T = 1000; fs = 1/T; t = (1:T)/T; freqs = 2*pi*(t-0.5-1/T)/(fs); % center frequencies of components f_1 = 2; f_2 = 24; f_3 = 288; % modes v_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 purposes wsub{1} = 2*pi*f_1; wsub{2} = 2*pi*f_2; wsub{3} = 2*pi*f_3; % composite signal, including noise f = v_1 + v_2 + v_3 + 0.1*randn(size(v_1)); % some sample parameters for VMD alpha = 2000; % moderate bandwidth constraint tau = 0; % noise-tolerance (no strict fidelity enforcement) K = 4; % 4 modes DC = 0; % no DC part imposed init = 1; % initialize omegas uniformly tol = 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)+1 subplot(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 code imf = 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)+1 subplot(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 signal save_T = length(signal); fs = 1/save_T; % 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; % 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
评论
    相关推荐
    • 视频会议系统.rar
      一套商用的视频会议程序(cool)。
    • 视频捕获-mpeg4压缩-网络传送源码(VC 6.0).rar
      视频捕获-mpeg4压缩-网络传送源码(VC 6.0).rar
    • VOIP_SRC_ALL.rar
      搞VOIP的一些原代码,包括MGCP,SIP,RTP/TRCP,SRTP(加密RTP)...等源码,参考价值无限!!这可是我所有的收藏阿,我都献给大家了:=)
    • mpeg.zip
      MPEG-2和MP3编解码
    • live.2008.01.04.tar.rar
      网络直播源代码,支持RTSP、RTP,支持的格式包括MPEG-1/2/4,MP3,H264等,功能强大,流媒体开发必须。
    • handvu-beta3.tar.gz
      结合摄像头使用,用于手部位置跟踪,精准度非常高,另可以识别六种手势。其手位跟踪模块可作为手势识别软件的一个组成部分,内含非常详细的说明文档。编译该代码,需要 OpenCV 0.97 或更高版本的支持。在Windows下的编译顺序为cubicles工程、handvu工程,之后根据视频数据来源要求编译其他代码。建议使用 Visual Studio 2003 ,如果需要编译 DirectX 相关模块,则请安装 DirectX SDK 9.0 或更高版本。
    • encore50src.zip
      mpeg4/divx编码程序,多谢qing_lee
    • 语音识别.rar
      老师布置的一个project,里面有题目和源程序,能够识别十个数字。
    • rtsp-1.0.1.tar.gz
      RTSP的实现代码,并附有RTSP、RTP、RTCP等协议的文档!
    • mpge4.rar
      一组关于流媒体、Mpeg4的技术文档,可供做流媒体的朋友参考。