信息與編碼實(shí)驗(yàn)二、容量代價(jià)函數(shù)的計(jì)算一、問題描述假定一個(gè)離散無記憶信源（DMC）的輸入符號(hào)集為AX={1,2,...,r},輸出符號(hào)集為AY={1,2,,s}，信道的轉(zhuǎn)移概率矩陣P(Y|X):{p(y|x)}；該信道不是免費(fèi)的，其輸入代價(jià)為b(x)。利用Matlab畫出容量代價(jià)函數(shù)C()的曲線圖。二、實(shí)驗(yàn)環(huán)境計(jì)算機(jī)、Windows2000或以上、Matlab6.5或以上三、實(shí)驗(yàn)?zāi)康?.了解容量代價(jià)函數(shù)性質(zhì)、意義。2.掌握簡(jiǎn)單的容量代價(jià)函數(shù)計(jì)算方法；3.掌握使用Matlab實(shí)現(xiàn)一般容量代價(jià)函數(shù)的計(jì)算方法；4.掌握Matlab求解非線性方程組的方法。四、實(shí)驗(yàn)原理1．C()的一般計(jì)算方法：計(jì)算min是很容易的。；當(dāng)minmax時(shí)，設(shè)試驗(yàn)信源X達(dá)到了C()，即C()=max{I(X;Y):E([b(x)]=}。在數(shù)學(xué)上，就是在約束條件下：(1)并且P(x)0(2)的約束下求平均互信息量的條件極大值。為此引入待定常數(shù)S和，并作輔助函數(shù)(3)其中由得，(4)上式兩邊同乘以p(x)并對(duì)x求和。C()=loge+S+(5)顯然(4)提供r個(gè)方程，再加上(1)(2)的方程，共r+2個(gè)方程；而有r個(gè)未知數(shù)p(x)、及S、，共r+2未知數(shù)，顯然可以求解。一般來說，S、是的函數(shù)。(1)(2)(4)組成的r+2個(gè)方程是一個(gè)非常復(fù)雜的非線性方程組。(4)式可以寫成矩陣形式:2．一般非線性方程組求解方法：下面介紹求非線性方程組一組是實(shí)根的梯度法：設(shè)非線性方程組為 并定義目標(biāo)函數(shù)為： 則梯度法的計(jì)算過程如下：選取一組初值x0,x1,...,xn-1計(jì)算目標(biāo)函數(shù)值F若F<，則X=(x0,x1,...,xn-1)即為方程組的一組實(shí)根，過程結(jié)束；否則繼續(xù)。計(jì)算目標(biāo)函數(shù)在(x0,x1,...,xn-1)的偏導(dǎo)數(shù)再計(jì)算計(jì)算其中=F/D。重復(fù)(2)~(5)直到滿足精度要求為止。在上述過程中，如果D=0，則說明遇到了目標(biāo)函數(shù)的局部極值點(diǎn)，此時(shí)可改變初值在試試。在我們的問題中，不會(huì)遇到該問題。另在計(jì)算中我們限制一個(gè)最大循環(huán)次數(shù)，以防D過小，使計(jì)算時(shí)間過長(zhǎng)的情況。五、實(shí)驗(yàn)過程（1）先從理論上計(jì)算r=s=2，傳輸錯(cuò)誤率為p，b(0)=0,b(1)=1的BSC信道的容量代價(jià)函數(shù)。這里min=0，簡(jiǎn)化信道只有一個(gè)輸入0，因此，有P{X=1}=，P{X=0}=1-，而C()=I(X;Y)=H(Y)-H(Y|X)，由于H（x）在x=0.5時(shí)達(dá)到極大值，故=0.5時(shí)H(Y)達(dá)到最大值，C()的公式如下。（2）對(duì)一般性的DMC信道，計(jì)算容量代價(jià)函數(shù)C()非常困難，我們借助matlab2012B來求解，具體代碼見附錄，這里只貼出代入數(shù)值后的結(jié)果。給定：r=s=2；p(y=0|x=0)=p(y=1|x=1)0.4,p=(y=1|x=0)=p(y=0|x=1)=0.6;b(0)=0,b(1)=1。由matlab計(jì)算繪制的圖像如下圖5.1Pe=0.6時(shí)容量代價(jià)函數(shù)圖像以下是誤碼率為0.8和0.55時(shí)的容量代價(jià)函數(shù)圖5.2容量代價(jià)函數(shù)Pe=0.85.3容量代價(jià)函數(shù)0.55從圖中不難的值C(β)的性質(zhì)：

在βmin≤β≤

βmax內(nèi)，是一個(gè)嚴(yán)格上升函數(shù)。

C(βmin)=Cmin，C(βmax)=Cmax。

當(dāng)β≥βmax時(shí)，C(β)是一個(gè)常數(shù)Cmax。

因此C(β)=max{I(X;Y):E[b(x)]=

β}，βmin≤β≤βmax。這說明，達(dá)到C(β)的試驗(yàn)信號(hào)，可以在E[b(x)]=β的區(qū)域中找到，實(shí)際上也只能在這個(gè)區(qū)域內(nèi)找到。（2）接下來測(cè)試程序運(yùn)行時(shí)間，誤碼率為0.55，代碼見附錄測(cè)試結(jié)果如下結(jié)果顯示，fsolve函數(shù)進(jìn)行了分次迭代，共循環(huán)了62次，耗時(shí)0.324402秒，那么平均每一次迭代耗費(fèi)時(shí)間0.005232秒。3.結(jié)果分析由上面的圖像可見，誤碼率越接近0.5，即輸入越接近等概，容量代價(jià)函數(shù)最大值越小，。根據(jù)上面的推論，當(dāng)p=0.5時(shí)，容量代價(jià)函數(shù)為0，這也符合實(shí)驗(yàn)的結(jié)果。六、附錄代碼（1）繪制容量代價(jià)函數(shù)圖像functionc=C(beta,b,Py_x)sizeB=size(b);sizePy_x=size(Py_x);if~(sizeB(2)==sizePy_x(1))disp('參數(shù)b的列數(shù)應(yīng)該與Py_x的行數(shù)相同');c=-1;return;endm=sizePy_x(1);n=sizePy_x(2);sumPy=sum(Py_x');if~(all(Py_x>=0)&all(sumPy==1))disp('參數(shù)Py_x每個(gè)數(shù)都不為負(fù)，且每一行之和應(yīng)為1');c=-1;return;endminBeta=min(b);ifbeta<minBetadisp('參數(shù)beta太小');c=-1;return;endx0=ones(1,m+2)./m;options=optimset('TolFun',0.000001,'Display','iter');%設(shè)置收斂條件及查看迭代次數(shù)r=fsolve(@c_bate,x0,options,beta,b,Py_x);s=r(m+1);u=r(m+2);c=log2(exp(1))+s*beta+u;EndfunctionF=c_bate(x,beta,b,Py_x)sizePy_x=size(Py_x);m=sizePy_x(1);n=sizePy_x(2);F=zeros(m+2,1);i=1:m;P_y=x(1,1:m)*Py_x;TP_y=zeros([m,n]);forj=1:mTP_y(j,:)=P_y;endTP_y=TP_y';t=Py_x(i,:)'./TP_y;t(find(t==0))=1;F(i)=sum(Py_x(i,:)'.*log2(t))-log2(exp(1))-x(m+1)*b(i)-x(m+2);F(m+1)=1-sum(x(1,1:m));F(m+2)=beta-sum(x(1,1:m).*b);（2）這里是畫容量代價(jià)函數(shù)的主體程序，也是運(yùn)行的主要函數(shù)。此函數(shù)里調(diào)用了以上兩個(gè)函數(shù)，所以此函數(shù)的運(yùn)行時(shí)間就是整體的運(yùn)行時(shí)間。function[minb,maxb]=plotCBeta(b,Py_x)tic;%開始計(jì)算程序運(yùn)行時(shí)間minBeta=min(b);maxBeta=max(b);beta=minBeta:0.1:maxBeta;l=length(beta);