mi.m

多变量的互信息的matlab的函数的源码

clear all
close all
clc
L=2^10;m=10;
f=2;fs=50;
for n=1:10000;
y(n)=sin(2*pi*f*n/fs);
u(n)=0.02*sin(2*pi*f*n/fs);
end
s=y+u;d=y+m*u;
S=paixu(s);D=paixu(d);    %对s和d进行从小到大排序
z1=(S(n)-S(1))/L;    %把S和D的值域平均分成L个区域，每个区域的宽度
z2=(D(n)-D(1))/L;     
%统计落在区域(i,j)内的点数分别为N(i,j)
N=zeros(L,L);
for k=1:n
    if s(k)==S(1)
        for r=1:L
            if d(k)==D(1)
                N(1,1)=N(1,1)+1;
            elseif D(1)+(r-1)*z2<d(k)&&d(k)<=D(1)+r*z2  
                N(1,r)=N(1,r)+1;
            end
        end
    elseif d(k)==D(1)
        for t=1:L
            if S(1)+(t-1)*z1<s(k)&&s(k)<=S(1)+t*z1
               N(t,1)=N(t,1)+1;
            end
        end
    end
end
for w=1:n
    for i=1:L
        if S(1)+(i-1)*z1<s(w)&&s(w)<=S(1)+i*z1
             for j=1:L
                 if D(1)+(j-1)*z2<d(w)&&d(w)<=D(1)+j*z2
                     N(i,j)= N(i,j)+1;
                 end
             end
        end
    end
end
%计算联合概率
p=zeros(L,L);
for i=1:L
    for j=1:L
        p(i,j)=N(i,j)/n;
    end 
end
%计算互信息
Is=entropy(s);Id=entropy(d);
disp(Is);
disp(Id);
Isd=0;
for i=1:L
    for j=1:L
        if p(i,j)==0
            continue;
        else
            Isd=Isd-p(i,j)*log2(p(i,j));
        end
    end
end
disp(Isd);
MI=Is+Id-Isd;
disp(MI);
NMI=(Is+Id)/Isd;
disp(NMI);