EMD与EEMD程序.doc

%%%%%%%%%%%载入信号
x=load('1.txt');%产生信号
N=length(x); %采样点数
fs=2000; %采样频率
dt=1/fs; %采样时间间隔
t=(0:N-1)*dt; %产生时间序列
%%%%%%%%EMD
imf=emd(x);
% EMD.M (EMD程序)
% G. Rilling, July 2002
%
% computes EMD (Empirical Mode Decomposition) according to:
%
% N. E. Huang et al., "The empirical mode decomposition and the
% Hilbert spectrum for non-linear and non stationary time series analysis,"
% Proc. Royal Soc. London A, Vol. 454, pp. 903-995, 1998
%
% with variations reported in:
%
% G. Rilling, P. Flandrin and P. Gonçalvès
% "On Empirical Mode Decomposition and its algorithms"
% IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing
% NSIP-03, Grado (I), June 2003
%
% stopping criterion for sifting :
% at each point : mean amplitude < threshold2*envelope amplitude
% &
% mean of boolean array ((mean amplitude)/(envelope amplitude) > threshold) < tolerance
% &
% |#zeros-#extrema|<=1
%
% inputs: - x : analysed signal (line vector)
% - t (optional) : sampling times (line vector) (default : 1:length(x))
% - stop (optional) : threshold, threshold2 and tolerance (optional)
% for sifting stopping criterion
% default : [0.05,0.5,0.05]
% - tst (optional) : if equals to 1 shows sifting steps with pause
% if equals to 2 no pause
%
% outputs: - imf : intrinsic mode functions (last line = residual)
% - ort : index of orthogonality
% - nbits : number of iterations for each mode
%
% calls: - extr (finds extrema and zero-crossings)
% - io : computes the index of orthogonality

function [imf,ort,nbits] = emd(x,t,stop,tst);

% default for stopping
defstop = [0.05,0.5,0.05];

if(nargin==1)
t = 1:length(x);
stop = defstop;
tst = 0;
end

if(nargin==2)
stop = defstop;
tst = 0;
end

if (nargin==3)
tst=0;
end

S = size(x);
if ((S(1) > 1) & (S(2) > 1)) | (length(S) > 2)
error('x must have only one row or one column')
end

if S(1) > 1
x = x';
end

S = size(t);
if ((S(1) > 1) & (S(2) > 1)) | (length(S) > 2)
error('t must have only one row or one column')
end

if S(1) > 1
t = t';
end

if (length(t)~=length(x))
error('x and t must have the same length')
end

S = size(stop);
if ((S(1) > 1) & (S(2) > 1)) | (S(1) > 3) | (S(2) > 3) | (length(S) > 2)
error('stop must have only one row or one column of max three elements')
end

if S(1) > 1
stop = stop';
S = size(stop);
end

if S(2) < 3
stop(3)=defstop(3);
end

if S(2) < 2
stop(2)=defstop(2);
end

sd = stop(1);
sd2 = stop(2);
tol = stop(3);

if tst
figure
end

% maximum number of iterations
MAXITERATIONS=2000;

% maximum number of symmetrized points for interpolations
NBSYM = 2;

lx = length(x);

sdt(lx) = 0;
sdt = sdt+sd;
sd2t(lx) = 0;
sd2t = sd2t+sd2;

% maximum number of extrema and zero-crossings in residual
ner = lx;
nzr = lx;

r = x;
imf = [];
k = 1;

% iterations counter for extraction of 1 mode
nbit=0;

% total iterations counter
NbIt=0;

while ner > 2

% current mode
m = r;

% mode at previous iteration
mp = m;

sx = sd+1;

% tests if enough extrema to proceed
test = 0;

[indmin,indmax,indzer] = extr(m);
lm=length(indmin);
lM=length(indmax);
nem=lm + lM;
nzm=length(indzer);

j=1;

% sifting loop
while ( mean(sx > sd) > tol | any(sx > sd2) | (abs(nzm-nem)>1)) & (test == 0) & nbit<MAXITERATIONS

