logo资料库

MatlabEMDHHTSVM-程序及图形.doc

第1页 / 共7页
第2页 / 共7页
第3页 / 共7页
第4页 / 共7页
第5页 / 共7页
第6页 / 共7页
第7页 / 共7页
资料共7页,全文预览结束
原始信号:(脑电图) EMD 分解图:
EMD 程序: % EMD_Gao: Emprical Mode Decomposition by Gao Qiang % % To get the EMD (imf) of signal X % % imf = emd_Gao(x_sig) % x_sig % % imf : % % % See: Huang, Royal Society Proceedings on Math, Physical, % % % % arrange in Matrix, in row, with residual in last row and Engineering Sciences, vol. 454, no. 1971, pp. 903-995, 8 March 1998 : input sinal (row vector) intrinsic mode function
% function imf = emd(x_sig) clc close all; % x=1:1000; % y=sin(x)+cos(x); % x_sig=y; % c=x_sig;% copy signal to c % load F:\杨件数据\jizhong.mat; % x=jizhong([1:5000]); % x=jizhong; x=x_sig; c=x'; N = length(x); % length of x % SDstop = 0.3; SDstop = 0.25; % standard deviation value when sift process stop imf = []; % Matrix which will contain the successive IMF, and the residue % loop to calculate the imf of signal x while (1) h = c; SD = 1; % criterion to stop sift process, % h is the data series to be processinged standard deviation while (SD > SDstop) maxVec = []; minVec = []; % look for max and min point for i = 2: N - 1 if h (i - 1) < h (i) & h (i) > h (i + 1) maxVec = [maxVec i]; % get max position end if h (i - 1) > h (i) & h (i) < h (i + 1) minVec = [minVec i]; % get min position end end % check if it is residual if (length (maxVec) + length (minVec)) < 2 break; end
% % % % % % % % % % % % handle end point if h (1) >h (2) maxVec = [1 maxVec]; else minVec = [1 minVec]; end if h (N) > h (N - 1) maxVec = [maxVec N]; else minVec = [minVec N]; end % left end point % right end point %handle end point lenmax=length(maxVec); lenmin=length(minVec); %left end point if h(1)>0 if(maxVec(1)0 if(maxVec(lenmax)
else if(maxVec(lenmax)
figure(i) plot(imf(i,:)); % for i=1:line % % % end subplot(line+1,1,1),plot(x); for i=1:line subplot(line+1,1,i+1),plot(imf(i,:)); end % figure(1) % plot(imf); % figure(2) % plot(imf'); return; HHT 程序: % emd_HHT : % % % imfX % imfHSp % % % to compute the Hilbert spectrum of every imf and then get instantaneous frequency of them : input data, imf : output data, the Hilbert spectrum of all imf function imfHSp = emd_HHT (imfX) aX = imfX; % get imf iNo = size (imfX, 1) - 1; % imf number dL = size (imfX, 2); imfHSp = []; % clear imf for i = 1 : iNo % get one imf imfA = aX (i, :); y = hilbert (imfA); % hilbert transform fai = angle(y); % get phase angle fai = unwrap(fai); omig = diff (fai); omig = [omig omig = omig / (2 * pi); % get angle speed omig(dL - 1) ];% add a data to omig at last position after diff
imfHSp = [imfHSp; omig]; end if iNo == 0 imfA = aX ; % get one imf y = hilbert (imfA); % hilbert transform fai = angle(y); % get phase angle fai = unwrap(fai); omig = diff (fai); omig = [omig omig = omig / (2 * pi); imfHSp = [imfHSp; omig]; % get angle speed omig(dL - 1) ];% add a data to omig at last position after diff end return;
分享到:
收藏