EMD算法的matlab程序介绍解析

合集下载

matlab练习程序（EM算法）

matlab练习程序（EM算法）最⼤期望算法（Expectation-maximization algorithm，⼜译为期望最⼤化算法），是在概率模型中寻找参数最⼤似然估计或者最⼤后验估计的算法，其中概率模型依赖于⽆法观测的隐性变量。

最⼤期望算法经过两个步骤交替进⾏计算：第⼀步是计算期望（E），利⽤对隐藏变量的现有估计值，计算其最⼤似然估计值；第⼆步是最⼤化（M），最⼤化在E步上求得的最⼤似然值来计算参数的值。

M步上找到的参数估计值被⽤于下⼀个E步计算中，这个过程不断交替进⾏。

下⾯⽤EM算法对⾼斯混合模型进⾏聚类，该算法不仅能够给出聚类后数据的均值，还能够给出协⽅差。

代码如下：clear all;close all;clc;mu1=[00];S1=[0.80.1];data1=mvnrnd(mu1,S1,1000);plot(data1(:,1),data1(:,2),'r.');hold on;mu2=[24];S2=[0.41.3];data2=mvnrnd(mu2,S2,1000);plot(data2(:,1),data2(:,2),'g.');mu3=[-23];S3=[2.41.3];data3=mvnrnd(mu3,S3,1000);plot(data3(:,1),data3(:,2),'b.');%利⽤EM算法对⾼斯混合模型聚类data=[data1;data2;data3];mu{1} = rand(1,2);mu{2} = rand(1,2);mu{3} = rand(1,2);sigma{1} = rand(1,2);sigma{2} = rand(1,2);sigma{3} = rand(1,2);p = [0.30.40.4];w=zeros(length(data),3);for i=1:1000%E-stepfor j=1:3w(:,j) = p(j)*mvnpdf(data,mu{j},sigma{j});endw = w./repmat(sum(w,2),[13]);%M-stepfor j=1:3mu{j} = w(:,j)'* data / sum(w(:,j));sigma{j} = sqrt(w(:,j)'*((data-mu{j}).*(data-mu{j})) / sum(w(:,j)));endp = sum(w) / length(data);endfigure;w = uint8(w);data1 = data(w(:,1)==1,:);data2 = data(w(:,2)==1,:);data3 = data(w(:,3)==1,:);plot(data1(:,1),data1(:,2),'r.');hold on;plot(data2(:,1),data2(:,2),'g.');plot(data3(:,1),data3(:,2),'b.');结果：原始⾼斯数据：聚类后：。

emd 算法原理

emd 算法原理
EMD算法，即经验模态分解算法，是一种能够将任意信号分解为一组固有振动模态的非平稳信号分解方法。

该算法的基本思想是将待分解信号视为一组固有振动模态的叠加，每个模态都是具有不同频率和振幅的信号。

通过不断迭代，可以逐步将信号分解为多个固有振动模态。

EMD算法的核心是求解局部极值点，从而确定每个固有振动模态的上下包络线。

具体而言，EMD算法分为以下几个步骤：
1. 将信号拟合为一条直线，并计算信号与该直线的差值。

2. 找到信号的所有局部极值点，包括极大值和极小值。

3. 将所有局部极值点连接成一组上下包络线，形成一个固有振动模态。

4. 将信号减去该固有振动模态，得到一个新的信号，并重复步骤1-3，直到该信号可以被分解为一组固有振动模态。

EMD算法的优点在于可以适应非线性和非平稳信号，但其缺点在于计算量较大，计算时间较长。

因此，在实际应用中需要谨慎选择算法参数，并注意算法的稳定性和可靠性。

- 1 -。

matlab源码-GMM的EM算法实现

法的实现与收敛性证明进行了详细说明。

本文主要针对如何用EM算法在混合高斯模型下进行聚类进行代码上的分析说明。

1. GMM模型：每个GMM 由K 个Gaussian 分布组成，每个Gaussian 称为一个“Comp onent〞，这些Component 线性加成在一起就组成了GMM 的概率密度函数：根据上面的式子，如果我们要从GMM 的分布中随机地取一个点的话，实际上可以分为两步：首先随机地在这K个Gaussian Component 之中选一个，每个Component 被选中的概率实际上就是它的系数pi(k) ，选中了Compone nt 之后，再单独地考虑从这个Component 的分布中选取一个点就可以了──这里已经回到了普通的Gaussian 分布，转化为了的问题。

