EMD算法的matlab程序介绍.docx

发布时间：2022-05-31 发布人：admin 分类：说明书资料大小：0.01M 资料格式：docx 举报版权申诉

qq_40757244-10209074-16359647595213999634.docx.pdf-第1页.png

第1页 / 共3页

qq_40757244-10209074-16359647595213999634.docx.pdf-第2页.png

第2页 / 共3页

qq_40757244-10209074-16359647595213999634.docx.pdf-第3页.png

第3页 / 共3页

文本预览

%此版本为 ALAN 版本的整合注释版 function imf = emd(x) % Empiricial Mode Decomposition (Hilbert-Huang Transform) % imf = emd(x) % Func : findpeaks x= transpose(x(:));%转置为行矩阵 imf = []; while ~ismonotonic(x) %当 x 不是单调函数，分解终止条件 x1 = x; sd = Inf;%均值 %直到 x1 满足 IMF 条件，得 c1 while (sd > 0.1) | ~isimf(x1) %当标准偏差系数 sd 大于 0.1 或 x1 不是固有模态函数时，分量终止条件 s1 = getspline(x1);%上包络线 s2 = -getspline(-x1);%下包络线 x2 = x1-(s1+s2)/2;%此处的 x2 为文章中的 h sd = sum((x1-x2).^2)/sum(x1.^2); x1 = x2; end imf{end+1} = x1; x end = x-x1; imf{end+1} = x; % FUNCTIONS function u = ismonotonic(x) %u=0 表示 x 不是单调函数，u=1 表示 x 为单调的 u1 = length(findpeaks(x))*length(findpeaks(-x)); if u1 > 0, u = 0; else, u = 1; end function u = isimf(x) %u=0 表示 x 不是固有模式函数，u=1 表示 x 是固有模式函数 N = length(x); u1 = sum(x(1:N-1).*x(2:N) < 0); u2 = length(findpeaks(x))+length(findpeaks(-x)); if abs(u1-u2) > 1, u = 0; else, function s = getspline(x) u = 1; end

%三次样条函数拟合成元数据包络线 N = length(x); p = findpeaks(x); s = spline([0 p N+1],[0 x(p) 0],1:N); ------------------------------------------------------------------------------- -------------------------------------------------------------------------------- function n = findpeaks(x) % Find peaks.找到极值 ,n 为极值点所在位置 % n = findpeaks(x) n u n(u) = n(u)+1; = find(diff(diff(x) > 0) < 0); = find(x(n+1) > x(n)); ------------------------------------------------------------------------------------------ ---------------------------------------------------------------------------------------- function plot_hht00(x,Ts) % 双边带调幅信号的 EMD 分解 % Plot the HHT. % plot_hht(x,Ts) % % :: Syntax % % % % Func : emd % Get HHT. clear all; close all; Ts=0.0005; t=0:Ts:10; % 采样率 2000HZ 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); % 调幅信号 %x=sin(2*pi*t).*sin(40*pi*t); x=sin(2*pi*t); s1 = getspline(x);%上包络线 s2 = -getspline(-x);%上包络线 x1 = (s1+s2)/2;%此处的 x2 为文章中的 h figure; plot(t,x);xlabel('Time'), ylabel('Amplitude');title('双边带调幅信号');hold on; plot(t,s1,'-r'); plot(t,s2,'-r'); plot(t,x1,'g');

imf = emd(x); for k = 1:length(imf) b(k) = sum(imf{k}.*imf{k}); th d{k} = diff(th)/Ts/(2*pi); = angle(hilbert(imf{k})); = 1-b/max(b); end [u,v] = sort(-b); 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

分享到：

赞收藏

资料库

EMD算法的matlab程序介绍.docx

相关推荐

人工智能

热门标签

最新资料