实验3_率失真函数的计算(1)

1、信息论与编码实验3率失真函数计算的程序设计学院：物理与电子学院班级：电信1105班姓名：学号：140411072一、实验问题假定一个DMS的信源符号集为 Au=1,2,r,其概率分布为 p(u);信宿符号集为 Av=1,2,s。而失真侧度矩阵为一个r s维矩阵D=dj。利用Matlab画出率失真函数 R()的曲线图。二、预备知识预备知识：信息论第二章的熵是针对不失真的情况，而在实际应用中只需要保留信息的主要特征即可，信号允许一定程度的失真，而率失真理论就是在这种情况下提出的。在实际问题中，信号有一定的失真是可以容忍的，但是当失真大于某个限度是，信号将会被严重损害，为此引入了失真函数d的概念。三

2、、实验目的1. 了解率失真函数性质、意义。2. 掌握简单的率失真函数计算方法；3. 掌握使用Matlab实现一般率失真函数的计算方法；4. 掌握Matlab求解非线性方程组的方法。四、实验要求1. 提前预习实验，认真阅读实验原理。2. 认真高效的完成实验，实验过程中服从实验室管理人员以 及实验指导老师的管理。五、实验内容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

3、给定 r=s=2。 p(u=1)=0.4,p=(u=2)=0.6;d=0,1;1,0 ，测试程序，即比较程序运行结果与理论计算结 果，R() H(P)H() 0 p0P6改变参数，画出函数图。7显示在计算精度为 0.000001 以及运行计算的配置 (CPU型号、 CPU 的频率、内存的 )的条件下，系统循环次数、累 计计算时间、平均每次循环所用时间等。六、实验原理1. R( (S)的表示方法计算 min 和 max 是 很 容 易 的 。 min p(u).mind(u,v) ； v u AUmaxmjnp(u)d(u,v)。当 max 时， R( )=0。u当 min max 时，R( )

4、=minl(U;V) : E(d)= o在数学上，就是在约束条件：rsE(d) p(u)p(v|u)d(u,v)u 1 v 1(1)sp(v|u) 1v1的约束下求平均无信息量|(U；V)p(u)p(v|u)logu 1 v 1勢(u，v)(9)的条件极小值。为此引入待定常数S和u(u=1,2,r),并作辅助函数Fp(v|u)p(u)p(v|u)logu 1 v 1p(v|u)p(v)r sS p(u)p(v|u)d(u,v)u 1 v 1u p(v|u)u 1 v 1其中p(v)p(u)p(v|u)u 1由喘荷0得，p(v|u)p(v)e初 Sd(u,v)up(u)为方便引入参数，u呗備则有

5、p(v|u) p(v) ueSd(u,v)显然提供r s个方程，提供r个方程，而(1)提供1各方 程，共rs+r+1个方程；而有r s个未知数p(v|u)、r个未知数u及 未知数S，共rs+r+1未知数，显然可以求解。为方便起见，我们 保留S作为参数。这样得到:Sd(u,v)p(v)eSd(u,v)p(u)eE(d) (S)p(u)p(v) ueSd(u，v)d(u,v)(7)u 1 v 1rR( (S) S (S) p(u)log uu 1(8)很容易得到S dR 0，即S是率是失真函数的导数。当S -d时，(S) min ；参量S是的递增函数，当从min到max逐渐增 大时，S将随 增大而

6、增大，当 =max时,S达到最大值Smax0。 对Smax的求解较麻烦，必须解非线性方程。为了简单我们不求Smax。如果r=s，即信源和信宿的符号集相同，则很容易通过(6)式求得u，进而通过 式求得p(v)。从而通过 (8)式划出率失真函 数曲线。2. R( (S)的迭代计算但一般情况下，r s,则只能通过 先求得p(v)，这是一个非 常复杂的方程。下面介绍 R( (S)的迭代方法计算方法和公式。首先假设p(v)固定，与信道传递概率p(v|u)无关，则求极值得:p*(v|u)Sd(u,v)p(v)esSd (u ,v)p(v)ev 1再假定p(v|u)不变，而把p(v)当成变量，则求极值得:r

7、p (v)p(u)p(v| u)u 1(10)具体算法为：选择绝对值相当大的负数Si。选定起始传递概率 p(v|u)=1/rs。通过(10)式求得P(v)，再通过(9)式求得p(v|u)。如此重复直到r sD(S)( n)p(u)p(n) (v|u)d(u,v)与 D(Si)( n+1)相差较小；并u 1 v 1且。r sp(n)( | )R(S)(n)p(u)pn)(v|u)log p 丄：1?与 R(S“(n+1)相差较小u 1 v 1p (v)再选择较大的S2直到Smax逼近于零为止。这样就可以画出R()曲线七、实验代码fun ctio nR,delta=R_delta(Pu,D)r,s

8、=size(D);eps=0.000001;delta_m in=0;for i=1:rdelta_min = delta_min + Pu(i)*m in( D(i,:);enddelta_max=zeros(1,s);for j=1:sfor i=1:rdelta_max(j) = delta_max(j) + Pu(i)*D(i,j);endend delta_max=min(delta_max);R=;delta=;P=ones(r,s);P=P/s;SS=100:-0.1:-100;SS=-exp(SS);for S=SSPv=Pu*P;Ed0=sum(Pu*(P.*D);Rs0=0

9、;for u=1:rfor v=1:sif P(u,v)=0 & Pu(u)=0Rs0 = Rs0 +Pu(u)*P(u,v)*log(P(u,v)/Pv(v);endendendP=exp(S*D);for i=1:sP(:,i)=P(:,i)*Pv(i);endfor i=1:rSumP=sum(P(i,:);P(i,:)=P(i,:)/SumP; endKm=50000;for k=1:KmPv=Pu*P;Edn=sum(Pu*(P.*D);Rsn=0;for u=1:rfor v=1:sif P(u,v)=0 & Pu(u)=0Rsn = Rsn +Pu(u)*P(u,v)*log(P

10、(u,v)/Pv(v);endendendP=exp(S*D);for i=1:sP(:,i)=P(:,i)*Pv(i);endfor i=1:rSumP=sum(P(i,:);P(i,:)=P(i,:)/SumP;end- Rs0)epsif abs(Edn - Ed0)eps & abs(Rsn break ;endEd0=Edn;Rs0=Rsn;endif kKm R=R,Rsn;delta=delta,Edn;endend function h=Hp(p) if p=0 & p=1h = -p*log(p)-(1-p)*log(1-p); elseh=0;end matlab 命令 clear all D=0,1;2,0;p=0.5;Pu=p,1-p; R,delta=R_delta(Pu,D);delta 仁0:p/1000