基于互信息的EEG脑网络情感识别（四）——建立MI矩阵

上一篇文章介绍了通过小波包分解，将数据按频率分成四个波段。其实（二）、（三）两篇都是对数据的预处理。

这篇文章主要讲一下互信息MI矩阵的建立。

用到的数据是在第三篇文章分完波段的数据，又将20s的数据分成了两个10s，然后进行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’来说，其每个波段的数据大小为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矩阵就建立好了

---

以上分析的运行结果是正确的，但是时间比较长

所以我对代码的一部分进行了修改

之前是通过两重循环（类似选择排序的思想）依次取出两行元素，现在用nchoosek(1 : phases, 2)，可以按照排列组合的形式取出元素（这里的phases是矩阵的行数）。
我把之前的代码改成了一个函数，直接调用函数就好。

    function MI = matrix(xx1)
    %计算MI矩阵的函数 
       phases = size(xx1, 1); 
       posiblecouple_phases = nchoosek(1 : phases, 2);
       MI11 = zeros(phases, phases);
       n=1280;
       for comb = 1 : size(posiblecouple_phases, 1)
           x = xx1(posiblecouple_phases(comb, :), :);
           x=x';
           [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);
           MIxy = Hx+Hy-Hxy;
           mi = MIxy/sqrt(Hx*Hy);
           r_c = posiblecouple_phases(comb, :);    %localizing row and column to storage R value
           MI11(r_c(1), r_c(2)) = mi;        
       end
       MI=MI11+MI11';
    end

关于代码，这里需要说明的是，MI11是一个上三角矩阵（不含对角线），所以要通过 MI=MI11+MI11’;将两个矩阵拼成一个MI矩阵，而对角线就为0.
下面附上矩阵的部分图：

---

该文章是学习笔记，后续会继续更新
下一篇文章的链接：基于互信息的EEG脑网络情感识别（五）