if(nbit>MAXITERATIONS/5 & mod(nbit,floor(MAXITERATIONS/10))==0)
disp(['mode ',int2str(k),' nombre d iterations : ',int2str(nbit)])
disp(['stop parameter mean value : ',num2str(s)])
end

% boundary conditions for interpolations :

if indmax(1) < indmin(1)
if m(1) > m(indmin(1))
lmax = fliplr(indmax(2:min(end,NBSYM+1)));
lmin = fliplr(indmin(1:min(end,NBSYM)));
lsym = indmax(1);
else
lmax = fliplr(indmax(1:min(end,NBSYM)));
lmin = [fliplr(indmin(1:min(end,NBSYM-1))),1];
lsym = 1;
end
else

if m(1) < m(indmax(1))
lmax = fliplr(indmax(1:min(end,NBSYM)));
lmin = fliplr(indmin(2:min(end,NBSYM+1)));
lsym = indmin(1);
else
lmax = [fliplr(indmax(1:min(end,NBSYM-1))),1];
lmin = fliplr(indmin(1:min(end,NBSYM)));
lsym = 1;
end
end

if indmax(end) < indmin(end)
if m(end) < m(indmax(end))
rmax = fliplr(indmax(max(end-NBSYM+1,1):end));
rmin = fliplr(indmin(max(end-NBSYM,1):end-1));
rsym = indmin(end);
else
rmax = [lx,fliplr(indmax(max(end-NBSYM+2,1):end))];
rmin = fliplr(indmin(max(end-NBSYM+1,1):end));
rsym = lx;
end
else
if m(end) > m(indmin(end))
rmax = fliplr(indmax(max(end-NBSYM,1):end-1));
rmin = fliplr(indmin(max(end-NBSYM+1,1):end));
rsym = indmax(end);
else
rmax = fliplr(indmax(max(end-NBSYM+1,1):end));
rmin = [lx,fliplr(indmin(max(end-NBSYM+2,1):end))];
rsym = lx;
end
end

tlmin = 2*t(lsym)-t(lmin);
tlmax = 2*t(lsym)-t(lmax);
trmin = 2*t(rsym)-t(rmin);
trmax = 2*t(rsym)-t(rmax);

% in case symmetrized parts do not extend enough
if tlmin(1) > t(1) | tlmax(1) > t(1)
if lsym == indmax(1)
lmax = fliplr(indmax(1:min(end,NBSYM)));
else
lmin = fliplr(indmin(1:min(end,NBSYM)));
end
if lsym == 1
error('bug')
end
lsym = 1;
tlmin = 2*t(lsym)-t(lmin);
tlmax = 2*t(lsym)-t(lmax);
end

if trmin(end) < t(lx) | trmax(end) < t(lx)
if rsym == indmax(end)
rmax = fliplr(indmax(max(end-NBSYM+1,1):end));
else
rmin = fliplr(indmin(max(end-NBSYM+1,1):end));
end
if rsym == lx
error('bug')
end
rsym = lx;
trmin = 2*t(rsym)-t(rmin);
trmax = 2*t(rsym)-t(rmax);
end

mlmax =m(lmax);
mlmin =m(lmin);
mrmax =m(rmax);
mrmin =m(rmin);

% definition of envelopes from interpolation

envmax = interp1([tlmax t(indmax) trmax],[mlmax m(indmax) mrmax],t,'spline');
envmin = interp1([tlmin t(indmin) trmin],[mlmin m(indmin) mrmin],t,'spline');

envmoy = (envmax + envmin)/2;

m = m - envmoy;

[indmin,indmax,indzer] = extr(m);
lm=length(indmin);
lM=length(indmax);
nem = lm + lM;
nzm = length(indzer);

% evaluation of mean zero
sx=2*(abs(envmoy))./(abs(envmax-envmin));
s = mean(sx);

% display

if tst
subplot(4,1,1)
plot(t,mp);hold on;
plot(t,envmax,'--k');plot(t,envmin,'--k');plot(t,envmoy,'r');

title(['IMF ',int2str(k),'; iteration ',int2str(nbit),' before sifting']);
set(gca,'XTick',[])
hold off

