数据转换基于时变滤波器的经验模态分解（TVF-EMD）附matlab代码.rar

% Time varying filtering based EMD function imf = tvf_emd(x) %input % x: input signal x(t) %output % imf: resulting intrinsic mode functions %%%%%%%%%%% preprocessing THRESH_BWR=0.1; % inst bandwidth threshold (stopping criterion) MAX_IMF=50; % max output imfs BSP_ORDER=26; % b-spline order end_flag=0; imf=zeros(MAX_IMF,numel(x)); if size(x,1) > 1 % convert to row vector x = x.'; end temp_x=x; t=1:numel(x); for nimf=1:MAX_IMF [indmin_x, indmax_x] = extr(temp_x); % is there enough extrema to continue if nimf==MAX_IMF imf(nimf,:)=temp_x; nimf=nimf+1; break; end if(numel([indmin_x, indmax_x])<4) imf(nimf,:)=temp_x; if numel(find(temp_x~=0))>0 nimf=nimf+1; end end_flag=1; end if end_flag==1 break; end num_padding = round(length(temp_x)*0.5);% padding number y = temp_x; flag_stopiter=0; for iter=1:100 y = [fliplr(y(2:2+num_padding-1)) y fliplr(y(end-num_padding:end-1))]; % padding to deal with boundary effect (symmetric) %y=extendsignal(y,num_padding); % padding to deal with boundary effect (match slope) tt=1:numel(y); ind_remov_pad=num_padding+1:numel(y)-num_padding; %%%%%%%%%%%%% extract inst amp and inst freq using hilbert transform [indmin_y, indmax_y] = extr(y); indexC_y=sort([indmin_y, indmax_y]); [instAmp0,instFreq0] = INST_FREQ_local(y); %%%%%%%%%%%%%% divide y(t) into two parts and compute bis freq (cutoff freq) [a1 f1 a2 f2 bis_freq instBWR avgFreq]=devide_y(y,instAmp0,instFreq0); instBWR2=instBWR; for j=1:2:numel(indexC_y)-2 ind=indexC_y(j):indexC_y(j+2); instBWR2(ind)=mean(instBWR(ind)); end bis_freq(instBWR2<THRESH_BWR)=1e-12; bis_freq(bis_freq>0.5)=0.45; bis_freq(bis_freq<=0)=1e-12; %%%%%%%%%%%% deal with mode mixing bis_freq = anti_modemixing(y,bis_freq,ind_remov_pad,num_padding); bis_freq=bis_freq(ind_remov_pad); bis_freq = [fliplr(bis_freq(2:2+num_padding-1)) bis_freq fliplr(bis_freq(end-num_padding:end-1))]; bis_freq = anti_modemixing(y,bis_freq,ind_remov_pad,num_padding); bis_freq=bis_freq(ind_remov_pad); bis_freq = [fliplr(bis_freq(2:2+num_padding-1)) bis_freq fliplr(bis_freq(end-num_padding:end-1))]; %%%%%%%%%%%%%% stopping criteria temp_instBWR=instBWR2(ind_remov_pad); ind_start=round(numel(temp_instBWR)*0.05);% incase boundary effect ind_end=round(numel(temp_instBWR)*0.95); if (iter>=2 && mean(temp_instBWR(ind_start:ind_end))<THRESH_BWR+THRESH_BWR/4*iter) || iter>=6 || (nimf>1 && mean(temp_instBWR(ind_start:ind_end))<THRESH_BWR+THRESH_BWR/4*iter) flag_stopiter=1; end if numel(find(temp_instBWR(ind_start:ind_end)>THRESH_BWR))/numel(instBWR2(ind_remov_pad))<0.2 % most of the signal are local narrowband flag_stopiter=1; end %%%%%%%%%%%%%% obtain local mean using time varying filtering phi=zeros(1,numel(bis_freq)); % build knots via h(t) for i=1:numel(bis_freq)-1 phi(i+1)=phi(i)+2*pi*bis_freq(i); end [indmin_knot, indmax_knot] = extr(cos(phi)); indexC_knot=sort([indmin_knot, indmax_knot]); if numel(indexC_knot)>2 pp_spline = splinefit(1:length(y),y,indexC_knot,BSP_ORDER,'p'); % TVF filtering with the extrema of h(t) as knots localmean = ppval(pp_spline,1:length(y)); else flag_stopiter=1; %not enough knots to perform filtering end if (max(abs(y(ind_remov_pad)-localmean(ind_remov_pad)))/min(abs(localmean(ind_remov_pad)))<1e-3) % prevent signal of very small amplitude to be extracted flag_stopiter=1; end temp_residual=y-localmean; temp_residual=temp_residual(ind_remov_pad); temp_residual=temp_residual(round(numel(temp_residual)*0.1):end-round(numel(temp_residual)*0.1)); localmean2=localmean(ind_remov_pad); localmean2=localmean2(round(numel(localmean2)*0.1):end-round(numel(localmean2)*0.1)); if abs(max(localmean2))/abs(max(instAmp0(ind_remov_pad)))<3.5e-2 || abs(max(temp_residual))/abs(max(instAmp0(ind_remov_pad)))<1e-2 flag_stopiter=1; end %%%%%%%% check stopping criteria if flag_stopiter imf(nimf,:)=y(ind_remov_pad); temp_x=temp_x-y(ind_remov_pad); break; end y=y-localmean; y=y(ind_remov_pad); end end imf(nimf:MAX_IMF,:)=[]; end function output_cutoff = anti_modemixing(y,bis_freq,ind_remov_pad,num_padding) org_bis_freq=bis_freq; output_cutoff=bis_freq; flag_intermitt=0; t=1:numel(bis_freq); intermitt=[]; %%%%% locate the maxima of the input signal [indmin_y, indmax_y] = extr(y); indexC_y=sort([indmin_y, indmax_y]); zero_span=[]; %%%% preprocessing for i=2:numel(indmax_y)-1 time_span=indmax_y(i-1):indmax_y(i+1); if (max(bis_freq(time_span))-min(bis_freq(time_span)))/min(bis_freq(time_span))>0.25 zero_span=[zero_span time_span]; end end bis_freq(zero_span)=0; %%%% find out all intermittences diff_bis_freq=zeros(size(bis_freq)); for i=1:numel(indmax_y)-1 time_span=indmax_y(i):indmax_y(i+1); if (max(bis_freq(time_span))-min(bis_freq(time_span)))/min(bis_freq(time_span))>0.25 intermitt=[intermitt indmax_y(i)]; diff_bis_freq(indmax_y(i))=bis_freq(indmax_y(i+1))-bis_freq(indmax_y(i)); end end ind_remov_pad([1:round(0.1*numel(ind_remov_pad)),round(0.9*numel(ind_remov_pad)):end])=[]; inters=intersect(ind_remov_pad,intermitt); if numel(inters) >0 flag_intermitt=1; end %%%% find out floors % plot(t(intermitt),bis_freq(intermitt),'r.',t,bis_freq,'b') for i=2:numel(intermitt)-1 u1=intermitt(i-1); u2=intermitt(i); u3=intermitt(i+1); % check the derivative of cutoff frequency if diff_bis_freq(u2)>0 % rising edge bis_freq(u1:u2)=0; end if diff_bis_freq(u2)<0 % falling edge bis_freq(u2:u3)=0; end end temp_bis_freq=bis_freq; temp_bis_freq(temp_bis_freq<1e-9)=0; %cutoff freq of very small value is considered to be zero (floor value) temp_bis_freq=temp_bis_freq(ind_remov_pad); temp_bis_freq = [fliplr(temp_bis_freq(2:2+num_padding-1)) temp_bis_freq fliplr(temp_bis_freq(end-num_padding:end-1))]; flip_bis_freq=fliplr(bis_freq); if numel(find(temp_bis_freq>1e-9))>0 && numel(find(flip_bis_freq>1e-9,1,'first'))>0 temp_bis_freq(1)=bis_freq(find(bis_freq>1e-9,1,'first'));temp_bis_freq(end)=flip_bis_freq(find(flip_bis_freq>1e-9,1,'first')); %padding boundary for interpolation else temp_bis_freq(1)=bis_freq(1); temp_bis_freq(end)=bis_freq(end); end bis_freq=temp_bis_freq; %%%% interpolate between peaks if numel(t(bis_freq~=0))<2 return; end bis_freq = interp1(t(bis_freq~=0),bis_freq(bis_freq~=0) ,t,'pchip'); flip_bis_freq=fliplr(org_bis_freq); if numel(find(org_bis_freq>1e-9))>0 && numel(find(flip_bis_freq>1e-9,1,'first'))>0 org_bis_freq(1)=org_bis_freq(find(org_b