MATLAB实现PCA去除眼电信号

一、去除伪迹的讨论

查阅文献[1]可知，所获取的低频脑电信号中，主要受心电（ECG）和眼电（EOG）的干扰较大。

但一般实验会忽略掉心电伪迹（为了简化，俺也是），若非要去除，可采用：在事件相关电位（ERP）研究中，通过设置刺激信号的随机出现，再对脑电信号叠加平均即可去除心电干扰。另外也可以通过成分分离的方法去除。
眼电的具体去除见后文，此处仅科普产生机理…
→ 眼电的产生机理：由于跨膜间存在静息电位，角膜侧为正，巩膜侧为负，这可以看作是在那里有一定大小的电偶极子存在。眼球转动过程中，偶极子的方向相对于记录电极发生改变。眼球初始位置与终止位置时的偶极子所产生的电位差，构成了眼电伪迹，眼球运动方向的电极记录到正的眼电伪迹，眼球运动相反方向的电极记录得到的眼电伪迹为负

除了上述两种信号伪迹，还有：
肌电信号EMG（解决办法：受试者尽量减少头颈部、手臂运动。加上肌电信号主要频率在100Hz以上，大部分已被之前的带通滤波去除）
呼吸（导致电极相对于头皮出现滑动）、皮肤电（被试出汗或者心理紧张所引起皮肤阻抗发生变化）
这类信号的频率较低，一般小于1Hz，脑电采集过程中将直接予以滤除。因此伪迹往往被忽略，不予考虑
非生理信号伪迹（在这些噪声中，50Hz的工频干扰尤为严重, 频谱不明显，脑电采集放大器已滤除大部分）

二、眼电去除

前处理：利用之前获取的低频脑电信号，plot作图发现由于信号采集过程中电极与皮肤产生相对位移，信号受干扰较大，故只截取有效部分。

按上述，应该不出现此种干扰信号

问题的产生：上述伪迹该忽略的被忽略，该滤除的被滤除。处理后的信号理应还含有周期性的高幅值眼电信号。
去除眼电的大概思路为：用信号分解方法，将采集到的数据进行成分分解。然后将分离出的眼电成分从采集到的数据中减掉，从而得到没有眼电干扰的脑电数据。此类方法中比较常见的是主成分分析(PCA，Principal Component Analysis)、独立成分分析 (ICA，Independent ComponentAnalysis)等方法。

主成分分析算法去除眼电伪迹的基本理论是：

假设脑电信号与眼电信号彼此正交，将原始n维（脑电信号的n为512）特征数据映射到彼此正交的k维特征空间↓，从而实现原始信号的特征分离。

第一个新坐标轴选择是原始数据中方差最大的方向，第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的，第三个轴是与第1,2个轴正交的平面中方差最大的。依次类推，可以得到n个这样的坐标轴。通过这种方式获得的新的坐标轴，我们发现，大部分方差都包含在前面k个坐标轴中，后面的坐标轴所含的方差几乎为0。于是，我们可以忽略余下的坐标轴，只保留前面k个含有绝大部分方差的坐标轴。事实上，这相当于只保留包含绝大部分方差的维度特征，而忽略包含方差几乎为0的特征维度，实现对数据特征的降维处理。

为什么一秒512个小包的数据是特征

识别并去除眼电伪迹后（PCA应用于眼电去噪时，由于通常来说眼电的能量远远大于脑电的能量，表现为眼电噪声的成分为第一个主成分，故通过舍弃第一主成分就能够达到去噪的效果。），逆向投影恢复原始数据。

数学实现（含代码）：
1）协方差生成 _____ X[ 512,512 ] , S0[ 512,512 ]

事先说明：没有必要对特征进行归一化。因为信号数据b的来源相同，量纲及量纲单位相同，已具备可比性，无需采用标准化方法消除由此带来的偏差。
若针对高维不同类型特征的数据，应标准化，具体可见关于归一化、标准化、中心化的解释

X=b;                       %脑电信号
avg=mean(X,2);             %返回各行平均值
X=bsxfun(@minus,X,avg);    %bsxfun（fun,A,B）对两个矩阵A和B之间的每一个元素进行指定的计算（函数fun指定）；并且具有自动扩维的作用%每信号的512个维度减去各自均值
[~,N]=size(X);
S0=1/N*(X'*X);             %X的协方差矩阵

b为 512维特征*86例样本的脑电信号数据 ↗

2)用特征值分解方法求协方差矩阵S的特征值与特征向量_______U0[ 512,512 ]

[U0,aaa]=eig(S0);          %计算矩阵S0的特征值aaa和S0的特征列向量U0
U0=fliplr(U0);             %U0原本是按特征值升序排列的，要调换顺序（对U0的列倒序排列:123→321

要去除各维度之间的相关性，相当于使协方差矩阵S0的非对角元素全变为0，使S0对角化为D。因此我们需要找到一个矩阵U，使S对角化：↗

正好，S的特征向量组成的矩阵能达到这一目的，因此把S的特征向量放入U的每一列中即可。U的列向量即为主成分向量。D为S的特征值组成的对角阵。

3）降维、去除眼电、信号复原

U0(:,1:5)=0;               %滤除前五个主成分（数字可改）
U=(U0(:,1:20)')*X;         % U[20,512]=U0[20,512]*X[512,512]   用前k（k=20）个主成分，通过降维矩阵将样本映射到低维空间上
U=U0(:,1:20)*U;            % U[512,512]=U0[512,20]*U[20,512]  还原数据，获得在原空间中的估计位置
U=bsxfun(@plus,U,avg);     % 最后加上均值

4）效果图

用前20个特征还原出的脑电信号↓

去除20个特征中的第一个主成分↓

去除20个特征中的前五个主成分↓

三、优化

贡献率（降维的k的值的选择）
前面所用特征值之和 / 总特征值之和 = 前面几个主成分的累计贡献率

sum = sum(diag(aaa));              %总特征值之和sum1 = 0;
for i = (512 - 20 + 1):512         sum1=sum1+aaa(i,i);
end                                %所用前20个特征值之和（倒序）contribution = sum1/sum            % ans = 0.95

线代知识—特征值、特征向量

感谢阅读~

[1]刘铁军.脑电信号中眼电伪迹去除方法研究[D].四川:电子科技大学,2008.