马老四
里论外几-
%此版本为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 = x1;
x = x-x1;
end
imf = 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, u = 1; end
function s = getspline(x)
%三次样条函数拟合成元数据包络线
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 = find(diff(diff(x) > 0) < 0);
u = find(x(n+1) > x(n));
n(u) = n(u)+1;
------------------------------------------------------------------------------------------
----------------------------------------------------------------------------------------
function plot_hht00(x,Ts)
% 双边带调幅信号的EMD分解
% 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.
clear all;
close all;
Ts=0.0005;
t=0:Ts:10; % 采样率2000HZ
% 调幅信号
%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.*imf);
th = angle(hilbert(imf));
d = 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.","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);
set(gca,"FontSize",8,"XLim",[0 c(end)]);
title("EMD分解结果");
end
xlabel("Time");
end