Matlab源程序.docx
- 文档编号:3648138
- 上传时间:2022-11-24
- 格式:DOCX
- 页数:22
- 大小:21.39KB
Matlab源程序.docx
《Matlab源程序.docx》由会员分享,可在线阅读,更多相关《Matlab源程序.docx(22页珍藏版)》请在冰豆网上搜索。
Matlab源程序
附录Matlab源程序
附录A信息熵
%函数说明:
%
%H=entropy(P,r)为信息熵函数%
%P为信源的概率矢量,r为进制数%
%H为信息熵%
%******************************%
functionH=entropy(P,r)
if(length(find(P<=0))~=0)
error('Notaprob.vector,negativecomponent');%判断是否符合概率分布条件
end
if(abs(sum(P)-1)>10e-10)
error('Notaprob.vector,componentdonotaddupto1');
end
H=(sum(-P.*log2(P)))/(log2(r)+eps);
附录B离散无记忆信道容量的迭代计算
%信道容量C的迭代算法%
%函数说明:
%
%[CC,Paa]=ChannelCap(P,k)为信道容量函数%
%变量说明:
%
%P:
输入的正向转移概率矩阵,k:
迭代计算精度%
%CC:
最佳信道容量,Paa:
最佳输入概率矩阵%
%Pa:
初始输入概率矩阵,Pba:
正向转移概率矩阵%
%Pb:
输出概率矩阵,Pab:
反向转移概率矩阵%
%C:
初始信道容量,r:
输入符号数,s:
输出符号数%
%**************************************************%
function[CC,Paa]=ChannelCap(P,k)
%提示错误信息
if(length(find(P<0))~=0)
error('Notaprob.vector,negativecomponent');%判断是否符合概率分布条件
end
if(abs(sum(P')-1)>10e-10)
error('Notaprob.vector,componentdonotaddupto1')%判断是否符合概率和为1
end
%1)初始化Pa
[r,s]=size(P);
Pa=(1/(r+eps))*ones(1,r);
sumrow=zeros(1,r);
Pba=P;
%2)进行迭代计算
n=0;
C=0;
CC=1;
whileabs(CC-C)>=k
n=n+1;
%
(1)先求Pb
Pb=zeros(1,s);
forj=1:
s
fori=1:
r
Pb(j)=Pb(j)+Pa(i)*Pba(i,j);
end
end
%
(2)再求Pab
suma=zeros(1,s);
forj=1:
s
fori=1:
r
Pab(j,i)=Pa(i)*Pba(i,j)/(Pb(j)+eps);
suma(j)=suma(j)+Pa(i)*Pba(i,j)*log2((Pab(j,i)+eps)/(Pa(i)+eps));
end
%3)求信道容量C
C=sum(suma);
%4)求下一次Pa,即Paa
L=zeros(1,r);
sumaa=0;
fori=1:
r
forj=1:
s
L(i)=L(i)+Pba(i,j)*log(Pab(j,i)+eps);
end
a(i)=exp(L(i));
end
sumaa=sum(a);
fori=1:
r
Paa(i)=a(i)/(sumaa+eps);
end
%5)求下一次C,即CC
CC=log2(sumaa);
Pa=Paa;
end
%打印输出结果
s0='很好!
输入正确,迭代结果如下:
';
s1='最佳输入概率分布Pa:
';
s2='信道容量C:
';
s3='迭代次数n:
';
s4='输入符号数r:
';
s5='输出符号数s:
';
s6='迭代计算精度k:
';
fori=1:
r
B{i}=i;
end
disp(s0);
disp(s1),disp(B),disp(Paa);
disp(s4),disp(r);
disp(s5),disp(s);
disp(s2),disp(CC);
disp(s6),disp(k);
disp(s3),disp(n);
附录CShannon编码
%函数说明:
%
%[p,x]=array(P,X)为按降序排序的函数%
%P为信源的概率矢量,X为概率元素的下标矢量%
%p为排序后返回的信源的概率矢量%
%x为排序后返回的概率元素的下标矢量%
%*******************************************%
function[p,x]=array(P,X)
P=[P;X];
[l,n]=size(P);
fori=1:
n
max=P(1,i);
maxN=i;
MAX=P(:
i);
forj=i:
n
if(max
MAX=P(:
j);
max=P(1,j);
maxN=j;
end
end
if(maxN>1)
if(i fork=(maxN-1): -1: i P(: k+1)=P(: k); end end end P(: i)=MAX; end p=P(1,: ); x=P(2,: ); %shannon编码生成器% %函数说明: % %[W,L,q]=shannon(p)为shannon编码函数% %p为信源的概率矢量,W为编码返回的码字% %L为编码返回的平均码字长度,q为编码效率% %*****************************************% function[W,L,q]=shannon(p) %提示错误信息 if(length(find(p<=0))~=0) error('Notaprob.vector,negativecomponent');%判断是否符合概率分布条件 end if(abs(sum(p)-1)>10e-10) error('Notaprob.vector,componentdonotaddupto1')%判断是否符合概率和为1 end %1)排序 n=length(p); x=1: n; [p,x]=array(p,x); %2)计算代码组长度l l=ceil(-log2(p)); %3)计算累加概率P P (1)=0; n=length(p); fori=2: n P(i)=P(i-1)+p(i-1); end %4)求得二进制代码组W %a)将十进制数转为二进制数 fori=1: n forj=1: l(i) temp(i,j)=floor(P(i)*2); P(i)=P(i)*2-temp(i,j); end end %b)给W赋ASCII码值,用于显示二进制代码组W fori=1: n forj=1: l(i) if(temp(i,j)==0) W(i,j)=48; else W(i,j)=49; end end end L=sum(p.*l);%计算平均码字长度 H=entropy(p,2);%计算信源熵 q=H/L;%计算编码效率 %打印输出结果 fori=1: n B{i}=i; end [n,m]=size(W); TEMP=32*ones(n,6); W=[W,TEMP]; W=W'; [n,m]=size(W); W=reshape(W,1,n*m); W=sprintf('%s',W); s0='很好! 输入正确,编码结果如下: '; s1='Shannon编码所得码字W: '; s2='Shannon编码平均码字长度L: '; s3='Shannon编码的编码效率q: '; disp(s0); disp(s1),disp(B),disp(W); disp(s2),disp(L); disp(s3),disp(q); 附录DFano编码 %函数说明: % %[next_P,next_index,code_num]=compare(current_P,current_index)% % 为比较函数,主要用于信源符号的分组% %current_P为当前分组的信源的概率矢量% %current_index为当前分组的信源的下标% %next_P为返回的下次分组的信源的概率矢量% %next_index为返回的下次分组的信源的下标% %code_num为返回的ASCII值% %*********************************************************************% function[next_P,code_num,next_index]=compare(current_P,current_index); n=length(current_P); add (1)=current_P (1); %1)求概率的依次累加和 fori=2: n add(i)=0; add(i)=add(i-1)+current_P(i); end %2)求概率和最接近的两小组 s=add(n); fori=1: n temp(i)=abs(s-2*add(i)); end [c,k]=min(temp); %3)对分组的信源赋ASCII值 if(current_index<=k) next_index=current_index; code_num=48; next_P=current_P(1: k); else next_index=current_index-k; code_num=49; next_P=current_P((k+1): n); end %fano编码生成器% %函数说明: % %[W,L,q]=fano(P)为fano编码函数% %P为信源的概率矢量,W为编码返回的码字% %L为编码返回的平均码字长度,q为编码效率% %*****************************************% function[W,L,q]=fano(P) %提示错误信息 if(length(find(P<=0))~=0) error('Notaprob.vector,negativecomponent');%判断是否符合概率分布条件 end if(abs(sum(P)-1)>10e-10) error('Notaprob.vector,componentdonotaddupto1')%判断是否符合概率和为1 end %1)排序 n=length(P); x=1: n; [P,x]=array(P,x); %2)将信源符号分组并得到对应的码字 fori=1: n current_index=i; j=1; current_P=P; while1 [next_P,code_num,next_index]=compare(current_P,current_index); current_index=next_index; current_P=next_P; W(i,j)=code_num; j=j+1; if(length(current_P)==1) break; end end l(i)=length(find(abs(W(i,: ))~=0));%得到各码字的长度 end L=sum(P.*l);%计算平均码字长度 H=entropy(P,2);%计算信源熵 q=H/L;%计算编码效率 %打印输出结果 fori=1: n B{i}=i; end [n,m]=size(W); TEMP=32*ones(n,5); W=[W,TEMP]; W=W'; [n,m]=size(W); W=reshape(W,1,n*m); W=sprintf('%s',W); s0='很好! 输入正确,编码结果如下: '; s1='Fano编码所得码字W: '; s2='Fano编码平均码字长度L: '; s3='Fano编码的编码效率q: '; disp(s0); disp(s1),disp(B),disp(W); disp(s2),disp(L); disp(s3),disp(q); 附录EHuffman编码 Huffman编码 (1) %huffman编码生成器% %函数说明: % %[W,L,q]=huffman(P)为huffman编码函数% %P为信源的概率矢量,W为编码返回的码字% %L为编码返回的平均码字长度,q为编码效率% %*****************************************% function[W,L,q]=huffman(P) if(length(find(P<=0))~=0) error('Notaprob.vector,negativecomponent');%判断是否符合概率分布条件 end if(abs(sum(P)-1)>10e-10) error('Notaprob.vector,componentdonotaddupto1')%判断是否符合概率和为1 end n=length(P);%计算输入元素个数 p=P; mark=zeros(n-1,n);%mark为n-1行、n列矩阵,用来记录每行最小两概率叠加后概率排列次序 %1)确定概率大小值的排列,得到mark矩阵。 fori=1: n-1 [p,num]=sort(p);%对输入元素排序并纪录 mark(i,: )=[num(1: n-i+1),zeros(1,i-1)]; p=[p (1)+p (2),p(3: n),1]; end %2)生成一个n-1行、n1(n×n)列矩阵table,每行可看做n个段, %每段长为n,记录一个码字(每个码字的长度不会超过n)。 fori=1: n-1 table(i,: )=blanks(n*n); end %3)计算各个元素码字,循环n-2次,决定矩阵table %从倒数第二行开始到第一行的每段的码字值,到编码表格table table(n-1,n)='1';%小值赋1 table(n-1,2*n)='0';%大值赋0 fori=2: n-1 table(n-i,1: n-1)=table(n-i+1,n*(find(mark(n-i+1,: )==1))... -(n-2): n*(find(mark(n-i+1,: )==1)));%按mark的记录依次赋值 table(n-i,n)='1'; table(n-i,n+1: 2*n-1)=table(n-i,1: n-1); table(n-i,2*n)='0'; forj=1: i-1 table(n-i,(j+1)*n+1: (j+2)*n)=table(n-i+1,... n*(find(mark(n-i+1,: )==j+1)-1)+1: n*find(mark(n-i+1,: )==j+1)); %按mark的记录依次赋值 end %4)得到编码后的码字 fori=1: n W(i,1: n)=table(1,n*(find(mark(1,: )==i)-1)+1: find(mark(1,: )==i)*n); l(i)=length(find(abs(W(i,: ))~=32)); end L=sum(P.*l);%计算平均码字长度 H=entropy(P,2);%计算信源熵 q=H/L;%计算编码效率 %打印输出结果 fori=1: n B{i}=i; end [m,n]=size(W); TEMP=blanks(m); W=[W,TEMP',TEMP',TEMP']; [m,n]=size(W); W=reshape(W',1,m*n); s0='很好! 输入正确,编码结果如下: '; s1='Huffman编码所得码字W: '; s2='Huffman编码平均码字长度L: '; s3='Huffman编码的编码效率q: '; disp(s0); disp(s1),disp(B),disp(W); disp(s2),disp(L); disp(s3),disp(q); Huffman编码 (2) %huffman编码生成器% %函数说明: % %[W,L,V,q]=huffman_better(P)为huffman_better编码函数% %P为信源的概率矢量,W为编码返回的码字,V为码字的方差% %L为编码返回的平均码字长度,q为编码效率% %*******************************************************% function[W,L,V,q]=huffman_better(P) if(length(find(P<=0))~=0) error('Notaprob.vector,negativecomponent');%判断是否符合概率分布条件 end if(abs(sum(P)-1)>10e-10) error('Notaprob.vector,componentdonotaddupto1')%判断是否符合概率和为1 end n=length(P);%计算输入元素个数 p=P; mark=zeros(n-1,n);%mark为n-1行、n列矩阵,用来记录每行最小两概率叠加后概率排列次序 %1)确定概率大小值的排列,得到mark矩阵。 t=1; fori=1: n-1 [p,num]=sort(p);%对输入元素排序并纪录 if(i~=1) if(count~=0) k=max(a(t,: )); fors=count: -1: 1 num(k-s)=num(k-s+1);%若有与新项相等的项,则将新项尽可能排列在其右侧 end num(k)=1; end t=t+1; end mark(i,: )=[num(1: n-i+1),zeros(1,i-1)]; p=[p (1)+p (2),p(3: n),1];%前两项求和得新项 count=0;%用于计数 forj=2: n-i if(p (1)==p(j))%判断p中是否有与求和后的新项相等的项 count=count+1; a(t,count)=j; end end end %2)生成一个n-1行、n1(n×n)列矩阵table,每行可看做n个段, %每段长为n,记录一个码字(每个码字的长度不会超过n)。 fori=1: n-1 table(i,: )=blanks(n*n); end %3)计算各个元素码字,循环n-2次,决定矩阵table %从倒数第二行开始到第一行的每段的码字值,到编码表格table table(n-1,n)='1';%小值赋1 table(n-1,2*n)='0';%大值赋0 fori=2: n-1 table(n-i,1: n-1)=table(n-i+1,n*(find(mark(n-i+1,: )==1))... -(n-2): n*(find(mark(n-i+1,: )==1)));%按mark的记录依次赋值 table(n-i,n)='1'; table(n-i,n+1: 2*n-1)=table(n-i,1: n-1); table(n-i,2*n)='0'; forj=1: i-1 table(n-i,(j+1)*n+1: (j+2)*n)=table(n-i+1,... n*(find(mark(n-i+1,: )==j+1)-1)+1: n*find(mark(n-i+1,: )==j+1)); %按mark的记录依次赋值 end end %4)得到编码后的码字 fori=1: n W(i,1: n)=table(1,n*(find(mark(1,: )==i)-1)+1: find(mark(1,: )==i)*n); l(i)=length(find(abs(W(i,: ))~=32)); end L=sum(P.*l);%计算平均码字长度 H=entropy(P,2);%计算信源熵 V=sum(P.*((l-L).^2));%计算码字的方差,以判断编码方法的优劣 q=H/L;%计算编码效率 %打印输出结果 fori=1: n B{i}=i; end [m,n]=size(W); TEMP=blanks(m); W=[W,TEMP',TEMP',TEMP']; [m,n]=size(W); W=reshape(W',1,m*n); s0='很好! 输入正确,编码结果如下: '; s1='Huffman编码所得码字W: '; s2='Huffman编码平均码字长度L: '; s3='Huffman编码所得码字W的方差V: '; s4='Huffman编码的编码效率q: '; disp(s0); disp(s1),disp(B),disp(W); disp(s2),disp(L); disp(s3),disp(V); disp(s4),disp(q) 附录F信息率失真函数的迭代计算 %信息率失真函数的迭代计算,迭代精度取为10^(-7)% %在信源的输入概率分布Pa和失真矩阵d已知的条件下求出信息率失真函数% %函数说明: % %[Pba,Rmin,Dmax,Smax]=RateDF(Pa,d,S)为信息率失真函数% %变量说明: % %Pa: 信源的输入概率矩阵,d: 失真矩阵,S: 拉氏乘子% %Pba: 最佳正向转移概率矩阵,Smax: 最大拉氏乘子% %Rmin: 最小信息率,Dmax: 允许的最大失真度% %Pb: 信源的输出概率矩阵,D: 允许的失真度,R: 信息率% %r: 输入信源数,s: 输出信源数% %*************************************************************
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Matlab 源程序
![提示](https://static.bdocx.com/images/bang_tan.gif)