EMD与EEMD程序.doc

%载入信号x=load(1.txt);%产生信号N=length(x); %采样点数fs=2000; %采样频率dt=1/fs; %采样时间间隔t=(0:N-1)*dt; %产生时间序列%EMDimf=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. Gonalvs% 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 threshold) tolerance% &% |#zeros-#extrema| 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 % 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) & nbitMAXITERATIONS/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) 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(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) 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 closeendEEMD程序% 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章节和内容）% % 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; endend 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 0X(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 0X(max)ans = 1 5 11 153）、函数极小值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 Y轴标签zlabel Z轴标签contour 等高线图contourc 等高线计算contourf 填充的等高线图hidden 网格线消影meshc 连接网格/等高线mesh 具有参考轴的3D网格peaks 具有两个变量的采样函数surf 3D阴影表面图su*ce 建立表面低层对象surfc 海浪和等高线的结合surfl 具有光照的3D阴影表面trimesh 三角网格图1 常用指令(General Purpose Commands)1.1 通用信息查询(General information)demo 演示程序help 在线帮助指令helpbrowser 超文本文档帮助信息helpdesk 超文本文档帮助信息helpwin 打开在线帮助窗info MATLAB 和MathWorks 公司的信息subscribe MATLAB 用户注册ver MATLAB 和TOOLBOX 的版本信息version MATLAB 版本whatsnew 显示版本新特征1.2 工作空间管理(Managing the workspace)clear 从内存中清除变量和函数exit 关闭MATLABload 从磁盘中调入数据变量pack 合并工作内存中的碎块quit 退出MATLABsave 把内存变量存入磁盘who 列出工作内存中的变量名whos 列出工作内存中的变量细节workspace 工作内存浏览器1.3 管理指令和函数(Managing commands and functions)edit 矩阵编辑器edit 打开M 文件inmem 查看内存中的P 码文件mex 创建MEX 文件open 打开文件pcode 生成P 码文件type 显示文件内容what 列出当前目录上的M、MAT、MEX 文件which 确定指定函数和文件的位置1.4 搜索路径的管理(Managing the seach patli)addpath 添加搜索路径rmpath 从搜索路径中删除目录path 控制MATLAB 的搜索路径pathtool 修改搜索路径1.5 指令窗控制(Controlling the command window)beep 产生beep 声echo 显示命令文件指令的切换开关diary 储存MATLAB 指令窗操作内容format 设置数据输出格式more 命令窗口分页输出的控制开关1.6 操作系统指令(Operating system commands)cd 改变当前工作目录computer 计算机类型copyfile 文件拷贝delete 删除文件dir 列出的文件dos 执行dos 指令并返还结果getenv 给出环境值ispc MATLAB 为PC(Windows)版本则为真isunix MATLAB 为Unix 版本则为真mkdir 创建目录pwd 改变当前工作目录unix 执行unix 指令并返还结果vms 执行vms dcl 指令并返还结果web 打开web 浏览器! 执行外部应用程序2 运算符和特殊算符(Operators and special characters)2.1 算术运算符(Arithmetic operators)+ 加- 减* 矩阵乘.* 数组乘 矩阵乘方. 数组乘方 反斜杠或左除/ 斜杠或右除./或. 数组除张量积注本表第三栏括号中的字符供在线救助时help 指令引述用2.2 关系运算符(Relational operators)= = 等号= 不等号 大于 = 大于或等于2.3 逻辑操作(Logical operators)& 逻辑与| 逻辑或 逻辑非xor 异或any 有非零元则为真all 所有元素均非零则为真2.4 特殊算符(Special characters)： 冒号( ) 圆括号 方括号 花括号 创建函数句柄. 小数点. 构架域的关节点. 父目录 续行号, 逗号; 分号% 注释号! 调用操作系统命令= 赋值符号 引号 复数转置号. 转置号, 水平串接; 垂直串接( ), ,. 下标赋值( ), ,. 下标标识subsindex 下标标识3 编程语言结构(Programming language constructs)3.1 控制语句(Control flow)break 终止最内循环case 同switch 一起使用catch 同try 一起使用continue 将控制转交给外层的for 或while 循环else 同if 一起使用elseif 同if 一起使用end 结束for,while,if 语句for 按规定次数重复执行语句if 条件执行语句otherwise 可同switch 一起使用return 返回switch 多个条件分支try try-cathch 结构while 不确定次数重复执行语句3.2 计算运行(Evaluation and execution)assignin 跨空间赋值builtin 执行内建的函数eval 字符串宏指令evalc 执行MATLAB 字符串evalin 跨空间计算串表达式的值feval 函数宏指令run 执行脚本文件3.3 脚本文件、函数及变量(Scripts,function,and variables)exist 检查变量或函数是否被定义function 函数文件头global 定义全局变量isglobal 若是全局变量则为真iskeyword 若是关键字则为真mfilename 正在执行的M 文件的名字persistent 定义永久变量script MATLAB 命令文件3.4 宗量处理(Augument handling)inputname 实际调用变量名nargchk 输入变量个数检查nargin 函数输入宗量的个数nargout 函数输出宗量的个数nargoutchk 输出变量个数检查varagin 输入宗量varagout 输出宗量3.5 信息显示(Message display)disp 显示矩阵和文字内容display 显示矩阵和文字内容的重载函数error 显示错误信息fprintf 把格式化数据写到文件或屏幕lasterr 最后一个错误信息lastwarn 最后一个警告信息sprintf 按格式把数字转换为串warning 显示警告信息3.6 交互式输入(Interactive input)input 提示键盘输入keyboard 激活键盘做为命令文件pause 暂停uicontrol 创建用户界面控制uimenu 创建用户界面菜单4 基本矩阵函数和操作(Elementary matrices and matrix manipulation)4.1 基本矩阵(Elementary matrices)eye 单位阵linspace 线性等分向量logspace 对数等分向量meshgrid 用于三维曲面的分格线坐标ones 全1 矩阵rand 均匀分布随机阵randn 正态分布随机阵repmat 铺放模块数组zeros 全零矩阵