subplot(4,1,2)
plot(t,sx)
hold on
plot(t,sdt,'--r')
plot(t,sd2t,':k')
title('stop parameter')
set(gca,'XTick',[])
hold off

subplot(4,1,3)
plot(t,m)
title(['IMF ',int2str(k),'; iteration ',int2str(nbit),' after sifting']);
set(gca,'XTick',[])

subplot(4,1,4);
plot(t,r-m)
title('residue');
disp(['stop parameter mean value : ',num2str(s)])
if tst == 2
pause(0.01)
else
pause
end

end

% end loop : stops if not enough extrema
if nem < 3
test = 1;
end

mp = m;
nbit=nbit+1;
NbIt=NbIt+1;

if(nbit==(MAXITERATIONS-1))
warning(['forced stop of sifting : too many iterations... mode ',int2str(k),'. stop parameter mean value : ',num2str(s)])
end

end
imf(k,:) = m;
if tst
disp(['mode ',int2str(k),' enregistre'])
end
nbits(k) = nbit;
k = k+1;
r = r - m;
[indmin,indmax,indzer] = extr(r);
ner = length(indmin) + length(indmax);
nzr = length(indzer);
nbit=1;

if (max(r) - min(r)) < (1e-10)*(max(x) - min(x))
if ner > 2
warning('forced stop of EMD : too small amplitude')
else

disp('forced stop of EMD : too small amplitude')
end
break
end

end

imf(k,:) = r;

ort = io(x,imf);

if tst
close
end

EEMD程序
% This is an EMD/EEMD program
%
% function allmode=eemd(Y,Nstd,NE)
%
% INPUT:
% Y: Inputted data;（输入数据）
% Nstd: ratio of the standard deviation of the added noise and that of Y;（添加噪声和Y的标准偏差的比率）
% NE: Ensemble number for the EEMD（EEMD集合的数量）
% OUTPUT:
% A matrix of N*(m+1) matrix, where N is the length of the input
% data Y, and m=fix(log2(N))-1. Column 1 is the original data, columns 2, 3, ...（一个N *(m + 1)矩阵,其中N是输入数据Y的长度，和m =fix(log2(N))-1。第1列是原始数据,列2,3,…）
% m are the IMFs from high to low frequency, and comlumn (m+1) is
% the residual (over all trend). （m是imf的由高到低的频率,和comlumn(m + 1是余量)）
%
% NOTE:
% It should be noted that when Nstd is set to zero and NE is set to
% 1, the program degenerates to a EMD program.（应该指出的是,当Nstd设置为零和NE设置为1，程序退化到一个EMD程序。）
%
% References can be found in the "Reference" section.（参考可以发现在“参考”一节。）
%
% The code is prepared by Zhaohua. For questions, please read the "Q&A"
% section or contact（该代码是Zhaohua Wu准备的，如果有问题，请读"Q&A"章节和内容）
% [e-mail]
%

function allmode=eemd(Y,Nstd,NE)
xsize=length(Y);
dd=1:1:xsize;
Ystd=std(Y);
Y=Y/Ystd;

TNM=fix(log2(xsize))-1;
TNM2=TNM+2;
for kk=1:1:TNM2,
for ii=1:1:xsize,
allmode(ii,kk)=0.0;
end
end

for iii=1:1:NE,
for i=1:xsize,
temp=randn(1,1)*Nstd;(产生随机数)
X1(i)=Y(i)+temp;
end

for jj=1:1:xsize,
mode(jj,1) = Y(jj);
end

xorigin = X1;
xend = xorigin;

nmode = 1;（节点数=1）
while nmode <= TNM,
xstart = xend;
iter = 1;（迭代次数=1）

while iter<=10,
[spmax, spmin, flag]=extrema(xstart);
upper= spline(spmax(:,1),spmax(:,2),dd);
lower= spline(spmin(:,1),spmin(:,2),dd);
mean_ul = (upper + lower)/2;（求最大最小包络的平均）
xstart = xstart - mean_ul;
iter = iter +1;
end
xend = xend - xstart;

nmode=nmode+1; （节点数+1）

