上一篇文章介绍了通过小波包分解，将数据按频率分成四个波段。其实（二）、（三）两篇都是对数据的预处理。
这篇文章主要讲一下互信息MI矩阵的建立。
用到的数据是在第二篇文章处理之后的数据。
下面先来介绍一下互信息的概念：
互信息（MI），是一种信息熵，是测量随机变量的信息量的常用方法之一，可以用来估计任意两个电极的时间序列之间的相互依赖关系，也常用来衡量电极间的信号同步性大小。假设一个离散的随机变量 X，其具有 N 个不同的状态，将这 N 个状态划分到 M个不同的区域空间，通过计算这 M 个区域中第 i 个区域的值的分布密度，得到相应的概率密度 pi ，则可得到 X 的信息熵为 HX:
同样地，通过计算两个时间序列的联合概率密度 pij，可以计算得到两个离散的时间序列
X，Y 之间的联合信息熵为 HXY：
则两个时间序列 X，Y 之间的互信息 MIXY 可以定义为如下：
即：
根据以上介绍以及公式，先给出大家一个代码
这个代码是计算两个列向量之间MI值的
u1 = rand(4,1); %随机生成一个列向量u1（4*1的列向量）
u2 = [2;32;6666;5]; %列向量u2
wind_size = size(u1,1); %返回u1的行数
n=wind_size; %n=4
x=[u1,u2];
[xrow, xcol] = size(x);
bin = zeros(xrow,xcol);
pmf = zeros(n, 2);
for i=1:2
minx=min(x(:,i));
maxx=max(x(:,i));
binwidth = (maxx - minx) / n;
edges = minx + binwidth*(0:n);
histcEdges = [-Inf edges(2:end-1) Inf];
[occur,bin(:,i)] = histc(x(:,i),histcEdges,1);
pmf(:,i) = occur(1:n)./xrow;
end
jointOccur = accumarray(bin,1,[n,n]); %（xi，yi）两个数据同时落入n*n等分方格中的数量即为联合概率密度
jointPmf = jointOccur./xrow;
Hx = -(pmf(:,1))'*log2(pmf(:,1)+eps);
Hy = -(pmf(:,2))'*log2(pmf(:,2)+eps);
Hxy = -(jointPmf(:))'*log2(jointPmf(:)+eps);
MI = Hx+Hy-Hxy;
mi = MI/sqrt(Hx*Hy);
最后的mi变量，即为两个列向量之间的MI值
下来分析我们的数据，对于文章（二）处理过的文件 ‘s01\s01-1.mat’来说，其data域的数据大小为32*2560（20s数据），为了提高实验的准确性，我们应该选取更小的时间间隔，这里选取10s的数据（其实应该选取更小一点的，但是方法相同）。
这里，只说一下思路，代码和文章（二）类似。
把1–10s的数据存到data1域，11–20s的数据存到data2域。即每部分的数据大小为32*1280。
处理好之后，我们正式开始建立MI矩阵。
思路如下：
对data1矩阵来说，32*1280，每次取出两行进行计算（取的方法有点类似于选择排序的两重循环），由于刚才给出的参考代码是两个列向量，所以，我们取出两行后，先对这两个行向量进行转置，成两个列向量后，再进行计算。
代码如下：
xx=load(E:\脑电数据集\再截取后数据\s01\s01-1.mat);
xx1=xx.data1;
xx2=xx.data2;
%处理前10秒的数据（data1），算MI1
n=1280;
MI1=zeros(32,32);
for p=1:32
for q=1:32
u1=xx1(p,:);
u2=xx1(q,:);
u1=u1';
u2=u2';
x=[u1,u2];
[xrow,xcol]=size(x);
bin=zeros(xrow,xcol);
pmf=zeros(n,2);
for i=1:2
minx=min(x(:,i));
maxx=max(x(:,i));
binwidth=(maxx-minx)/n;
edges=minx+binwidth*(0:n);%
histcEdges=[-Inf edges(2:end-1) Inf];
[occur,bin(:,i)]=histc(x(:,i),histcEdges,1);
pmf(:,i) = occur(1:n)./xrow;
end
jointOccur = accumarray(bin,1,[1280 1280]);
jointPmf = jointOccur./xrow;
Hx = -(pmf(:,1))'*log2(pmf(:,1)+eps);
Hy = -(pmf(:,2))'*log2(pmf(:,2)+eps);
Hxy = -(jointPmf(:))'*log2(jointPmf(:)+eps);
MIxy = Hx+Hy-Hxy;
mi = MIxy/sqrt(Hx*Hy);
MI1(p,q)=mi;
end
end
%处理后10秒的数据（data2），算MI2
n=1280;
MI2=zeros(32,32);
for p=1:32
for q=1:32
u1=xx2(p,:);
u2=xx2(q,:);
u1=u1';
u2=u2';
x=[u1,u2];
[xrow,xcol]=size(x);
bin=zeros(xrow,xcol);
pmf=zeros(n,2);
for i=1:2
minx=min(x(:,i));
maxx=max(x(:,i));
binwidth=(maxx-minx)/n;
edges=minx+binwidth*(0:n);%
histcEdges=[-Inf edges(2:end-1) Inf];
[occur,bin(:,i)]=histc(x(:,i),histcEdges,1);
pmf(:,i) = occur(1:n)./xrow;
end
jointOccur = accumarray(bin,1,[1280 1280]);
jointPmf = jointOccur./xrow;
Hx = -(pmf(:,1))'*log2(pmf(:,1)+eps);
Hy = -(pmf(:,2))'*log2(pmf(:,2)+eps);
Hxy = -(jointPmf(:))'*log2(jointPmf(:)+eps);
MIxy = Hx+Hy-Hxy;
mi = MIxy/sqrt(Hx*Hy);
MI2(p,q)=mi;
end
end
save('E:\脑电数据集\MI矩阵\s01\s01-1','MI1','MI2');
运行之后就可以得到32*32的MI矩阵啦
附一张运行之后的MI矩阵截图（只截一部分）
MI矩阵就建立好了
该文章是学习笔记，后续会继续更新