实验3-率失真函数计算的程序设计

实验3 率失真函数计算的程序设计

一、实验问题

假定一个离散无记忆信源（DMS ）的信源符号集为A u ={1,2,...,r},其概率分布为p(u)；信宿符号集为A v ={1,2,....,s}。而失真侧度矩阵为一个r ⨯s 维矩阵D=[d ij ]。利用Matlab 画出率失真函数R(δ)的曲线图。

二、实验环境

计算机、Windows 2000 或以上、Matlab 6.5或以上

三、实验目的

1．了解率失真函数性质、意义。

2．掌握简单的率失真函数计算方法；

3．掌握使用Matlab 实现一般率失真函数的计算方法；

4．掌握Matlab 求解非线性方程组的方法。

四、实验内容

1．从理论上计算r=s=2。p(u=1)=p,p (u=2)=1-p;d=[0,1;1,0]的率失真函数R(δ)。

2．对一般性的DMS 信源，计算率失真函数R(δ)的理论公式进行推导。

3．找出比较合适的方程求解方法。

4．使用编制Matlab 编制程序求解一般的率失真函数R(δ)。

5．给定r=s=2。p(u=1)=0.4,p=(u=2)=0.6;d=[0,1;1,0]，测试程序，即比较程序运行结果与理论计算结果,

⎩⎨⎧≥≤≤-=p p H p H R δδδδ

00)()()( 6．改变参数，画出函数图。

7．显示在计算精度为0.000001以及运行计算的配置(CPU 型号、CPU 的频率、内存的)的条件下，系

统循环次数、累计计算时间、平均每次循环所用时间等。

五、实验要求

1. 提前预习实验，认真阅读实验原理。

2. 认真高效的完成实验，实验过程中服从实验室管理人员以及实验指导老师的管理。

3. 将实验报告写成论文的形式。要求有：

● 问题的提出：包括R(δ)的物理意义、用途（可以举出具体的用途）、计算的困难性等。

● 解决问题的原理方法：包括所有的公式推导的细节。

● 解决问题的具体方法：包括程序框图及Matlab 源程序。

● 实验结果：利用你的程序给出不同参数得到的实验结果，包括实验曲线图、程序循环次数、

累计计算时间、平均每次循环所用时间等。

● 结果分析：包括R(δ)的性质、程序收敛情况、程序改进的方向等。

4. 每个同学必须独立完成实验（不能抄袭，否则两人均为零分），实验成绩是该门课程成绩的主要

依据。

六、实验原理

1．R(δ(S))的表示方法

计算δmin 和δmax 是很容易的。

∑∈=

U A u v v u d u p ),(min ).(min δ； ∑=u

v v u d u p ),()(min max δ。当δ≥δmax 时，R(δ)=0。 当δmin ≤δ≤δmax 时，R(δ)=min{I(U;V)：E(d)= δ}。

在数学上，就是在约束条件：

δ==∑∑==),()|()()(11v u d u v p u p d E r u s

v (1)

1)|(1=∑=s

v u v p (2)

的约束下求平均信息量

∑∑===r u s v v u d v p u v p u v p u p V U I 11),()

()|(log

)|()();( 的条件极小值。为此引入待定常数S 和μu (u=1,2,...,r)，并作辅助函数

∑∑∑∑∑∑======+-=r u s

v u r u s v r u s v u v p v u d u v p u p S v p u v p u v p u p u v p F 1

11111

)|(),()|()()()|(log

)|()()]|([μ (3) 其中∑==r

u u v p u p v p 1)|()()( 由0)|(=∂∂u v p F 得，})

(),(ex p{)()|(u p v u Sd v p u v p u μ+= 为方便引入参数，})(ex p{u p u

u μλ=则有

),()()|(v u Sd u e v p u v p λ= (4)

显然(4)提供r ⨯s 个方程，(2)提供r 个方程，而(1)提供1各方程，共rs+r+1个方程；而有r ⨯s 个未知数p(v|u)、r 个未知数λu 及未知数S ，共rs+r+1未知数，显然可以求解。为方便起见，我们保留S 作为参数。这样得到：

1)(1),(=∑=s

v v u Sd u e v p λ (5)

1)(1),(=∑=r u v u Sd u e u p λ

(6)

∑∑====r u s

v v u Sd u v u d e v p u p S d E 11

),(),()()()()(λδ (7)

∑=+=r

u u u p S S S R 1log )()())((λδδ (8) 很容易得到0<=δ

d dR S ，即S 是率是失真函数的导数。当S →-∞时，δ(S)→ δmin ；参量S 是δ的递增函数，当δ从δmin 到δmax 逐渐增大时，S 将随δ增大而增大，当δ=δmax 时,S 达到最大值S max <0。对S max 的求解较

麻烦，必须解非线性方程。为了简单我们不求S max 。

如果r=s ，即信源和信宿的符号集相同，则很容易通过(6)式求得λu ，进而通过(5)式求得p(v)。从而通过(7)(8)式划出率失真函数曲线。

2．R(δ(S))的迭代计算

但一般情况下，r ≠s ，则只能通过(6)先求得p(v)，这是一个非常复杂的方程。下面介绍R(δ(S))的迭代方法计算方法和公式。

首先假设p(v)固定，与信道传递概率p(v|u)无关，则求极值得： ∑==s v v u Sd v u Sd e

v p e v p u v p 1)

,()

,(*)()()|( (9)

再假定p(v|u)不变，而把p(v)当成变量，则求极值得：

∑==r u u v p u p v p 1*)|()()( (10) 具体算法为：

选择绝对值相当大的负数S 1。选定起始传递概率p (1)(v|u)=1/rs 。

通过(10)式求得P (1)(v)，再通过(9)式求得p (2)(v|u)。如此重复直到

∑∑===r u s

v n v u d u v p u p n S D 11)(1),()|()())((与D(S 1)(n+1)相差较小；并且。

∑∑===r u s v n n n v p u v p u v p u p n S R 11)()()

(1)()|(log )|()())((与R(S 1)(n+1)相差较小 再选择较大的S 2直到S max 逼近于零为止。这样就可以画出R(δ)曲线.

七。实验步骤：

1.建立熵函数求解matlab 文件 Hp.m

function h=Hp(p)

%熵函数计算，输入是概率p （标量或一维矢量）

%输出是h,与p 同维度

h=zeros(1,length(p))

index=find(p>0&p<1);

P=p(index);

h(index) = -P.*log(P)-(1-P).*log(1-P);

end

2.根据实验原理，很容易知道，delta 的最小值和最大植：

delta_min=sum(Pu.*min(D',[],1));delta_max=min(Pu*D);

建立求率失真度的文件：R_delta.m

%率失真函数计算

function [R,delta]=R_delta(Pu,D,eps)

% Pu 信源概率矢量，D 失真测度矩阵（rxs 阶）

% 信源Ur ，信宿Vs