基于二进制信源N次扩展熵率的MATLAB仿真.doc

基于二进制信源N次扩展熵率的MATLAB仿真

摘要：本文通过MATLAB分析一个二进制的二符号信源，分别采用了滑窗法、卷积法、变换形法求出了该信源的N次（其中N最大为8）扩展信源的熵率，得出了熵率随着N的增加而逐渐减小的结论。此外，本文还简单比较了三种算法的计算效率与优缺点。
关键词：MATLAB 二进制信源 熵率

一、问题的重述和分析
给定一个二进制信源，共包含n个样本，发出0、1的概率分别为p、1-p。对给定的二进制信源做N次扩展后，分析扩展后符号的概率统计特性，计算其熵率HN(XN)并分析H(XN)与N的关系。
本文中，令N=840，p=0.6，N最大取值为8。下面，通过三种办法计算该二进制信源N次扩展的熵率。
（一）滑窗法
该算法的基本思想是设定一个宽度为N的窗口，由符号序列的第1个符号滑动到倒数第N个符号，通过进制转换记录这几组数据中每个窗口消息组成的序列出现的频度，近似N次扩展信源每个消息发出的概率，最终求得其熵率HN(XN)并绘出HN(XN)-N图（注：计算熵率时可以过滤掉概率为零的分量，因为x→0时，x*lbx→0）。
MATLAB算法代码：
clc
sample=840;%给定独立二进制信源[0,1]的样本数
step=8; %设定最大扩展阶数
probability_0=0.6;%设定符号0的概率
H=zeros(1,step)建立熵率向量
x=randsrc(1,sample,[0,1;0.6,0.4]);
for N=1:step
n=zeros(1,2^N);
order=2.^([0:N-1]);
for M=0:(2^N-1)
for i=1:(sample+1-N)
if(order*x(i:N+i-1)'==M)
n(M+1)=n(M+1)+1;%记录信源发出的每种符号序列的个数
end
end
end
p=n/(sample+1-N);%N次扩展的概率向量
p=p(find(p~=0));%过滤掉概率为0的分量，其原理是x→0时，x*lbx→0
H(N)=-p*log2(p')/N;%输出每种扩展下的熵率
end
H%输出熵率向量
stem([1:step],H)
xlabel('扩展阶数N=1:8');
ylabel('熵率');
输出：
H =

0.9756 0.9749 0.9742 0.9715 0.9685 0.9654 0.9587 0.9438
HN(XN)-N图：

（二）卷积法
该算法的基本思想同滑窗法基本相同，不同的是通过符号序列与2的幂次向量做卷积实现窗口的滑动与进制的转换。记录每一个数据出现的频度，近似N次扩展信源每个消息的发出概率，最终求得其熵率HN(XN)并绘出HN(XN)-N图（注：计算熵率时可以过滤掉概率为零的分量，因为x→0时，x*lbx→0）。
MATLAB算法代码：
clc
sample=840;%给定独立二进制信源[0,1]的样本数
step=8; %设定最大扩展阶数
probability_0=0.6;%设定符号0的概率
x=randsrc(1,sample,[0,1;0.6,0.4]);
H=zeros(1,step);%建立熵率向量
n=zeros(1,2^step);%建立向量n记录每种符号序列出现的个数
for N=1:step
order=2.^([0:N-1]);
y=conv(x,order);
y1=y(step:(length(y)+1-step));%只截取卷积的有效部分
for i=1:2^step
n(i)=length(find(y1==i-1));%记录信源发出的每种符号序列的个数
end
p=n/length(y1);%N次扩展的概率向量
p=p(find(p~=0));%过滤掉概率为0的分量，其原理是x→0时，x*lbx→0
H(N)=-p*log2(p')/N;%输出每种扩展下的熵率
end
H%输出熵率向量
plot([1:step],H,'b-',[1:step],H,'bo')
title('卷积法');
xlabel('扩展阶数N=1:8');
ylabel('熵率');
输出：
H =

0.9714 0.9699 0.9689 0.9682 0.9665 0.9625 0.9545 0.9384
HN(XN)-N图：

（三）变换形法
该算法的基本思想是将原序列重组为一个N行、n/N列的一个新序列。将每一列数据看做一个窗口，仍然通过进制转换记录每个窗口消息组成的序列出现的频度，近似N次扩展信源每个消息的发出概率，最终求得其熵率HN(XN)并绘出HN(XN)-N图（注：计算熵率时可以过滤掉概率为零的分量，因为x→0时，x*lbx→0）。
MATLAB算法代码：
clc
sample=840;%给定独立二进制信源[0,1]的样本数
step=8; %设定最大扩展阶数
probability_0=0.6;%设定符号0的概率
x=randsrc(1,sample,[0,1;0.6,0.4]);
H=zeros(1,step);%建立熵率向量
for N=1:step
x1=reshape(x,N,[]);%经过reshape后的x
n=zeros(1,2^N);
order=2.^([0:N-1]);
for M=0:(2^N-1)
for i=1:(sample/N)
if(order*x1(:,i)==M)
n(M+1)=n(M+1)+1;%记录信源发出的每种符号序列的个数
end
end
end
p=n*N/sample;%N次扩展的概率向量
p=p(find(p~=0));%过滤掉概率为0的分量，其原理是x→0时，x*lbx→0
H(N)=-p*log2(p')/N;%输出每种扩展下的熵率
end
H%输出熵率向量
plot([1:step],H,'b-',[1:step],H,'bo')
title('变换形法');
xlabel('扩展阶数N=1:8');
ylabel('熵率');
输出：
H =

0.9636 0.9627 0.9568 0.9472 0.9364 0.9018 0.8423 0.7561
HN(XN)-N图：

二、三种算法效率的比较
在上文提到的三种算法中，每种算法都有其优点和局限性，下面就三种算法的特点进行一些较为简单粗略的比较。
对于第一种方法——滑窗法，相比变换形法，该方法能够记录较多的消息序列组数。故在样本总数一定时，该方法求得结果的可靠性更高。其不足是算法的时间复杂度较卷积法与变换形法较大，在样本数为840的条件下，若扩展阶数大于10，响应时间很长，甚至无法响应。
对于第二种方法——卷积法，其算法效率是三种方法中最高的。该算法兼具有多数据处理即高可靠性等特点，利用卷积运算降低了算法的时间复杂度。在样本数为840的条件下，可实现扩展阶数为16的熵率计算(如下图)。

对于第三种方法——变换形法，虽然该算法的数据处理量可以做到较大，但由于该算法只考察了在每种重组变换下每一列组成序列的统计特性，相比另两种算法，丢失的信息量较大，可靠性不高。由于其反映的信息量是最少的，所以计算出的熵率是三种方法中最小的，且下降速率是最大的。
三种方法熵率计算的比较图像如下图所示（其中蓝线为卷积法，黑线为滑窗法，红线为变换形法），由该图像进一步的说明了三种算法的计算效率由高到低分别为卷积法、滑窗法、变换形法。

三、参考资料
[1]陈运.信息论与编码.2版.北京：电子工业出版社，2011.
9