一、EMD简介
1 经验模态分解
EMD的本质是由数据的特征时间尺度来获得数量不同的本征模函数（intrinsic mode function,IMF），不同的本征模分量IMF代表不同的特征波动序列，使原始数据的波动特征在不同时间尺度下突显出来，由于5种环境时间序列具有一定的随机性和间断性，通过EMD分解，可在丰富输入变量多样性的同时，根据得到的IMF分量，突出环境序列在不同时间尺度下的局部特性，反映出原始环境序列的波动性、周期性和趋势变化，其具体分解过程如下：
1）对于一个原始数据序列x(t)，找到它所有的极大值点确定为上包络线，所有极小值点确定为下包络线，m(t)表示上包络线和下包络线的均值，第1个分量h1(t)=x(t)-m(t);
2）在第2个筛选过程中，h1(t)被视作原始数据，m1(t)是h1(t)的上下包络线的均值，如步骤1）所示确定第2个分量h2(t);
3）筛选过程重复n次，直到是hn(t)是一个本征模态函数或剩余分量rn(t)为一个单调函数，终止分解过程；
4）综上，指定q1(t)=h1(t),q2(t)=h2(t)，···，qk(t)=hn(t),x(t)最终分解为n个本征模分量qi(t)和一个剩余分量rn(t)，如式（2）所示。
2 主成分分析
经过EMD分解得到的数据序列充实了特征序列的数量，但是输入变量的维数也随之增多。为了在提高预测精度的同时，保持LSTM网络模型的计算速度，同时克服过拟合的问题，需通过PCA对输入变量进行降维处理，在保证信息有效性和代表性的前提下，提升模型的计算效率和精度。
PCA作为一种经典的数据降维方法，通过线性变换将原始数据转换到新的特征空间中，以此来提取数据的最主要线性分量，算法步骤如下:
1）得到原始数据的标准化矩阵；
2）根据式（1）计算标准化矩阵的相关系数矩阵R=(rij)m×n;
3）计算特征值矩阵和特征值对应的特征向量；
4）计算贡献率τi和累计贡献率ηi，计算方法如式（3）和式（4）所示，根据ηi确定需选取的主成分的个数；
5）通过选取的主成分的特征值和所对应的特征向量，最终可得到降维后的数据。
二、 LSTM简介
1 LSTM控制流程
LSTM的控制流程：是在前向传播的过程中处理流经细胞的数据，不同之处在于 LSTM 中细胞的结构和运算有所变化。
这一系列运算操作使得 LSTM具有能选择保存信息或遗忘信息的功能。咋一看这些运算操作时可能有点复杂，但没关系下面将带你一步步了解这些运算操作。
2 核心概念
LSTM 的核心概念在于细胞状态以及“门”结构。细胞状态相当于信息传输的路径，让信息能在序列连中传递下去。你可以将其看作网络的“记忆”。理论上讲，细胞状态能够将序列处理过程中的相关信息一直传递下去。
因此，即使是较早时间步长的信息也能携带到较后时间步长的细胞中来，这克服了短时记忆的影响。信息的添加和移除我们通过“门”结构来实现，“门”结构在训练过程中会去学习该保存或遗忘哪些信息。
3 Sigmoid
门结构中包含着 sigmoid 激活函数。Sigmoid 激活函数与 tanh 函数类似，不同之处在于 sigmoid 是把值压缩到 0~1 之间而不是 -1~1 之间。这样的设置有助于更新或忘记信息，因为任何数乘以 0 都得 0，这部分信息就会剔除掉。同样的，任何数乘以 1 都得到它本身，这部分信息就会完美地保存下来。这样网络就能了解哪些数据是需要遗忘，哪些数据是需要保存。
4LSTM门结构
LSTM 有三种类型的门结构：遗忘门、输入门和输出门。
4.1 遗忘门
遗忘门的功能是决定应丢弃或保留哪些信息。来自前一个隐藏状态的信息和当前输入的信息同时传递到 sigmoid 函数中去，输出值介于 0 和 1 之间，越接近 0 意味着越应该丢弃，越接近 1 意味着越应该保留。
4.2 输入门
输入门用于更新细胞状态。首先将前一层隐藏状态的信息和当前输入的信息传递到 sigmoid 函数中去。将值调整到 0~1 之间来决定要更新哪些信息。0 表示不重要，1 表示重要。
其次还要将前一层隐藏状态的信息和当前输入的信息传递到 tanh 函数中去，创造一个新的侯选值向量。最后将 sigmoid 的输出值与 tanh 的输出值相乘，sigmoid 的输出值将决定 tanh 的输出值中哪些信息是重要且需要保留下来的。
4.3 细胞状态
下一步，就是计算细胞状态。首先前一层的细胞状态与遗忘向量逐点相乘。如果它乘以接近 0 的值，意味着在新的细胞状态中，这些信息是需要丢弃掉的。然后再将该值与输入门的输出值逐点相加，将神经网络发现的新信息更新到细胞状态中去。至此，就得到了更新后的细胞状态。
4.4 输出门
输出门用来确定下一个隐藏状态的值，隐藏状态包含了先前输入的信息。首先，我们将前一个隐藏状态和当前输入传递到 sigmoid 函数中，然后将新得到的细胞状态传递给 tanh 函数。
最后将 tanh 的输出与 sigmoid 的输出相乘，以确定隐藏状态应携带的信息。再将隐藏状态作为当前细胞的输出，把新的细胞状态和新的隐藏状态传递到下一个时间步长中去。
让我们再梳理一下。遗忘门确定前一个步长中哪些相关的信息需要被保留；输入门确定当前输入中哪些信息是重要的，需要被添加的；输出门确定下一个隐藏状态应该是什么。
三、部分源代码
clear;
close all;
clc;
%% 参数设置
%% 数据加载
f=xlsread('215-风速数据.xls','B2706:B3449');
u=emd(f);%emd
% 如果用eemd ceemd 则还有参数需要设置
Nstd = 0.02;
NR = 1;
MaxIter = 200;
% [u ,~]=eemd(f,Nstd,NR,MaxIter);%eemd
% [u ,~]=ceemd(f,Nstd,NR,MaxIter);%ceemd
K=size(u,1);
figure
subplot(K+1,1,1)
plot(f)
ylabel('原始')
for i=2:K+1
subplot(K+1,1,i)
plot(u(i-1,:))
ylabel(['IMF',num2str(i-1)])
end
ylabel('res')
save emd_data u
function [modes its]=ceemd(x,Nstd,NR,MaxIter)
% WARNING: for this code works it is necessary to include in the same
%directoy the file emd.m developed by Rilling and Flandrin.
%This file is available at %http://perso.ens-lyon.fr/patrick.flandrin/emd.html
%We use the default stopping criterion.
%We use the last modification: 3.2007
%
% This version was run on Matlab 7.10.0 (R2010a)
%----------------------------------------------------------------------
% INPUTs
% x: signal to decompose
% Nstd: noise standard deviation
% NR: number of realizations
% MaxIter: maximum number of sifting iterations allowed.
%
% OUTPUTs
% modes: contain the obtained modes in a matrix with the rows being the modes
% its: contain the sifting iterations needed for each mode for each realization (one row for each realization)
% -------------------------------------------------------------------------
% Syntax
%
% modes=ceemdan(x,Nstd,NR,MaxIter)
% [modes its]=ceemdan(x,Nstd,NR,MaxIter)
%
%--------------------------------------------------------------------------
% This algorithm was presented at ICASSP 2011, Prague, Czech Republic
% Plese, if you use this code in your work, please cite the paper where the
% algorithm was first presented.
% If you use this code, please cite:
%
% M.E.TORRES, M.A. COLOMINAS, G. SCHLOTTHAUER, P. FLANDRIN,
% "A complete Ensemble Empirical Mode decomposition with adaptive noise,"
% IEEE Int. Conf. on Acoust., Speech and Signal Proc. ICASSP-11, pp. 4144-4147, Prague (CZ)
%
% -------------------------------------------------------------------------
% Date: June 06,2011
% Authors: Torres ME, Colominas MA, Schlotthauer G, Flandrin P.
% For problems with the code, please contact the authors:
% To: macolominas(AT)bioingenieria.edu.ar
% CC: metorres(AT)santafe-conicet.gov.ar
% -------------------------------------------------------------------------
x=x(:)';
desvio_x=std(x);
x=x/desvio_x;
modes=zeros(size(x));
temp=zeros(size(x));
aux=zeros(size(x));
acum=zeros(size(x));
iter=zeros(NR,round(log2(length(x))+5));
for i=1:NR
white_noise{i}=randn(size(x));%creates the noise realizations
end;
for i=1:NR
modes_white_noise{i}=emd(white_noise{i});%calculates the modes of white gaussian noise
end;
for i=1:NR %calculates the first mode
temp=x+Nstd*white_noise{i};
[temp, o, it]=emd(temp,'MAXMODES',1,'MAXITERATIONS',MaxIter);
temp=temp(1,:);
aux=aux+temp/NR;
iter(i,1)=it;
end;
modes=aux; %saves the first mode
k=1;
aux=zeros(size(x));
acum=sum(modes,1);
while nnz(diff(sign(diff(x-acum))))>2 %calculates the rest of the modes
for i=1:NR
tamanio=size(modes_white_noise{i});
if tamanio(1)>=k+1
noise=modes_white_noise{i}(k,:);
noise=noise/std(noise);
noise=Nstd*noise;
try
[temp, o, it]=emd(x-acum+std(x-acum)*noise,'MAXMODES',1,'MAXITERATIONS',MaxIter);
temp=temp(1,:);
catch
it=0;
temp=x-acum;
end;
else
[temp, o, it]=emd(x-acum,'MAXMODES',1,'MAXITERATIONS',MaxIter);
temp=temp(1,:);
end;
aux=aux+temp/NR;
iter(i,k+1)=it;
end;
modes=[modes;aux];
aux=zeros(size(x));
acum=zeros(size(x));
acum=sum(modes,1);
k=k+1;
end;
modes=[modes;(x-acum)];
[a b]=size(modes);
iter=iter(:,1:a);
modes=modes*desvio_x;
its=iter; 四、运行结果
五、matlab版本及参考文献
1 matlab版本
2014a
2 参考文献
[1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例（第2版）[M].电子工业出版社，2016.
[2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社，2017.
[3]周品.MATLAB 神经网络设计与应用[M].清华大学出版社，2013.
[4]陈明.MATLAB神经网络原理与实例精解[M].清华大学出版社，2013.
[5]方清城.MATLAB R2016a神经网络设计与应用28个案例分析[M].清华大学出版社，2018.