那么如何用GMM 来做clustering 呢？其实很简单，现在我们有了数据，假定它们是由GMM 生成出来的，那么我们只要根据数据推出GMM 的概率分布来就可以了，然后GMM 的K 个Component 实际上就对应了K 个clus ter 了。

根据数据来推算概率密度通常被称作density estimation ，特别地，当我们在〔或假定〕了概率密度函数的形式，而要估计其中的参数的过程被称作“参数估计〞。

2. 参数与似然函数：现在假设我们有N 个数据点，并假设它们服从某个分布〔记作p(x) 〕，现在要确定里面的一些参数的值，例如，在GMM 中，我们就需要确定影响因子pi (k)、各类均值pMiu(k) 和各类协方差pSigma(k) 这些参数。

我们的想法是，找到这样一组参数，它所确定的概率分布生成这些给定的数据点的概率最大，而这个概率实际上就等于，我们把这个乘积称作似然函数(Likelihood Function)。

通常单个点的概率都很小，许多很小的数字相乘起来在计算机里很容易造成浮点数下溢，因此我们通常会对其取对数，把乘积变成加和，得到log-likelihood function 。

(完整word版)EM算法matlab程序

X=zeros（600，2);X(1：200,：）= normrnd（0,1，200,2）；X(201:400,:）= normrnd(0，2,200，2);X（401:600，:）= normrnd（0,3,200,2）；[W，M,V，L] = EM_GM(X，3,［]，[],1，[］）下面是程序源码:打印帮助function[W，M,V，L］= EM_GM（X，k,ltol,maxiter,pflag,Init)％[W，M,V,L］= EM_GM(X，k，ltol,maxiter，pflag,Init）％% EM algorithm for k multidimensional Gaussian mixture estimation％％Inputs：％X(n，d）- input data，n=number of observations，d=dimension of variable % k - maximum number of Gaussian components allowed％ltol - percentage of the log likelihood difference between 2 iterations ([] for none) % maxiter — maximum number of iteration allowed (［] for none）％pflag — 1 for plotting GM for 1D or 2D cases only，0 otherwise （[］for none) % Init — structure of initial W，M, V: Init。

W，Init。

M, Init。

V （［] for none）％％Ouputs：% W(1,k）— estimated weights of GM% M（d，k) — estimated mean vectors of GM％V（d,d,k) — estimated covariance matrices of GM％L — log likelihood of estimates％％Written by％Patrick P。

经验模态分解 (emd) 方法划分层序

经验模态分解(emd) 方法划分层序摘要：1.经验模态分解（EMD）简介2.EMD方法在划分层序中的应用3.具体实施步骤与案例分析4.总结与展望正文：一、经验模态分解（EMD）简介经验模态分解（Empirical Mode Decomposition，简称EMD）是一种自适应的信号分解方法，由Norden E.Huang等人于1998年首次提出。

该方法主要通过对信号进行局部均值拟合，将原始信号分解为多个本征模态函数（Intrinsic Mode Functions，简称IMFs）。

本征模态函数代表了信号在不同时间尺度上的特征，从而实现了信号的时频分析。

二、EMD方法在划分层序中的应用1.地质勘探：EMD方法在地质勘探领域具有广泛应用，如地层划分、岩性识别等。

通过对地震、测井等原始信号进行经验模态分解，可以获取各个本征模态函数，进一步分析地层的结构和成分。

2.工程监测：在工程领域，EMD方法可用于结构健康监测、故障诊断等。

例如，对桥梁、建筑物等结构物的振动信号进行经验模态分解，可以识别出结构的损伤程度和位置。

3.生物医学：EMD方法在生物医学领域也有广泛应用，如心电信号分析、脑电信号分析等。

通过对生物信号进行经验模态分解，可以获取有价值的信息，有助于疾病的诊断和治疗。

4.金融分析：EMD方法在金融领域也有显著的应用，如股票价格预测、汇率预测等。

通过对金融时间序列数据进行经验模态分解，可以分析市场的波动特征，为投资者提供参考。

三、具体实施步骤与案例分析1.数据预处理：对原始信号进行去噪、滤波等预处理，以消除信号中的噪声和干扰。

2.经验模态分解：利用EMD方法将预处理后的信号分解为多个本征模态函数。

3.划分层序：根据本征模态函数的特性，对信号进行分层。

例如，可以按照频率、能量等特征将本征模态函数划分为不同层次。

4.分析与诊断：对划分的层次进行进一步分析，提取有价值的信息，实现信号的诊断和分析。

案例分析：以地质勘探为例，经验模态分解可以应用于地震信号的处理，划分出不同频率的本征模态函数。

经验模态分解算法

经验模态分解算法
EMD算法的步骤如下：
1.将要分解的信号称为原始信号，记为x(t)。

2.寻找x(t)的极大值点和极小值点，这些点将原始信号分为一系列小段。

3.对每个小段进行插值，使均匀分布的数据点可以拟合出这个小段。

4. 利用Cubic Spline插值法或其他插值方法找到一个包络线，该包络线连接这些插值点的极大值点和极小值点。

即为信号中的一条上包络线和一条下包络线。

5.计算出平均值函数m(t)=(上包络线+下包络线)/2
6.计算x(t)与m(t)的差值d(t)=x(t)-m(t)。

7.如果d(t)是一条IMF，则终止算法；否则将d(t)作为新的原始信号，重复步骤2-6
8.将计算出的IMF组合起来，得到原始信号x(t)的EMD分解结果。

EMD算法的特点是对信号进行自适应分解，能够捕捉到不同频率的局部特征。

它不需要提前设定基函数或者滤波器，而是根据信号中的局部特征自动适应地生成各个IMF。

因此，EMD算法在信号处理领域中得到了广泛应用，如地震信号分析、生物信号处理等。

然而，EMD算法也存在一些问题。

其中最主要的问题是固有模态函数的提取过程中可能出现模态混叠的情况，即两个或多个IMF的频率相似且在一些区间内相互重叠，使得提取的IMF不纯粹。

为了克服这个问题，研
究者们提出了一些改进的EMD算法，如快速EMD、改进的EMD等。

这些改进方法在一定程度上提高了EMD算法的可靠性和稳定性。

总之，经验模态分解算法是一种有效的信号分解方法，能够提供信号的局部特征表示。

它在很多领域有广泛的应用，但仍然需要进一步的研究和改进，以提高其分解效果和精度。

EM算法matlab程序

X=zeros(600,2);X(1:200,:) = normrnd(0,1,200,2);X(201:400,:) = normrnd(0,2,200,2);X(401:600,:) = normrnd(0,3,200,2);[W,M,V,L] = EM_GM(X,3,[],[],1,[])下面是程序源码：打印帮助function[W,M,V,L] = EM_GM(X,k,ltol,maxiter,pflag,Init)% [W,M,V,L] = EM_GM(X,k,ltol,maxiter,pflag,Init)%% EM algorithm for k multidimensional Gaussian mixture estimation%% Inputs:% X(n,d) - input data, n=number of observations, d=dimension of variable% k - maximum number of Gaussian components allowed% ltol - percentage of the log likelihood difference between 2 iterations ([] for none) % maxiter - maximum number of iteration allowed ([] for none)% pflag - 1 for plotting GM for 1D or 2D cases only, 0 otherwise ([] for none)% Init - structure of initial W, M, V: Init.W, Init.M, Init.V ([] for none)%% Ouputs:% W(1,k) - estimated weights of GM% M(d,k) - estimated mean vectors of GM% V(d,d,k) - estimated covariance matrices of GM% L - log likelihood of estimates%% Written by% Patrick P. C. Tsui,% PAMI research group% Department of Electrical and Computer Engineering% University of Waterloo,% March, 2006%%%%% Validate inputs %%%%ifnargin <= 1,disp('EM_GM must have at least 2 inputs: X,k!/n')returnelseifnargin == 2,ltol = 0.1; maxiter = 1000; pflag = 0; Init = [];err_X = Verify_X(X);err_k = Verify_k(k);iferr_X | err_k,return;endelseifnargin == 3,maxiter = 1000; pflag = 0; Init = [];err_X = Verify_X(X);err_k = Verify_k(k);[ltol,err_ltol] = Verify_ltol(ltol);iferr_X | err_k | err_ltol,return;endelseifnargin == 4,pflag = 0; Init = [];err_X = Verify_X(X);err_k = Verify_k(k);[ltol,err_ltol] = Verify_ltol(ltol);[maxiter,err_maxiter] = Verify_maxiter(maxiter);iferr_X | err_k | err_ltol | err_maxiter,return;endelseifnargin == 5,Init = [];err_X = Verify_X(X);err_k = Verify_k(k);[ltol,err_ltol] = Verify_ltol(ltol);[maxiter,err_maxiter] = Verify_maxiter(maxiter);[pflag,err_pflag] = Verify_pflag(pflag);iferr_X | err_k | err_ltol | err_maxiter | err_pflag,return;end elseifnargin == 6,err_X = Verify_X(X);err_k = Verify_k(k);[ltol,err_ltol] = Verify_ltol(ltol);[maxiter,err_maxiter] = Verify_maxiter(maxiter);[pflag,err_pflag] = Verify_pflag(pflag);[Init,err_Init]=Verify_Init(Init);iferr_X | err_k | err_ltol | err_maxiter | err_pflag | err_Init,return;end elsedisp('EM_GM must have 2 to 6 inputs!');returnend%%%% Initialize W, M, V,L %%%%t = cputime;ifisempty(Init),[W,M,V] = Init_EM(X,k); L = 0;elseW = Init.W;M = Init.M;V = Init.V;endLn = Likelihood(X,k,W,M,V);% Initialize log likelihoodLo = 2*Ln;%%%% EM algorithm %%%%niter = 0;while(abs(100*(Ln-Lo)/Lo)>ltol) & (niter<=maxiter),E = Expectation(X,k,W,M,V);% E-step[W,M,V] = Maximization(X,k,E); % M-stepLo = Ln;Ln = Likelihood(X,k,W,M,V);niter = niter + 1;endL = Ln;%%%% Plot 1D or 2D %%%%ifpflag==1,[n,d] = size(X);ifd>2,disp('Can only plot 1 or 2 dimensional applications!/n');elsePlot_GM(X,k,W,M,V);endelapsed_time = sprintf('CPU time used for EM_GM: %5.2fs',cputime-t); disp(elapsed_time);disp(sprintf('Number of iterations: %d',niter-1));end%%%%%%%%%%%%%%%%%%%%%%%%%% End of EM_GM %%%%%%%%%%%%%%%%%%%%%%%%%%functionE = Expectation(X,k,W,M,V)[n,d] = size(X);a = (2*pi)^(0.5*d);S = zeros(1,k);iV = zeros(d,d,k);forj=1:k,ifV(:,:,j)==zeros(d,d), V(:,:,j)=ones(d,d)*eps;endS(j) = sqrt(det(V(:,:,j)));iV(:,:,j) = inv(V(:,:,j));endE = zeros(n,k);fori=1:n,forj=1:k,dXM = X(i,:)'-M(:,j);pl = exp(-0.5*dXM'*iV(:,:,j)*dXM)/(a*S(j));E(i,j) = W(j)*pl;endE(i,:) = E(i,:)/sum(E(i,:));end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End of Expectation %%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%function[W,M,V] = Maximization(X,k,E)[n,d] = size(X);W = zeros(1,k); M = zeros(d,k);V = zeros(d,d,k);fori=1:k, % Compute weightsforj=1:n,W(i) = W(i) + E(j,i);M(:,i) = M(:,i) + E(j,i)*X(j,:)';endM(:,i) = M(:,i)/W(i);endfori=1:k,forj=1:n,dXM = X(j,:)'-M(:,i);V(:,:,i) = V(:,:,i) + E(j,i)*dXM*dXM';endV(:,:,i) = V(:,:,i)/W(i);endW = W/n;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End of Maximization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%functionL = Likelihood(X,k,W,M,V)% Compute L based on K. V. Mardia, "Multivariate Analysis", Academic Press, 1979, PP. 96-97 % to enchance computational speed[n,d] = size(X);U = mean(X)';S = cov(X);L = 0;fori=1:k,iV = inv(V(:,:,i));L = L + W(i)*(-0.5*n*log(det(2*pi*V(:,:,i))) ...-0.5*(n-1)*(trace(iV*S)+(U-M(:,i))'*iV*(U-M(:,i))));end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End of Likelihood %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%functionerr_X = Verify_X(X)err_X = 1;[n,d] = size(X);ifn<d,disp('Input data must be n x d!/n');returnenderr_X = 0;%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End of Verify_X %%%%%%%%%%%%%%%%%%%%%%%%%%%%%functionerr_k = Verify_k(k)err_k = 1;if~isnumeric(k) | ~isreal(k) | k<1,disp('k must be a real integer >= 1!/n');returnenderr_k = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%% End of Verify_k %%%% %%%%%%%%%%%%%%%%%%%%%%%%%function[ltol,err_ltol] = Verify_ltol(ltol)err_ltol = 1;ifisempty(ltol),ltol = 0.1;elseif~isreal(ltol) | ltol<=0,disp('ltol must be a positive real number!');returnenderr_ltol = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End of Verify_ltol %%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%function[maxiter,err_maxiter] = Verify_maxiter(maxiter) err_maxiter = 1;ifisempty(maxiter),maxiter = 1000;elseif~isreal(maxiter) | maxiter<=0,disp('ltol must be a positive real number!');returnenderr_maxiter = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End of Verify_maxiter %%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[pflag,err_pflag] = Verify_pflag(pflag)err_pflag = 1;ifisempty(pflag),pflag = 0;elseifpflag~=0 & pflag~=1,disp('Plot flag must be either 0 or 1!/n');returnenderr_pflag = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End of Verify_pflag %%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[Init,err_Init] = Verify_Init(Init)err_Init = 1;ifisempty(Init),% Do nothing;elseifisstruct(Init),[Wd,Wk] = size(Init.W);[Md,Mk] = size(Init.M);[Vd1,Vd2,Vk] = size(Init.V);ifWk~=Mk | Wk~=Vk | Mk~=Vk,disp('k in Init.W(1,k), Init.M(d,k) and Init.V(d,d,k) must equal!/n')returnendifMd~=Vd1 | Md~=Vd2 | Vd1~=Vd2,disp('d in Init.W(1,k), Init.M(d,k) and Init.V(d,d,k) must equal!/n')returnendelsedisp('Init must be a structure: W(1,k), M(d,k), V(d,d,k) or []!');returnenderr_Init = 0;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End of Verify_Init %%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%function[W,M,V] = Init_EM(X,k)[n,d] = size(X);[Ci,C] = kmeans(X,k,'Start','cluster', ...'Maxiter',100, ...'EmptyAction','drop', ...'Display','off');% Ci(nx1) - cluster indeices; C(k,d) - cluster centroid (i.e. mean) whilesum(isnan(C))>0,[Ci,C] = kmeans(X,k,'Start','cluster', ...'Maxiter',100, ...'EmptyAction','drop', ...'Display','off');endM = C';Vp = repmat(struct('count',0,'X',zeros(n,d)),1,k);fori=1:n,% Separate cluster pointsVp(Ci(i)).count = Vp(Ci(i)).count + 1;Vp(Ci(i)).X(Vp(Ci(i)).count,:) = X(i,:);endV = zeros(d,d,k);fori=1:k,W(i) = Vp(i).count/n;V(:,:,i) = cov(Vp(i).X(1:Vp(i).count,:));end%%%%%%%%%%%%%%%%%%%%%%%%%%%% End of Init_EM %%%%%%%%%%%%%%%%%%%%%%%%%%%%functionPlot_GM(X,k,W,M,V)[n,d] = size(X);ifd>2,disp('Can only plot 1 or 2 dimensional applications!/n'); returnendS = zeros(d,k);R1 = zeros(d,k);R2 = zeros(d,k);fori=1:k, % Determine plot range as 4 x standard deviations S(:,i) = sqrt(diag(V(:,:,i)));R1(:,i) = M(:,i)-4*S(:,i);R2(:,i) = M(:,i)+4*S(:,i);endRmin = min(min(R1));Rmax = max(max(R2));R = [Rmin:0.001*(Rmax-Rmin):Rmax];clf, hold onifd==1,Q = zeros(size(R));fori=1:k,P = W(i)*normpdf(R,M(:,i),sqrt(V(:,:,i)));Q = Q + P;plot(R,P,'r-'); grid on,endplot(R,Q,'k-');xlabel('X');ylabel('Probability density');else% d==2plot(X(:,1),X(:,2),'r.');fori=1:k,Plot_Std_Ellipse(M(:,i),V(:,:,i));endxlabel('1^{st} dimension');ylabel('2^{nd} dimension');axis([Rmin Rmax Rmin Rmax])endtitle('Gaussian Mixture estimated by EM'); %%%%%%%%%%%%%%%%%%%%%%%%%%%% End of Plot_GM %%%% %%%%%%%%%%%%%%%%%%%%%%%%functionPlot_Std_Ellipse(M,V)[Ev,D] = eig(V);d = length(M);ifV(:,:)==zeros(d,d),V(:,:) = ones(d,d)*eps;endiV = inv(V);% Find the larger projectionP = [1,0;0,0]; % X-axis projection operatorP1 = P * 2*sqrt(D(1,1)) * Ev(:,1);P2 = P * 2*sqrt(D(2,2)) * Ev(:,2);ifabs(P1(1)) >= abs(P2(1)),Plen = P1(1);elsePlen = P2(1);endcount = 1;step = 0.001*Plen;Contour1 = zeros(2001,2);Contour2 = zeros(2001,2);forx = -Plen:step:Plen,a = iV(2,2);b = x * (iV(1,2)+iV(2,1));c = (x^2) * iV(1,1) - 1;Root1 = (-b + sqrt(b^2 - 4*a*c))/(2*a);Root2 = (-b - sqrt(b^2 - 4*a*c))/(2*a);ifisreal(Root1),Contour1(count,:) = [x,Root1] + M';Contour2(count,:) = [x,Root2] + M';count = count + 1;endendContour1 = Contour1(1:count-1,:);Contour2 = [Contour1(1,:);Contour2(1:count-1,:);Contour1(count-1,:)]; plot(M(1),M(2),'k+');plot(Contour1(:,1),Contour1(:,2),'k-');plot(Contour2(:,1),Contour2(:,2),'k-'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End of Plot_Std_Ellipse %%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%。

EMD分解HHT变化matlab源代码

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方法概述经验模态分解(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)能够处理非线性和非平稳信号。

EMD算法

EMD算法
（一）数据信号处理方面
对数据信号进行EMD分解就是为了获得本征模函数。
一个本征模函数必须满足以下两个条件： ⑴函数在整个时间范围内，局部极值点和过零点的数目必须相等，或最多相差一个； ⑵在任意时刻点，局部最大值的包络（上包络线）和局部最小值的包络（下包络线）平均必须为零。
二、基本原理
一、简介
EMD算法
经验模态分解（Empirical Mode Decomposition，简称EMD )方法是由美国NASA的黄锷博士提出的一种信号分析方法.它依据数据自身的时间尺度特征来进行信号分解，无须预先设定任何基函数。
适合于分析非线性、非平稳信号序列，具有很高的信噪比。
二、基本原理
接收端提取隐秘信息时,仅需将(g1,g2)代入到公式(2-2)中计算出f,f 的值就是嵌入的隐秘信息。
EMD算法
EMD算法虽然可以高效地嵌入信息,简单地提取信息,但是对于像素灰度值的修改是永久性的。因此,提取出隐秘信息后,无法还原原始图像。而Chang在文献《Reversible Data Hiding Scheme Using Two Steganographic Images》中提出的算法很好的解决了这个问题,首次将EMD算法应用到可逆信息隐藏中。
emd算法二基本原理二信息隐藏方面的应用在文献efficientsteganographicembeddingbyexploitingmodificationdirection中zhangxinpeng首次提出了emd算法emd算法仅通过修改一组像素中一个像素的灰度值就可以嵌入一个2nl进制的隐秘信息具有图像像素值改变量小嵌入信息量大的特点
EMD算法
分解过程是：找出原数据序列X(t）所有的极大值点并用三次样条插值函数拟合形成原数据的上包络线；同样，找出所有的极小值点，并将所有的极小值点通过三次样条插值函数拟合形成数据的下包络线，上包络线和下包络线的均值记作ml，将原数据序列X(t）减去该平均包络 ml，得到一个新的数据序列h1：

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

%此版本为 ALAN 版本的整合注释版
function imf =emd(x
%Empiricial Mode Decomposition (Hilbert-HuangTransform
%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 =x-x1;
end
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, u =1; end
function s =getspline(x
%三次样条函数拟合成元数据包络线
N =length(x;
p =findpeaks(x;
s =spline([0p N+1],[0x(p0],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{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([kk k],'MarkerSize',3;hold on;
set(gca,'FontSize',8,'XLim',[0c(end],'YLim',[050];%设置 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',[0c(end];
title('EMD分解结果 ';
end
xlabel('Time';
end