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

实验3率失真函数计算的程序设计一、实验问题假定一个离散无记忆信源（DMS）的信源符号集为Au={1,2,...,r},其概率分布为p(u)；信宿符号集为Av={1,2,....,s}。而失真侧度矩阵为一个rs维矩阵D=[dij]。利用Matlab画出率失真函数R()的曲线图。二、实验环境计算机、Windows2000或以上、Matlab6.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]，测试程序，即比较程序运行结果与理论计算结果,ppHpHR　　　　　　　　　00)()()(6．改变参数，画出函数图。7．显示在计算精度为0.000001以及运行计算的配置(CPU型号、CPU的频率、内存的)的条件下，系统循环次数、累计计算时间、平均每次循环所用时间等。五、实验要求1.提前预习实验，认真阅读实验原理。2.认真高效的完成实验，实验过程中服从实验室管理人员以及实验指导老师的管理。3.将实验报告写成论文的形式。要求有：问题的提出：包括R()的物理意义、用途（可以举出具体的用途）、计算的困难性等。解决问题的原理方法：包括所有的公式推导的细节。解决问题的具体方法：包括程序框图及Matlab源程序。实验结果：利用你的程序给出不同参数得到的实验结果，包括实验曲线图、程序循环次数、累计计算时间、平均每次循环所用时间等。结果分析：包括R()的性质、程序收敛情况、程序改进的方向等。4.每个同学必须独立完成实验（不能抄袭，否则两人均为零分），实验成绩是该门课程成绩的主要依据。六、实验原理1．R((S))的表示方法计算min和max是很容易的。UAuvvudup),(min).(min；uvvudup),()(minmax。当max时，R()=0。当minmax时，R()=min{I(U;V)：E(d)=}。在数学上，就是在约束条件：),()|()()(11vuduvpupdErusv(1)1)|(1svuvp(2)的约束下求平均信息量rusvvudvpuvpuvpupVUI11),()()|(log)|()();(的条件极小值。为此引入待定常数S和u(u=1,2,...,r)，并作辅助函数rusvurusvrusvuvpvuduvpupSvpuvpuvpupuvpF111111)|(),()|()()()|(log)|()()]|([(3)其中ruuvpupvp1)|()()(由0)|(uvpF得，})(),(exp{)()|(upvuSdvpuvpu为方便引入参数，})(exp{upuu则有),()()|(vuSduevpuvp(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),(svvuSduevp(5)1)(1),(ruvuSdueup(6)rusvvuSduvudevpupSdE11),(),()()()()((7)ruuupSSSR1log)()())(((8)很容易得到0ddRS，即S是率是失真函数的导数。当S-时，(S)min；参量S是的递增函数，当从min到max逐渐增大时，S将随增大而增大，当=max时,S达到最大值Smax0。对Smax的求解较麻烦，必须解非线性方程。为了简单我们不求Smax。如果r=s，即信源和信宿的符号集相同，则很容易通过(6)式求得u，进而通过(5)式求得p(v)。从而通过(7)(8)式划出率失真函数曲线。2．R((S))的迭代计算但一般情况下，rs，则只能通过(6)先求得p(v)，这是一个非常复杂的方程。下面介绍R((S))的迭代方法计算方法和公式。首先假设p(v)固定，与信道传递概率p(v|u)无关，则求极值得：svvuSdvuSdevpevpuvp1),(),(*)()()|((9)再假定p(v|u)不变，而把p(v)当成变量，则求极值得：ruuvpupvp1*)|()()((10)具体算法为：选择绝对值相当大的负数S1。选定起始传递概率p(1)(v|u)=1/rs。通过(10)式求得P(1)(v)，再通过(9)式求得p(2)(v|u)。如此重复直到rusvnvuduvpupnSD11)(1),()|()())((与D(S1)(n+1)相差较小；并且。rusvnnnvpuvpuvpupnSR11)()()(1)()|(log)|()())((与R(S1)(n+1)相差较小再选择较大的S2直到Smax逼近于零为止。这样就可以画出R()曲线.七。实验步骤：1.建立熵函数求解matlab文件Hp.mfunctionh=Hp(p)%熵函数计算，输入是概率p（标量或一维矢量）%输出是h,与p同维度h=zeros(1,length(p))index=find(p0&p1);P=p(index);h(index)=-P.*log(P)-(1-P).*log(1-P);end2.根据实验原理，很容易知道，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[r,s]=size(D);%delta的最小值和最大植delta_min=sum(Pu.*min(D',[],1));delta_max=min(Pu*D);R=[];delta=[];%给定初始值%P(u,v)=P(V|U)P=ones(r,s)/s;%信道模型SS=100:-0.1:-100;SS=-exp(SS);forS=SSPv=Pu*P;Ed0=sum(Pu*(P.*D));%Rs0=sum(Pu*(P.*log(P)))-sum(Pv.*log(Pv));Rs0=0;foru=1:rforv=1:sifP(u,v)~=0&Pu(u)~=0Rs0=Rs0+Pu(u)*P(u,v)*log(P(u,v)/Pv(v));endendendP=exp(S*D);fori=1:sP(:,i)=P(:,i)*Pv(i);endfori=1:rSumP=sum(P(i,:));P(i,:)=P(i,:)/SumP;endKm=50000;fork=1:KmPv=Pu*P;Edn=sum(Pu*(P.*D));%Rsn=sum(Pu*(P.*log(P)))-sum(Pv.*log(Pv));Rsn=0;foru=1:rforv=1:sifP(u,v)~=0&Pu(u)~=0Rsn=Rsn+Pu(u)*P(u,v)*log(P(u,v)/Pv(v));endendendP=exp(S*D);fori=1:sP(:,i)=P(:,i)*Pv(i);endfori=1:rSumP=sum(P(i,:));P(i,:)=P(i,:)/SumP;endifabs(Edn-Ed0)eps&abs(Rsn-Rs0)epsbreak;endEd0=Edn;Rs0=Rsn;endifkKmR=[R,Rsn];delta=[delta,Edn];endend3.测试，以汉明失真测度矩阵为失真测度，分别绘制不同精度和不同输入概率下的率失真函数曲线测试1：对于不同的精度测试程序如下clearallalpha=2;p=0.5;D=[0,alpha;alpha,0];%汉明失真矩阵Pu=[p,1-p];%信源概率统计delta=[0:p/1000:p];R=Hp(p)-Hp(delta/alpha);%汉明失真，率失真函数理论值eps1=0.001;%精度为eps1=0.001[R1,delta1]=R_delta(Pu,D,eps1);eps2=0.000000001;%精度为eps2=0.000000001[R2,delta2]=R_delta(Pu,D,eps2);plot(delta,R,'b',delta1,R1,'r-.',delta2,R2,'m-*');legend('理论值','eps1=0.001','eps2=0.000000001')测试2：对于不同的概率输入，测试程序如下：clearallalpha=2;p1=0.5;p2=0.2;D=[0,alpha;alpha,0];%汉明失真矩阵Pu1=[p1,1-p1];%信源概率统计Pu2=[p2,1-p2];%信源概率统计delta1=[0:p1/1000:p1];delta2=[0:p2/1000:p2];eps=0.00001;%精度为eps=0.00001[R1,delta1]=R_delta(Pu1,D,eps);[R2,delta2]=R_delta(Pu2,D,eps);plot(delta1,R1,'b-*',delta2,R2,'r-^');legend('p1=0.5','p2=0.2')Figure1-测试1图形Figure2-测试2图形