EMD经验模式算法去杂波、滤波处理程序，已调试MATLAB.zip

%% EMD分解&重构 20150710-ByZXQ function Signal=Funplot_hht00(signal,Ts,imfk) % Plot the HHT. % plot_hht(x,Ts) % % :: Syntax % The array x is the input signal and Ts is the sampling period. % Example on use: [x,Fs] = wavread('Hum.wav'); % plot_hht(x(1:6000),1/Fs); % Func : emd % Get HHT.imfk为重构信号中需要挖除的imf分量 %% 数据分割 x=real(signal); y=imag(signal); % x=real(x); %%%%%%%%%%%%%%%%% real EMD %%%%%%%%%%%%%%%%%%%% % s1 = Fungetspline(x);%上包络线 % s2 = -Fungetspline(-x);%上包络线 % x1 = (s1+s2)/2;%此处的x2为文章中的hhel % figure; % % plot(N,x);xlabel('sample'), ylabel('Amplitude');title('双边带调幅信号');hold on; % % plot(N,s1,'-r'); % % plot(N,s2,'-r'); % % plot(N,x1,'g'); % plot(t,x);xlabel('Times'), ylabel('Amplitude');title('直达波时域波形');hold on; % plot(t,s1,'-r'); % plot(t,s2,'-r'); % plot(t,x1,'g'); imf = FunEMD(x); % for k = 1:length(imf) % b(k) = sum(imf{k}.*imf{k}); % th = angle(hilbert(imf{k})); % d{k} = diff(th)/Ts/(2*pi); % end % [u,v] = sort(-b); % b = 1-b/max(b); % % Set time-frequency plots. % N = length(x); % c = linspace(0,(N-2)*Ts,N-1); % figure; % for k = v(1:2) % plot(c,d{k},'k.','Color',b([k k k]),'MarkerSize',3); hold on; % set(gca,'FontSize',8,'XLim',[0 c(end)],'YLim',[0 50]);%设置x、y轴句柄 % xlabel('Time'), ylabel('Frequency');title('原信号时频图'); % end % % Set IMF plots. M = length(imf); N = length(x); c = linspace(0,(N-1)*Ts,N); % for k1 = 0:4:M-1 % figure % for k2 = 1:min(4,M-k1), subplot(4,1,k2), plot(c,imf{k1+k2}); set(gca,'FontSize',8,'XLim',[0 c(end)]);title('实部EMD分解结果');end % xlabel('Time'); % end %%%%%%%%%%%%%%%%% imag EMD %%%%%%%%%%%%%%%%%%%% % s1 = Fungetspline(y);%上包络线 % s2 = -Fungetspline(-y);%上包络线 % y1 = (s1+s2)/2;%此处的x2为文章中的hhel imf = FunEMD(y); % for k = 1:length(imf) % b(k) = sum(imf{k}.*imf{k}); % th = angle(hilbert(imf{k})); % d{k} = diff(th)/Ts/(2*pi); % end % [u,v] = sort(-b); % b = 1-b/max(b); % % Set time-frequency plots. % N = length(y); % c = linspace(0,(N-2)*Ts,N-1); % figure; % for k = v(1:2) % plot(c,d{k},'k.','Color',b([k k k]),'MarkerSize',3); hold on; % set(gca,'FontSize',8,'XLim',[0 c(end)],'YLim',[0 50]);%设置x、y轴句柄 % xlabel('Time'), ylabel('Frequency');title('原信号时频图'); % end % % Set IMF plots. M = length(imf); N = length(y); c = linspace(0,(N-1)*Ts,N); % for k1 = 0:4:M-1 % figure % for k2 = 1:min(4,M-k1), subplot(4,1,k2), plot(c,imf{k1+k2}); set(gca,'FontSize',8,'XLim',[0 c(end)]);title('虚部EMD分解结果');end % xlabel('Time'); % end % %%%%%%%%%%%%%%%%%%% 确定干扰位置 %%%%%%%%%%%%%%%%%% % var=std(imf{1});%取第一个imf分量确定干扰位置 % win=1:10;%滑窗大小 % par=imf{1}(win); % num=1e4/50; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%总长为1e4 % for i=1:1:num-1 % var(i)=std(par); % win=win+50; % % if win<=9990 % par=imf{1}(win); % % end % end % for i=1:1:num-1 % [x,index]=find( var(i)<var); % end % INdex=[index]; %满足条件的位置 % %% 抑制干扰 % m=length(INdex); % for i=1:1:m % % win=1:100;%滑窗大小 % Win=INdex(i)*50+1:(INdex(i)+1)*50; % imf{1}(Win)=0; % imf{2}(Win)=0;%置为0 % % end % figure % plot(c,imf{1}); %% 重构信号 X=0; %若imf{1}为挖除重构，则X=imf{1}；直接去除X=0 for k1 = imfk:M-1 %k1=1为imf{1}，k=0则需要挖除重构 X=X+imf{k1}; end Y=0; %若imf{1}为挖除重构，则X=imf{1} for k1 = imfk:M-1 %起始k1=1为imf{1}，k=0则需要挖除重构 Y=Y+imf{k1}; end % figure; plot(c,X); % figure; subplot(2,1,1); % plot(real(y_BF_1)); % subplot(2,1,2); % % plot(c,X); % figure; plot(c,Y); % figure; subplot(2,1,1); % plot(imag(y_BF_1)); % subplot(2,1,2); % plot(c,Y); %% 重构信号 Signal=X+1i*Y; %% 模糊函数结果 % Funambiguity(y_BF,y_Scan) % save([Out_p,'y_Scan_20.mat'],'y_Scan_20');