for jj=1:1:xsize,
mode(jj,nmode) = xstart(jj);
end
end

for jj=1:1:xsize,
mode(jj,nmode+1)=xend(jj);
end

allmode=allmode+mode;

end

allmode=allmode/NE;
allmode=allmode*Ystd;

findpeaks 求极值
diff(f)求f的导数
峰值点可以用imregionalmax, 零点用find就可以了
峰值通过find（x==max(x))就可以解决
max只能找到最高的峰值，不是楼主所需，举个例子
Y=[10 0 10 20 30 20 10 0 30 50 70 50 30 0 50 0];
X=1:size(Y,2)；
max=imregionalmax(Y)

max =

1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0
X(max)

ans =

1 5 11 15
极值点可以用imregionalmax和imregionalmin，零点可以用find，这里如果用find（x==max(x))这条命令只能找到一个值
举个例子
Y=[10 0 10 20 30 20 10 0 30 50 70 50 30 0 50 0];
X=1:size(Y,2)；
max=imregionalmax(Y)

max =

1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0
X(max)

ans =

1 5 11 15
3）、函数极小值
A、一元函数的极小值
fminbnd：求得函数在给定区间内的局部极小值。
该函数的调用格式为
[x ,fval]= fminbnd(fun,x1,x2)
输入变量：fun 为函数；x1 和 x2 分别用于指定区间的左右边界。
输出变量：x为极小值位置，fval—极小值。
例：求humps函数在[0.4 0.8]区间上的极小值
[x ,fval]= fminbnd('humps',0.4,0.8)
思考：如何求humps函数在[ 0.8，1]区间上的极大值？
%函数
function T=Nearest(a)
syms t;
y=(a-t).^2;
Y=sum(y);
f=diff(Y,t);
T=solve(f);

%验证
a=[1,2,3,4,5,6,78,45,7,8,9,0];
t=Nearest(a)

%结果
t=14

解释一下，可以先列个函数，y=sum(a(i)-t)^2),要求的就是y取极值（即最值）时对应的t值，那么对y求导数，求出导函数的零点，极为所求的t值。

%结果

[x,fval] = fminsearch(函数, 自变量初始值)
注意：函数表达式的自变量为一向量。

matlab常用函数与常用指令大全
matlab常用函数- -
1、特殊变量与常数
ans 计算结果的变量名
computer 确定运行的计算机
eps 浮点相对精度
Inf 无穷大
I 虚数单位
inputname 输入参数名
NaN 非数
nargin 输入参数个数
nargout 输出参数的数目
pi 圆周率
nargoutchk 有效的输出参数数目
realmax 最大正浮点数
realmin 最小正浮点数
varargin 实际输入 的参量
varargout 实际返回的参量
操作符与特殊字符
+ 加 - 减
* 矩阵乘法 .* 数组乘（对应元素相乘）
^ 矩阵幂 .^ 数组幂（各个元素求幂）
\ 左除或反斜杠 / 右除或斜面杠
./ 数组除（对应元素除）

kron Kronecker张量积
: 冒号 () 圆括
[] 方括 . 小数点
.. 父目录 ... 继续
, 逗号（分割多条命令） ; 分号（禁止结果显示）
% 注释 ! 感叹号
' 转置或引用 = 赋值
== 相等 <> 不等于
& 逻辑与 | 逻辑或
~ 逻辑非 xor 逻辑异或

2、基本数学函数
abs 绝对值和复数模长
acos,acodh 反余弦，反双曲余弦
acot,acoth 反余切，反双曲余切
acsc,acsch 反余割，反双曲余割
angle 相角
asec,asech 反正割，反双曲正割
secant 正切
asin,asinh 反正弦，反双曲正弦
atan,atanh 反正切，双曲正切
tangent 正切
atan2 四象限反正切
ceil 向着无穷大舍入
complex 建立一个复数
conj 复数配对
cos,cosh 余弦，双曲余弦
csc,csch 余切，双曲余切
cot,coth 余切，双曲余切
exp 指数
fix 朝0方向取整
floor 朝负无穷取整
*** 最大公因数
imag 复数值的虚部
lcm 最小公倍数
log 自然对数
log2 以2为底的对数
log10 常用对数
mod 有符号的求余

