EMD样例.txt
- 格式:doc
- 大小:48.00 KB
- 文档页数:7
(一)简单的EMD程序function imf = emd(x)% Empiricial Mode Decomposition (Hilbert-Huang Transform)% imf = emd(x)% Func : findpeaksx = transpose(x(:));%转置imf = [];while ~ismonotonic(x) %当x不是单调函数,分解终止条件x1 = x;sd = Inf;%均值%直到x1满足IMF条件,得c1while (sd > 0.1) | ~isimf(x1) %当标准偏差系数sd大于0.1或x1不是固有模态函数时,分量终止条件s1 = getspline(x1);%上包络线s2 = -getspline(-x1);%下包络线x2 = x1-(s1+s2)/2;%此处的x2为文章中的hsd = sum((x1-x2).^2)/sum(x1.^2);x1 = x2;endimf{end+1} = x1;x = x-x1;endimf{end+1} = x;% FUNCTIONSfunction u = ismonotonic(x)%u=0表示x不是单调函数,u=1表示x为单调的u1 = length(findpeaks(x))*length(findpeaks(-x));if u1 > 0, u = 0;else, u = 1; endfunction 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; endfunction 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 = 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:1; % 采样率2000HZ% 调幅信号x=sin(2*pi*t).*sin(40*pi*t);s1 = getspline(x);%上包络线s2 = -getspline(-x);%上包络线x1 = (s1+s2)/2;%此处的x2为文章中的hfigure;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 = 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-1figurefor 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分解结果');endxlabel('Time');end(二)function imf = emd(x)% Empiricial Mode Decomposition (Hilbert-Huang Transform)% imf = emd(x)% Func : findpeaksx = transpose(x(:));%转置imf = [];while ~ismonotonic(x) %当x不是单调函数,分解终止条件x1 = x;sd = Inf;%均值%直到x1满足IMF条件,得c1while (sd > 0.1) | ~isimf(x1) %当标准偏差系数sd大于0.1或x1不是固有模态函数时,分量终止条件s1 = getspline(x1);%上包络线s2 = -getspline(-x1);%下包络线x2 = x1-(s1+s2)/2;%此处的x2为文章中的hsd = sum((x1-x2).^2)/sum(x1.^2);x1 = x2;endimf{end+1} = x1;x = x-x1;endimf{end+1} = x;% FUNCTIONSfunction u = ismonotonic(x)%u=0表示x不是单调函数,u=1表示x为单调的u1 = length(findpeaks(x))*length(findpeaks(-x));if u1 > 0, u = 0;else, u = 1; endfunction 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; endfunction 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 = 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:1; % 采样率2000HZ% 调幅信号x=sin(2*pi*t).*sin(40*pi*t);s1 = getspline(x);%上包络线s2 = -getspline(-x);%上包络线x1 = (s1+s2)/2;%此处的x2为文章中的hfigure;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 = 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-1figurefor 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分解结果');endxlabel('Time');endplot_hht00(x,Ts)只画出了时频图,没有三维图。
EMDEmpirical Mode Decomposition经验模态分解美国工程院院士黄锷1998年提出一种自适应数据处理或挖掘方法,适用于非线性、非平稳时间序列的处理。
1.什么是平稳和非平稳时间序列的平稳,一般是宽平稳,即时间序列的方差和均值是和时间无关的常数,协方差与与时间间隔有关、与时间无关。
未来样本时间序列,其均值、方差、协方差必定与已经获得的样本相同,理解为平稳的时间序列是有规律且可预测的,样本拟合曲线的形态具有“惯性”。
而非平稳信号样本的本质特征只存在于信号所发生的当下,不会延续到未来,不可预测。
严格来说实际上不存在理想平稳序列,实际情况下都是非平稳。
2.什么是EMD经验模态分解方法?EMD理论上可以应用于任何类型时间序列信号的分解,在实际工况中大量非平稳信号数据的处理上具有明显优势。
这种优势是相对于建立在先验性假设的谐波基函数上的傅里叶分解和小波基函数上的小波分解而言的。
EMD分解信号不需要事先预定或强制给定基函数,而是依赖信号本身特征自适应地进行分解。
相对于小波分解:EMD克服了基函数无自适应性的问题,小波分析需要选定一个已经定义好的小波基,小波基的选择至关重要,一旦选定,在整个分析过程中无法更换。
这就导致全局最优的小波基在局部的表现可能并不好,缺乏适应性。
而EMD不需要做预先的分析与研究,可以直接开始分解,不需要人为的设置和干预。
相对于傅里叶变换:EMD克服了传统傅里叶变换中用无意义的谐波分量来表示非线性、非平稳信号的缺点,并且可以得到极高的时频分辨率。
EMD方法的关键是将复杂信号分解为有限个本征模函数IMF,Intrinsic Mode Function。
分解出来的IMF分量包含了原信号的不同时间尺度上的局部特征信号。
这句话中:不同时间尺度=局部平稳化,通过数据的特征时间尺度来获得本征波动模式,然后分解or筛选数据。
本质上,EMD将一个频率不规则的波化为多个单一频率的波+残波的形式。
EMD(经验模态分解)算法三EMD(经验模态分解)算法三经验模态分解(EMD)算法是一种用于信号和数据分解的信号处理方法,用于提取信号中的本征模态函数(IMFs)。
其主要思想是将信号分解为一系列本征模态函数,每个本征模态函数代表一个具有特定频率和幅值的本征振动模式。
该算法已被广泛应用于信号处理、图像处理、数据分析等领域。
EMD算法的基本步骤如下:1.将待分解的信号表示为一个局部极值点的峰谷序列。
2.通过连接相邻局部极值点,构建一系列包络线。
3.将原始信号与包络线之差作为细节信号,重复步骤1和步骤2,直到细节信号达到其中一种停止条件。
4.将分解出的所有细节信号相加得到分解后的信号。
具体来说,EMD算法的主要步骤如下:1.初始化。
将原始信号记为x(t),并设置初始模态函数集合为空。
令h(t)=x(t)。
2.局部极值点提取。
在h(t)中寻找所有局部极大值点和局部极小值点,记为m(t)和n(t)。
3.插值。
通过对局部极大值点和局部极小值点之间的过零点进行三次样条插值,得到包络线e(t)。
4.分离。
将原始信号x(t)减去包络线e(t),得到细节信号d(t)。
令h(t)=d(t)。
5.判断停止条件。
判断细节信号d(t)是否满足其中一种停止条件,如果满足则停止分解,否则返回步骤26.更新模态函数集合。
将e(t)添加到模态函数集合中。
7.分解。
将细节信号d(t)作为新的原始信号,重复步骤2至步骤6EMD算法的优点是不依赖于模型假设,能够适应多种类型的信号和数据。
它能够在时域和频域上对信号进行分解,提取信号中的局部特征,具有较好的局部适应性和高精度。
然而,EMD算法也存在一些问题。
首先,EMD算法对噪声非常敏感,在存在较高噪声的情况下,容易产生过分分解和模态混叠的问题。
其次,EMD算法的计算复杂度较高,随着信号长度的增加,计算时间也会增加。
为了解决EMD算法存在的问题,研究者提出了许多改进算法,如快速EMD算法(FEMD)、改进的EMD算法(CEEMD)等。
主题:EMD算法Python实现一、介绍EMD算法的概念1. EMD算法全称为Earth Mover's Distance,中文意为“地球移动距离”,是一种用来衡量两个分布之间的相似性的算法。
2. EMD算法最早由Y. Rubner等人在1998年提出,是一种基于距离度量的非线性的分布对齐方法。
3. EMD算法被广泛应用于图像处理、信号处理、文本分析等领域,具有很好的实际应用价值。
4. EMD算法的核心思想是通过将一个分布转换为另一个分布的最小代价来计算两个分布之间的距离。
二、EMD算法Python实现的基本原理1. EMD算法的实现需要解决一个最小化问题,即寻找两个分布之间的最小代价。
2. 在Python中,可以使用scipy包中的optimize模块来实现EMD算法,该模块提供了优化算法的实现,可以直接调用进行分布对齐计算。
3. EMD算法的实现可以分为以下几个步骤:1)将两个分布表示为柱状图;2)计算每个柱状图之间的距离矩阵;3)使用optimize模块中的线性规划函数来求解最小代价;4)根据求解结果得到两个分布之间的EMD距离。
三、EMD算法Python实现的具体步骤1. 导入必要的Python库import numpy as npfrom scipy.optimize import linear_sum_assignment2. 定义两个分布distribution1 = np.array([0.3, 0.5, 0.2])distribution2 = np.array([0.4, 0.4, 0.2])3. 计算距离矩阵distance_matrix = np.zeros((len(distribution1),len(distribution2)))for i in range(len(distribution1)):for j in range(len(distribution2)):distance_matrix[i][j] = abs(i - j)4. 使用线性规划函数求解最小代价row_ind, col_ind = linear_sum_assignment(distance_matrix)5. 根据求解结果计算EMD距离emd_distance = sum([distance_matrix[row_ind[i]][col_ind[i]] * distribution1[row_ind[i]] for i in range(len(row_ind))])四、实例演示假设有两个分布分别为distribution1 = [0.3, 0.5, 0.2]和distribution2 = [0.4, 0.4, 0.2],我们可以利用上述Python实现的EMD算法来计算它们之间的距禿:distribution1 = np.array([0.3, 0.5, 0.2])distribution2 = np.array([0.4, 0.4, 0.2])distance_matrix = np.zeros((len(distribution1),len(distribution2)))for i in range(len(distribution1)):for j in range(len(distribution2)):distance_matrix[i][j] = abs(i - j)row_ind, col_ind = linear_sum_assignment(distance_matrix) emd_distance = sum([distance_matrix[row_ind[i]][col_ind[i]] * distribution1[row_ind[i]] for i in range(len(row_ind))])print("The EMD distance between distribution1 and distribution2 is:", emd_distance)五、总结1. EMD算法是一种用于计算两个分布之间距离的算法,具有广泛的应用价值。
%%%%%%%%%%%载入信号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,July2002%%computesEMD(EmpiricalModeDecomposition)accordingto:%%N.E.Huangetal.,"Theempiricalmodedecompositionandthe%Hilbertspectrumfornon-linearandnonstationarytimeseriesanalysis,"%Proc.RoyalSoc.LondonA,Vol.454,pp.903-995,1998%%withvariationsreportedin:%%G.Rilling,P.FlandrinandP.Gon?alvès %"OnEmpiricalModeDecompositionanditsalgorithms"%IEEE-EURASIPWorkshoponNonlinearSignalandImageProcessing%NSIP-03,Grado(I),June2003%%stoppingcriterionforsifting:%ateachpoint:meanamplitude<threshold2*envelopeamplitude%&%meanofbooleanarray((meanamplitude)/(envelopeamplitude)>threshold)<tolerance %&%|#zeros-#extrema|<=1%%inputs:-x:analysedsignal(linevector)%-t(optional):samplingtimes(linevector)(default:1:length(x))%-stop(optional):threshold,threshold2andtolerance(optional)%forsiftingstoppingcriterion%default:[0.05,0.5,0.05]%-tst(optional):ifequalsto1showssiftingstepswithpause%ifequalsto2nopause%%outputs:-imf:intrinsicmodefunctions(lastline=residual)%-ort:indexoforthogonality%-nbits:numberofiterationsforeachmode%%calls:-extr(findsextremaandzero-crossings)%-io:computestheindexoforthogonalityfunction[imf,ort,nbits]=emd(x,t,stop,tst);%defaultforstoppingdefstop=[0.05,0.5,0.05];if(nargin==1)t=1:length(x);stop=defstop;tst=0;endif(nargin==2)stop=defstop;tst=0;endif(nargin==3)tst=0;endS=size(x);if((S(1)>1)&(S(2)>1))|(length(S)>2)error('xmusthaveonlyoneroworonecolumn')endif S(1)>1x=x';endS=size(t);if((S(1)>1)&(S(2)>1))|(length(S)>2)error('tmusthaveonlyoneroworonecolumn')endif S(1)>1t=t';endif(length(t)~=length(x))error('xandtmusthavethesamelength')endS=size(stop);if((S(1)>1)&(S(2)>1))|(S(1)>3)|(S(2)>3)|(length(S)>2)error('stopmusthaveonlyoneroworonecolumnofmaxthreeelements') endif S(1)>1stop=stop';S=size(stop);endif S(2)<3stop(3)=defstop(3);endif S(2)<2stop(2)=defstop(2);endsd=stop(1);sd2=stop(2);tol=stop(3);if tstfigureend%maximumnumberofiterationsMAXITERATIONS=2000;%maximumnumberofsymmetrizedpointsforinterpolationsNBSYM=2;lx=length(x);sdt(lx)=0;sdt=sdt+sd;sd2t(lx)=0;sd2t=sd2t+sd2;%maximumnumberofextremaandzero-crossingsinresidualner=lx;nzr=lx;r=x;imf=[];k=1;%iterationscounterforextractionof1modenbit=0;%totaliterationscounterNbIt=0;while ner>2%currentmodem=r;%modeatpreviousiterationmp=m;sx=sd+1;%testsifenoughextrematoproceedtest=0;[indmin,indmax,indzer]=extr(m);lm=length(indmin);lM=length(indmax);nem=lm+lM;nzm=length(indzer);j=1;%siftingloopwhile(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),'nombrediterations:',int2str(nbit)])disp(['stopparametermeanvalue:',num2str(s)])end%boundaryconditionsforinterpolations: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);elselmax=fliplr(indmax(1:min(end,NBSYM)));lmin=[fliplr(indmin(1:min(end,NBSYM-1))),1];lsym=1;endelseif m(1)<m(indmax(1))lmax=fliplr(indmax(1:min(end,NBSYM)));lmin=fliplr(indmin(2:min(end,NBSYM+1)));lsym=indmin(1);elselmax=[fliplr(indmax(1:min(end,NBSYM-1))),1];lmin=fliplr(indmin(1:min(end,NBSYM)));lsym=1;endendif 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);elsermax=[lx,fliplr(indmax(max(end-NBSYM+2,1):end))]; rmin=fliplr(indmin(max(end-NBSYM+1,1):end));rsym=lx;endelseif 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);elsermax=fliplr(indmax(max(end-NBSYM+1,1):end));rmin=[lx,fliplr(indmin(max(end-NBSYM+2,1):end))]; rsym=lx;endendtlmin=2*t(lsym)-t(lmin);tlmax=2*t(lsym)-t(lmax);trmin=2*t(rsym)-t(rmin);trmax=2*t(rsym)-t(rmax);%incasesymmetrizedpartsdonotextendenoughif tlmin(1)>t(1)|tlmax(1)>t(1)if lsym==indmax(1)lmax=fliplr(indmax(1:min(end,NBSYM)));elselmin=fliplr(indmin(1:min(end,NBSYM)));endif lsym==1error('bug')endlsym=1;tlmin=2*t(lsym)-t(lmin);tlmax=2*t(lsym)-t(lmax);endif trmin(end)<t(lx)|trmax(end)<t(lx)if rsym==indmax(end)rmax=fliplr(indmax(max(end-NBSYM+1,1):end));elsermin=fliplr(indmin(max(end-NBSYM+1,1):end));endif rsym==lxerror('bug')endrsym=lx;trmin=2*t(rsym)-t(rmin);trmax=2*t(rsym)-t(rmax);endmlmax=m(lmax);mlmin=m(lmin);mrmax=m(rmax);mrmin=m(rmin);%definitionofenvelopesfrominterpolationenvmax=interp1([tlmaxt(indmax)trmax],[mlmaxm(indmax)mrmax],t,'spline'); envmin=interp1([tlmint(indmin)trmin],[mlminm(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);%evaluationofmeanzerosx=2*(abs(envmoy))./(abs(envmax-envmin));s=mean(sx);%displayif tstsubplot(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),'beforesifting']);set(gca,'XTick',[])hold offsubplot(4,1,2)plot(t,sx)hold onplot(t,sdt,'--r')plot(t,sd2t,':k')title('stopparameter')set(gca,'XTick',[])hold offsubplot(4,1,3)plot(t,m)title(['IMF',int2str(k),';iteration',int2str(nbit),'aftersifting']);set(gca,'XTick',[])subplot(4,1,4);plot(t,r-m)title('residue');disp(['stopparametermeanvalue:',num2str(s)])if tst==2pause(0.01)elsepauseendend%endloop:stopsifnotenoughextremaif nem<3test=1;endmp=m;nbit=nbit+1;NbIt=NbIt+1;if(nbit==(MAXITERATIONS-1))warning(['forcedstopofsifting:toomanyiterations...mode',int2str(k),'.stopparametermea nvalue:',num2str(s)])endendimf(k,:)=m;if tstdisp(['mode',int2str(k),'enregistre'])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>2warning('forcedstopofEMD:toosmallamplitude')elsedisp('forcedstopofEMD:toosmallamplitude')endbreakendendimf(k,:)=r;ort=io(x,imf);if tstcloseendEEMD程序%ThisisanEMD/EEMDprogram%%functionallmode=eemd(Y,Nstd,NE)%%INPUT:%Y:Inputteddata;(输入数据)%Nstd:ratioofthestandarddeviationoftheaddednoiseandthatofY;(添加噪声和Y的标准偏差的比率)%NE:EnsemblenumberfortheEEMD(EEMD集合的数量)%OUTPUT:%AmatrixofN*(m+1)matrix,whereNisthelengthoftheinput%dataY,andm=fix(log2(N))-1.Column1istheoriginaldata,columns2,3,...(一个N*(m+1)矩阵,其中N是输入数据Y的长度,和m=fix(log2(N))-1。
winver---------检查Windows版本wmimgmt.msc----打开windows管理体系结构(WMI)write----------写字板mspaint--------画图板mstsc----------远程桌面连接dfrg.msc-------磁盘碎片整理程序notepad--------打开记事本ddeshare-------打开DDE共享设置cmd命令全集winver---------检查Windows版本wmimgmt.msc----打开windows管理体系结构(WMI) wupdmgr--------windows更新程序wscript--------windows脚本宿主设置write----------写字板winmsd---------系统信息wiaacmgr-------扫描仪和照相机向导winchat--------XP自带局域网聊天mem.exe--------显示内存使用情况Msconfig.exe---系统配置实用程序mplayer2-------简易widnows media player mspaint--------画图板mstsc----------远程桌面连接mplayer2-------媒体播放机magnify--------放大镜实用程序mmc------------打开控制台mobsync--------同步命令dxdiag---------检查DirectX信息drwtsn32------ 系统医生devmgmt.msc--- 设备管理器dfrg.msc-------磁盘碎片整理程序diskmgmt.msc---磁盘管理实用程序dcomcnfg-------打开系统组件服务ddeshare-------打开DDE共享设置dvdplay--------DVD播放器net stop messenger-----停止信使服务net start messenger----开始信使服务notepad--------打开记事本nslookup-------网络管理的工具向导ntbackup-------系统备份和还原narrator-------屏幕“讲述人”ntmsmgr.msc----移动存储管理器ntmsoprq.msc---移动存储管理员操作请求netstat -an----(TC)命令检查接口syncapp--------创建一个公文包sysedit--------系统配置编辑器sigverif-------文件签名验证程序sndrec32-------录音机shrpubw--------创建共享文件夹secpol.msc-----本地安全策略syskey---------系统加密,一旦加密就不能解开,保护windows xp系统的双重密码services.msc---本地服务设置Sndvol32-------音量控制程序sfc.exe--------系统文件检查器sfc /scannow---windows文件保护tsshutdn-------60秒倒计时关机命令tourstart------xp简介(安装完成后出现的漫游xp程序)taskmgr--------任务管理器eventvwr-------事件查看器eudcedit-------造字程序explorer-------打开资源管理器packager-------对象包装程序perfmon.msc----计算机性能监测程序progman--------程序管理器regedit.exe----注册表rsop.msc-------组策略结果集regedt32-------注册表编辑器rononce -p ----15秒关机regsvr32 /u *.dll----停止dll文件运行regsvr32 /u zipfldr.dll------取消ZIP支持cmd.exe--------CMD命令提示符chkdsk.exe-----Chkdsk磁盘检查certmgr.msc----证书管理实用程序calc-----------启动计算器charmap--------启动字符映射表cliconfg-------SQL SERVER 客户端网络实用程序Clipbrd--------剪贴板查看器conf-----------启动netmeetingcompmgmt.msc---计算机管理cleanmgr-------垃圾整理ciadv.msc------索引服务程序osk------------打开屏幕键盘odbcad32-------ODBC数据源管理器oobe/msoobe /a----检查XP是否激活lusrmgr.msc----本机用户和组logoff---------注销命令iexpress-------木马捆绑工具,系统自带Nslookup-------IP地址侦测器fsmgmt.msc-----共享文件夹管理器utilman--------辅助工具管理器sndrec32-------录音机Nslookup-------IP地址侦测器explorer-------打开资源管理器logoff---------注销命令tsshutdn-------60秒倒计时关机命令lusrmgr.msc----本机用户和组services.msc---本地服务设置notepad--------打开记事本cleanmgr-------垃圾整理 netstart messenger----开始信使服务compmgmt.msc---计算机管理 netstop messenger-----停止信使服务conf-----------启动netmeetingdvdplay--------DVD播放器charmap--------启动字符映射表diskmgmt.msc---磁盘管理实用程序calc-----------启动计算器dfrg.msc-------磁盘碎片整理程序chkdsk.exe-----Chkdsk磁盘检查devmgmt.msc--- 设备管理器regsvr32 /u *.dll----停止dll文件运行drwtsn32------ 系统医生regedt32-------注册表编辑器Msconfig.exe---系统配置实用程序rsop.msc-------组策略结果集mem.exe--------显示内存使用情况regedit.exe----注册表winchat--------XP自带局域网聊天progman--------程序管理器winmsd---------系统信息perfmon.msc----计算机性能监测程序winver---------检查Windows版本sfc /scannow-----扫描错误并复原taskmgr-----任务管理器-------------------------------------- wmimgmt.msc----打开windows管理体系结构wupdmgr--------windows更新程序winver---------检查Windows版本wmimgmt.msc----打开windows管理体系结构wupdmgr--------windows更新程序wscript--------windows脚本宿主设置write----------写字板winmsd-----系统信息wiaacmgr-------扫描仪和照相机向导winchat--------XP自带局域网聊天mem.exe--------显示内存使用情况Msconfig.exe---系统配置实用程序mplayer2-------简易widnows media playermspaint--------画图板mstsc----------远程桌面连接mplayer2-------媒体播放机magnify--------放大镜实用程序mmc------------打开控制台mobsync--------同步命令dxdiag---------检查DirectX信息drwtsn32------ 系统医生devmgmt.msc--- 设备管理器dfrg.msc-------磁盘碎片整理程序diskmgmt.msc---磁盘管理实用程序dcomcnfg-------打开系统组件服务ddeshare-------打开DDE共享设置dvdplay--------DVD播放器net stop messenger-----停止信使服务net start messenger----开始信使服务notepad--------打开记事本nslookup-------网络管理的工具向导ntbackup-------系统备份和还原narrator-------屏幕"讲述人"ntmsmgr.msc----移动存储管理器ntmsoprq.msc---移动存储管理员操作请求netstat -an----(TC)命令检查接口syncapp--------创建一个公文包sysedit--------系统配置编辑器sigverif-------文件签名验证程序sndrec32-------录音机shrpubw--------创建共享文件夹secpol.msc-----本地安全策略syskey---------系统加密,一旦加密就不能解开,保护windows xp系统的双重密码services.msc---本地服务设置Sndvol32-------音量控制程序sfc.exe--------系统文件检查器sfc /scannow---windows文件保护tsshutdn-------60秒倒计时关机命令tourstart------xp简介(安装完成后出现的漫游xp程序)taskmgr--------任务管理器eventvwr-------事件查看器eudcedit-------造字程序explorer-------打开资源管理器packager-------对象包装程序perfmon.msc----计算机性能监测程序progman--------程序管理器regedit.exe----注册表rsop.msc-------组策略结果集regedt32-------注册表编辑器rononce -p ----15秒关机regsvr32 /u *.dll----停止dll文件运行regsvr32 /u zipfldr.dll------取消ZIP支持cmd.exe--------CMD命令提示符chkdsk.exe-----Chkdsk磁盘检查certmgr.msc----证书管理实用程序calc-----------启动计算器charmap--------启动字符映射表cliconfg-------SQL SERVER 客户端网络实用程序Clipbrd--------剪贴板查看器conf-----------启动netmeetingcompmgmt.msc---计算机管理cleanmgr-------**整理ciadv.msc------索引服务程序osk------------打开屏幕键盘odbcad32-------ODBC数据源管理器oobe/msoobe /a----检查XP是否激活lusrmgr.msc----本机用户和组logoff---------注销命令iexpress-------木马捆绑工具,系统自带Nslookup-------IP地址侦测器fsmgmt.msc-----共享文件夹管理器utilman--------辅助工具管理器gpedit.msc-----组策略有几个不在不同的系统不能用。
EMD方法介绍及实证分析目录1.总体经验模式分解方法介绍 (1)1.1 EMD方法的引入 (1)1.2 EMD的基本理论和方法 (2)1.3 EEMD (3)2.实证分析 (4)2.1汇率算例分析 (4)2.2 基于EMD和GARCH模型的股价预测分析 (10)2.2.1 研究对象与数据选取 (10)2.2.2 EMD分解及分析 (10)2.2.3 自回归模型的拟合和预测 (15)2.2.4 GARCH模型的拟合和预测 (18)2.2.5 预测数据重组 (23)参考文献 (24)1.总体经验模式分解方法介绍1.1 EMD方法的引入近年来,小波变换(Wavelet Transformation , WT)理论在股票市场系统变量的多时间尺度分析与建模中取得了丰富的成果。
小波变换在时域和频域都具有良好的多分辨率分析能力,被誉为数学显微镜。
但小波变换实质上是一种窗口可调的傅立叶变换,其小波窗内的信号必须是平稳的,因而没有从根本上摆脱傅立叶分析的局限,小波变换虽然能够在频域和时域内同时得到较高的分辨率,但仍然存在一定的限制,这种限制通常会造成很多虚假的谐波,且小波基函数的选择对小波分解结果有显著的影响。
针对小波变换的不足,1998年,Huang等人提出来一种全新的多分辨率信号分析方法——经验模态分解(Empirical Mode Decomposition , EMD)。
EMD是基于信号局部特征时间尺度,从原信号中提取本征模态函数(Intrinsic Mode Function , IMF)。
在线性框架下基于EMD得到的Hilbert谱与小波谱具有相同的表现特性,而Hilbert谱在频域和时域内的分辨率都远高于小波谱,依此得到的分析结果可以更准确地反映系统原有的物理特性。
由于EMD方法比小波变换有更强的局部表现力,所以在处理非线性、非平稳信号时,EMD方法是一种更有效的方法,而金融时间序列(如股价、股价指数、收益率等)就是一类典型的非线性、非平稳时间序列。
EMD分解HHT变化matlab源代码第一篇:EMD分解HHT变化matlab源代码%function [] = UseEMD(x,t)function [] = UseEMD()%UNTITLED2 Summary of this function goes here %Detailed explanation goes hereN = 39;% # of data samplest = linspace(0,1,N);x=[20.3421.2520.6219.317.615.6815.4616.2717.9418.9718.7118.4719.1119.3118.6717.3216.2116.0916.0417.5419.320.1120.4221.3321.9422.1722.9723.0923.3923.9726.7129.3129.8230.3529.7827.9427.1727.8327.36];y=fft(x);plot(y,x);[imf,ort,nbits] = emd(x);emd_visu(x,t,imf,1);%¾ùÖµµÄƽ·½imfp2=mean(imf,2).^2%ƽ·½µÄ¾ùÖµimf2p=mean(imf.^2,2)%¸÷¸öIMFµÄ·½²îmse=imf2p-imfp2%·½²î°Ù·Ö±È£¬Ò²¾ÍÊÇ·½²î¹±Ï×ÂÊmseb=mse/sum(mse)*100%ÏÔʾ¸÷¸öIMFµÄ·½²îºÍ¹±Ï×ÂÊ[mse mseb]%¼ÆËãimfµÄÐÐÁÐάÊý[m,n] = size(imf)hx = hilbert(x)xf = real(hx);xi = imag(hx);%¼ÆËã˲ʱÕñ·ùA=sqrt(xf.^2 + xi.^2);figure(4)plot(t,A);title('˲ʱÕñ·ù')%¼ÆËã˲ʱÏàλp=angle(hx);figure(5);plot(t,p);title('˲ʱÏàλ')%¼ÆËã˲ʱƵÂÊ%dt = diff(t);%dx=diff(p);%sp = dx./dt;%figure(6);%plot(t(1:N-1),sp);title('˲ʱƵÂÊ') %imf1µÄhilbert±ä»»xn1 = hilbert(imf(1,:));xr1 = real(xn1);xi1 = imag(xn1);%imf1µÄ˲ʱÕñ·ùA1=sqrt(xr1.^2+xi1.^2);figure(7);subplot(2,1,1);plot(t,A1);xlabel('ʱ¼äª¡tª¢');ylabel('˲ʱÕñ·ù');title('imf1')%imf2µÄHilbert±ä»»xn2=hilbert(imf(2,:));xr2=real(xn2);xi2=imag(xn2);%imf2µÄ˲ʱÕñ·ùA2=sqrt(xr2.^2+xi2.^2); subplot(2,1,2);plot(t,A2);xlabel('ʱ¼äª¡tª¢');ylabel('˲ʱÕñ·ù');title('imf2')%imf1µÄ˲ʱÏàλP1=atan2(xi1,xr1);figure(8);subplot(2,1,1);plot(t,P1);xlabel('ʱ¼äª¡tª¢');ylabel('˲ʱÏàλ');title('imf1')%imf2µÄ˲ʱÏàλP2=atan2(xi2,xr2);subplot(2,1,2);plot(t,P2);xlabel('ʱ¼äª¡tª¢');ylabel('˲ʱÏàλ');title('imf2')%imf1˲ʱƵÂÊxh1=unwrap(P1);fs=1000;xhd1=fs*diff(xh1)/(2*pi);figure(9);subplot(2,1,1);plot(t(1:38),xhd1);xlabel('ʱ¼äª¡tª¢');ylabel('˲ʱƵÂÊ');title('imf1')%imf2˲ʱƵÂÊxh2=unwrap(P2);fs=1000;xhd2=fs*diff(xh2)/(2*pi);subplot(2,1,2);plot(t(1:38),xhd2);xlabel('ʱ¼äª¡tª¢');ylabel('˲ʱƵÂÊ');title('imf2')%¼ÆËãHHTµÄʱƵÆ×[A, fa, tt] = hhspectrum(imf);if size(imf,1)> 1[A,fa,tt] = hhspectrum(imf(1:end-1, :));else [A,fa,tt] = hhspectrum(imf);end[E,tt1]=toimage(A,fa,tt,length(tt));disp_hhs(E,[],fs)%¶þάͼÐÎÏÔʾHHTÊÓÆµÆ×ª£EÊÇÇóµÃµÄHHTÆ×for i=1:m-1faa=fa(i,:);[FA,TT1]=meshgrid(faa,tt1);%ÈýάͼÏÔʾHHTʱƵͼfigur e(11);surf(FA,TT1,E)title('HHTʱƵÆ×ÈýάÏÔʾ')end%»-±ß¼ÊÆ×E=flipud(E);for k=1:size(E,1)bjp(k)=sum(E(k,:))*1/fs;endf =(0:N-3)/N*(fs/2);figure(12)plot(f,bjp);end注意:matlab中需加载instfreq.m文件,从网上可下到第二篇:LU分解MatLab算法分析最近矩阵分析老师出了一道题目作为作业,是一道程序题,题目是对A矩阵做LU分解,要求能对A实现PA=LU的分解,并最终输出L,U,P矩阵。
经验模态分解 (emd)方法一、EMD方法概述经验模态分解(EMD)是一种用于信号分解和特征提取的自适应方法,它可以将一个复杂的信号分解为一系列本征模态函数(IMF)的叠加。
IMF是具有自适应频率的函数,它们能够准确地描述信号的局部特征。
EMD方法不需要先验知识和基函数的选择,因此在信号分析和图像处理领域中得到了广泛应用。
二、EMD方法的基本原理EMD方法的基本原理是将信号分解为一组IMF,并且每个IMF均满足以下两个条件:1)在整个信号上,它的正负波动次数应该相等或相差不超过一个;2)在任意一点上,它的均值应该为零。
通过迭代处理,可以得到一系列IMF,并且每一次迭代都能更好地逼近原始信号。
三、EMD方法的步骤EMD方法的具体步骤如下:1)将原始信号进行局部极大值和极小值的插值,得到上、下包络线;2)计算信号的局部均值;3)将信号减去局部均值,得到一次IMF分量;4)判断分量是否满足IMF的两个条件,如果满足则停止,否则将分量作为新的信号进行迭代处理,直到满足条件为止。
四、EMD方法在信号分析中的应用EMD方法在信号分析中有着广泛的应用。
例如,在地震学中,可以利用EMD方法对地震信号进行分解,提取出不同频率范围的地震波,从而对地震波进行特征提取和识别。
另外,在生物医学信号处理中,EMD方法可以应用于心电图信号的分解和特征提取,有助于对心脏疾病进行诊断和监测。
五、EMD方法在图像处理中的应用EMD方法在图像处理中也有着广泛的应用。
例如,在图像压缩领域,可以利用EMD方法对图像进行分解,提取出不同频率的图像分量,从而实现对图像的压缩和重构。
此外,在图像去噪和边缘检测中,EMD方法也能够有效地提取出图像的局部特征信息,有助于准确地去除噪声和检测图像边缘。
六、EMD方法的优缺点EMD方法具有以下优点:1)能够自适应地分解信号,无需先验知识和基函数的选择;2)能够准确地描述信号的局部特征;3)能够处理非线性和非平稳信号。
EEMD:总体经验模态分解,一种改进的EMD方法。
步骤如下:1)通过给目标信号加上一组白噪声来获得一个总体;2)对进行EMD分解,得到各个IMF分量;3)给目标信号加入不同的白噪,重复以上两步;4)分解后得到各自的IMF分量组;5)取相应IMF的均值作为最终的IMF组;6)取相应的剩余分量的均值作为最终的IMF组。
一.EEMD模块%function allmode=eemd(Y,Nstd,NE)%% This is an EMD/EEMD program%% INPUT:% Y: Inputted data;1-d data only% Nstd: ratio of the standard deviation of the added noise and that of Y;% NE: Ensemble number for the 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, ...% m are the IMFs from high to low frequency, and comlumn (m+1) is the% residual (over all trend).%% 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.(for EMD Nstd=0,NE=1)% This code limited sift number=10 ,the stoppage criteria can't change.应该注意,当Nstd设置为零,NE设置为1时,该项目退化EMD计划。
(EMD Nstd = 0 NE = 1) 这段代码筛选有限数量= 10,中止条件不能改变。
%% References:% Wu, Z., and N. E Huang (2008),% Ensemble Empirical Mode Decomposition: a noise-assisted data analysis method.% Advances in Adaptive Data Analysis. V ol.1, No.1. 1-41.%% code writer: Zhaohua Wu.% footnote:S.C.Su 2009/03/04%% There are three loops in this code coupled together.% 1.read data, find out standard deviation ,devide all data by std% 2.evaluate TNM as total IMF number--eq1.% TNM2=TNM+2,original data and residual included in TNM2% assign 0 to TNM2 matrix% 3.Do EEMD NE times-------------------------------------------------------------loop EEMD start % 4.add noise% 5.give initial values before sift% 6.start to find an IMF------------------------------------------------IMF loop start% 7.sift 10 times to get IMF--------------------------sift loop start and end% 8.after 10 times sift --we got IMF% 9.subtract IMF from data ,and let the residual to find next IMF by loop% 6.after having all the IMFs---------------------------------------------IMF loop end% 9.after TNM IMFs ,the residual xend is over all trend% 3.Sum up NE decomposition result-------------------------------------------------loop EEMD end % 10.Devide EEMD summation by NE,std be multiply back to data%% Association: no% this function ususally used for doing 1-D EEMD with fixed% stoppage criteria independently.%% Concerned function: extrema.m% above mentioned m file must be put togetherfunction allmode=eemd(Y,Nstd,NE)%part1.read data, find out standard deviation ,devide all data by stdxsize=length(Y);dd=1:1:xsize;Ystd=std(Y);Y=Y/Ystd;%part2.evaluate TNM as total IMF number,ssign 0 to TNM2 matrixTNM=fix(log2(xsize))-1;TNM2=TNM+2;for kk=1:1:TNM2,for ii=1:1:xsize,allmode(ii,kk)=0.0;endend%part3 Do EEMD -----EEMD loop startfor iii=1:1:NE, %EEMD loop -NE times EMD sum together%part4 --Add noise to original data,we have X1for i=1:xsize,temp=randn(1,1)*Nstd;X1(i)=Y(i)+temp;end%part4 --assign original data in the first columnfor jj=1:1:xsize,mode(jj,1) = Y(jj);end%part5--give initial 0 to xorigin and xendxorigin = X1;xend = xorigin;%part6--start to find an IMF-----IMF loop startnmode = 1;while nmode <= TNM,xstart = xend; %last loop value assign to new iteration loop%xstart -loop start dataiter = 1; %loop index initial value%part7--sift 10 times to get IMF---sift loop startwhile iter<=10,[spmax, spmin, flag]=extrema(xstart); %call function extrema%the usage of spline ,please see part11.upper= spline(spmax(:,1),spmax(:,2),dd); %upper spline bound of this siftlower= spline(spmin(:,1),spmin(:,2),dd); %lower spline bound of this siftmean_ul = (upper + lower)/2;%spline mean of upper and lowerxstart = xstart - mean_ul;%extract spline mean from Xstartiter = iter +1;end%part7--sift 10 times to get IMF---sift loop end%part8--subtract IMF from data ,then let the residual xend to start to find next IMF xend = xend - xstart;nmode=nmode+1;%part9--after sift 10 times,that xstart is this time IMFfor jj=1:1:xsize,mode(jj,nmode) = xstart(jj);endend%part6--start to find an IMF-----IMF loop end%part 10--after gotten all(TNM) IMFs ,the residual xend is over all trend% put them in the last columnfor jj=1:1:xsize,mode(jj,nmode+1)=xend(jj);end%after part 10 ,original +TNM-IMF+overall trend ---those are all in modeallmode=allmode+mode;end%part3 Do EEMD -----EEMD loop end%part10--devide EEMD summation by NE,std be multiply back to dataallmode=allmode/NE;allmode=allmode*Ystd;%part11--the syntax of the matlab function spline%yy= spline(x,y,xx); this means%x and y are matrixs of n1 points ,use n1 set (x,y) to form the cubic spline%xx and yy are matrixs of n2 points,we want know the spline value yy(y-axis) in the xx (x-axis)position%after the spline is formed by n1 points ,find coordinate value on the spline for [xx,yy] --n2 position.二.Extrema模块% This is a utility program for significance test.%% function [spmax, spmin, flag]= extrema(in_data)%% INPUT:% in_data: Inputted data, a time series to be sifted;% OUTPUT:% spmax: The locations (col 1) of the maxima and its corresponding % values (col 2)% spmin: The locations (col 1) of the minima and its corresponding % values (col 2)%% References can be found in the "Reference" section.%% The code is prepared by Zhaohua Wu. For questions, please read the "Q&A"section or% contact% zhwu@%function [spmax, spmin, flag]= extrema(in_data)flag=1;dsize=length(in_data);spmax(1,1) = 1;spmax(1,2) = in_data(1);jj=2;kk=2;while jj<dsize,if ( in_data(jj-1)<=in_data(jj) & in_data(jj)>=in_data(jj+1) ) spmax(kk,1) = jj;spmax(kk,2) = in_data (jj);kk = kk+1;endjj=jj+1;endspmax(kk,1)=dsize;spmax(kk,2)=in_data(dsize);if kk>=4slope1=(spmax(2,2)-spmax(3,2))/(spmax(2,1)-spmax(3,1));tmp1=slope1*(spmax(1,1)-spmax(2,1))+spmax(2,2);if tmp1>spmax(1,2)spmax(1,2)=tmp1;endslope2=(spmax(kk-1,2)-spmax(kk-2,2))/(spmax(kk-1,1)-spmax(kk-2,1)); tmp2=slope2*(spmax(kk,1)-spmax(kk-1,1))+spmax(kk-1,2);if tmp2>spmax(kk,2)spmax(kk,2)=tmp2;endelseflag=-1;endmsize=size(in_data);dsize=max(msize);xsize=dsize/3;xsize2=2*xsize;spmin(1,1) = 1;spmin(1,2) = in_data(1);jj=2;kk=2;while jj<dsize,if ( in_data(jj-1)>=in_data(jj) & in_data(jj)<=in_data(jj+1))spmin(kk,1) = jj;spmin(kk,2) = in_data (jj);kk = kk+1;endjj=jj+1;endspmin(kk,1)=dsize;spmin(kk,2)=in_data(dsize);if kk>=4slope1=(spmin(2,2)-spmin(3,2))/(spmin(2,1)-spmin(3,1));tmp1=slope1*(spmin(1,1)-spmin(2,1))+spmin(2,2);if tmp1<spmin(1,2)spmin(1,2)=tmp1;endslope2=(spmin(kk-1,2)-spmin(kk-2,2))/(spmin(kk-1,1)-spmin(kk-2,1)); tmp2=slope2*(spmin(kk,1)-spmin(kk-1,1))+spmin(kk-1,2);if tmp2<spmin(kk,2)spmin(kk,2)=tmp2;endelseflag=-1;endflag=1;三.例子clc;clear;t=0:0.000325:0.32;x=10*sin(2*pi*25*t)+sin(2*pi*75*t);%原始信号两组分能量相差太大发生模态混叠dt=0.000325;imfe_x=eemd(x,0.1,100); subplot(411)plot(t,x);ylabel('x(t)');subplot(412)plot(t,imfe_x(:,1)); ylabel('c1')subplot(413)plot(t,imfe_x(:,2)); ylabel('c2')subplot(414)plot(t,imfe_x(:,3)); ylabel('c3')。