[zz]求一維序列的信息熵（香濃熵）的matlab程序實例

對於一個二維信號，好比灰度圖像，灰度值的範圍是0-255，所以只要根據像素灰度值(0-255)出現的機率，就能夠計算出信息熵。
    可是，對於一個一維信號，好比說心電信號，數據值的範圍並非肯定的，不會是(0-255)這麼肯定，若是進行域值變換，使其轉換到一個整數範圍的話，就會丟失數據，請高手指點，怎麼計算。

好比數字信號是x(n),n=1~N
(1)先用Hist函數對x(n)的賦值範圍進行分塊，好比賦值範圍在0~10的對應第  
    一塊，10~20的第二塊，以此類推。這以前須要對x(n)作一些歸一化處理
(2)統計每一塊的數據個數，並求出相應的機率
(3)用信息熵公式求解

以上求解方法得到的雖然是近似的信息熵，可是通常認爲，這麼作是沒有問題的

求一維序列的信息熵的matlab程序代碼以下：（已寫成調用的函數形式）

測試程序：

fs=12000;
N=12000;
T=1/fs;
t=(0:N-1)*T;
ff=104;
sig=0.5*(1+sin(2*pi*ff*t)).*sin(2*pi*3000*t)+rand(1,length(t));

Hx=yyshang(sig,10)

 

%———————求一維離散序列信息熵matlab代碼

function Hx=yyshang(y,duan)
%不以原信號爲參考的時間域的信號熵
%輸入：maxf:原信號的能量譜中能量最大的點
%y:待求信息熵的序列
%duan:待求信息熵的序列要被分塊的塊數
%Hx:y的信息熵
%duan=10;%將序列按duan數等分，若是duan=10,就將序列分爲10等份
x_min=min(y);
x_max=max(y);
maxf(1)=abs(x_max-x_min);
maxf(2)=x_min;
duan_t=1.0/duan;
jiange=maxf(1)*duan_t;
% for i=1:10
% pnum(i)=length(find((y_p>=(i-1)*jiange)&(y_p<i*jiange)));
% end

pnum(1)=length(find(y<maxf(2)+jiange));
for i=2:duan-1
    pnum(i)=length(find((y>=maxf(2)+(i-1)*jiange)&(y<maxf(2)+i*jiange)));
end
pnum(duan)=length(find(y>=maxf(2)+(duan-1)*jiange));
%sum(pnum)
ppnum=pnum/sum(pnum);%每段出現的機率
%sum(ppnum)
Hx=0;
for i=1:duan
    if ppnum(i)==0
        Hi=0;
    else
        Hi=-ppnum(i)*log2(ppnum(i));
    end
    Hx=Hx+Hi;
end
end

%----------------

擴展閱讀：

 

 

實驗一：計算離散信源的熵

1、實驗設備:

一、計算機

二、軟件：Matlab

2、實驗目的:

  一、熟悉離散信源的特色；

  二、學習仿真離散信源的方法

  三、學習離散信源平均信息量的計算方法

  四、熟悉 Matlab 編程；

3、實驗內容:

   一、寫出計算自信息量的Matlab 程序

   二、寫出計算離散信源平均信息量的Matlab 程序。

   三、掌握二元離散信源的最大信息量與機率的關係。

   四、將程序在計算機上仿真實現，驗證程序的正確性並完成習題。

4、實驗報告要求

簡要總結離散信源的特色及離散信源平均信息量的計算,寫出習題的MATLAB實現語句。

 

信息論基礎：

自信息的計算公式

              Matlab實現：I=log2(1/p) 或I=-log2(p)

熵（平均自信息）的計算公式

  

Matlab實現：HX=sum(-x.*log2(x))；或者h=h-x(i)*log2(x(i));

習題：

1. 甲地天氣預報構成的信源空間爲：

 

乙地信源空間爲：

 

求此兩個信源的熵。求各類天氣的自信息量。

案：

運行程序：

p1=[1/2,1/4,1/8,1/8];%p1表明甲信源對應的機率

p2=[7/8,1/8];%p2表明乙信源對應的機率

H1=0.0;

H2=0.0;

I=[];

J=[];

for i=1:4

    H1=H1+p1(i)*log2(1/p1(i));

    I(i)=log2(1/p1(i));

end

disp('自信息量分別爲：');

I

disp('H1信源熵爲：');

H1

for j=1:2

    H2=H2+p2(j)*log2(1/p2(j));

     J(j)=log2(1/p2(j));

end

disp('自信息量分別爲：');

J

disp('H2信源熵爲：');

H2