nchoosek 二项式系数和全部组合数
real 复数的实部
rem 相除后求余
round 取整为最近的整数
sec,sech 正割，双曲正割
sign 符号数
sin,sinh 正弦，双曲正弦
sqrt 平方根
tan,tanh 正切，双曲正切
3、基本矩阵和矩阵操作
blkding 从输入参量建立块对角矩阵
eye 单位矩阵
linespace 产生线性间隔的向量
logspace 产生对数间隔的向量
numel 元素个数
ones 产生全为1的数组
rand 均匀颁随机数和数组
randn 正态分布随机数和数组
zeros 建立一个全0矩阵 colon) 等间隔向量
cat 连接数组
diag 对角矩阵和矩阵对角线
fliplr 从左自右翻转矩阵
flipud 从上到下翻转矩阵
repmat 复制一个数组

reshape 改造矩阵
roy90 矩阵翻转90度

tril 矩阵的下三角
triu 矩阵的上三角

dot 向量点集
cross 向量叉集

ismember 检测一个集合的元素
intersect 向量的交集

setxor 向量异或集
setdiff 向是的差集

union 向量的并集
数值分析和傅立叶变换
cumprod 累积

cumsum 累加
cumtrapz 累计梯形法计算数值微分

factor 质因子
inpolygon 删除多边形区域内的点

max 最大值
mean 数组的均值

mediam 中值
min 最小值

perms 所有可能的转换
polyarea 多边形区域

primes 生成质数列表
prod 数组元素的乘积

rectint 矩形交集区域
sort 按升序排列矩阵元素

sortrows 按升序排列行
std 标准偏差

sum 求和
trapz 梯形数值积分

var 方差
del2 离散拉普拉斯

diff 差值和微分估计
gradient 数值梯度

cov 协方差矩阵
corrcoef 相关系数

conv2 二维卷积
conv 卷积和多项式乘法

filter IIR或FIR滤波器
deconv 反卷积和多项式除法

filter2 二维数字滤波器
cplxpair 将复数值分类为共轭对

fft 一维的快速傅立叶变换
fft2 二维快速傅立叶变换

fftshift 将FFT的DC分量移到频谱中心
ifft 一维快速反傅立叶变换

ifft2 二维傅立叶反变换
ifftn 多维快速傅立叶变换

ifftshift 反FFT偏移
nextpow2 最靠近的2的幂次

unwrap 校正相位角

多项式与插值
conv 卷积和多项式乘法

roots 多项式的根
poly 具有设定根的多项式

polyder 多项式微分
polyeig 多项式的特征根

polyfit 多项式拟合
polyint 解析多项式积分

polyval 多项式求值
polyvalm 矩阵变量多项式求值

residue 部分分式展开
interp1 一维插值

interp2 二维插值
interp3 三维插值

interpft 使用FFT的一维插值
interpn 多维插值

meshgrid 为3维点生成x和y的网格
ndgrid 生成多维函数和插值的数组

pchip 分段3次Hermite插值多项式
ppval 分段多项式的值

spline 3次样条数据插值

绘图函数
bar 竖直条图

barh 水平条图
hist 直方图

histc 直方图计数
hold 保持当前图形

loglog x,y对数坐标图
pie 饼状图

plot 绘二维图
polar 极坐标图

semilogy y轴对数坐标图
semilogx x轴对数坐标

subplot 绘制子图
bar3 数值3D竖条图

bar3h 水平3D条形图
comet3 3D慧星图

cylinder 圆柱体
fill3 填充的3D多边形

plot3 3维空间绘图
quiver3 3D震动（速度）图

slice 体积薄片图
sphere 球

stem3 绘制离散表面数据
wate***ll 绘制瀑布

trisurf 三角表面
clabel 增加轮廓标签到等高线图中

datetick 数据格式标记
grid 加网格线

gtext 用鼠标将文本放在2D图中
legend 图注

plotyy 左右边都绘Y轴
title 标题

xlabel X轴标签
ylabel