Matlab源程序

function [p,x]=array(P,X)
P=[P;X];
[l,n]=size(P);
for i=1:n
max=P(1,i);
maxN=i;
MAX=P(:,i);
for j=i:n
if (max<P(1,j))
MAX=P(:,j);
max=P(1,j);
maxN=j;
end
end
if (maxN>1)
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); %计算平均码字长度
% Pb：输出概率矩阵，Pab：反向转移概率矩阵%
% C：初始信道容量，r：输入符号数，s：输出符号数%
%**************************************************%
function [CC,Paa]=ChannelCap(P,k)
%提示错误信息
if (length(find(P<0))~=0)
P(i)=P(i)*2-temp(i,j);
end
end
% b)给W赋ASCII码值，用于显示二进制代码组W
for i=1:n
for j=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); %计算信源熵
附录B 离散无记忆信道容量的迭代计算
%信道容量C的迭代算法%
%函数说明：%
% [CC,Paa]=ChannelCap(P,k)为信道容量函数%
%变量说明：%
% P：输入的正向转移概率矩阵，k：迭代计算精度%
% CC：最佳信道容量，Paa：最佳输入概率矩阵%
% Pa：初始输入概率矩阵，Pba：正向转移概率矩阵%
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);
附录D Fano编码
%函数说明：%
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;
while abs(CC-C)>=k
n=n+1;
% (1)先求Pb
Pb=zeros(1,s);
for j=1:s
for i=1:r
Pb(j)=Pb(j)+Pa(i)*Pba(i,j);
end
end
% (2)再求Pab
suma=zeros(1,s);
for j=1:s
for i=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));
function [W,L,q]=huffman(P)
if (length(find(P<=0))~=0)
error('Not a prob.vector,negative component'); %判断是否符合概率分布条件
next_P=current_P((k+1):n);
end
% fano编码生成器%
%函数说明：%
% [W,L,q]=fano(P)为fano编码函数%
% P为信源的概率矢量，W为编码返回的码字%
% L为编码返回的平均码字长度，q为编码效率%
%*****************************************%
% L为编码返回的平均码字长度，q为编码效率%
%*****************************************%
function [W,L,q]=shannon(p)
%提示错误信息
if (length(find(p<=0))~=0)
error('Not a prob.vector,negative component'); %判断是否符合概率分布条件
end
% 1)排序
n=length(P);
x=1:n;
[P,x]=array(P,x);
% 2)将信源符号分组并得到对应的码字
for i=1:n
current_index=i;
j=1;
current_P=P;
while 1
[next_P,code_num,next_index]=compare(current_P,current_index);
error('Not a prob.vector,negative component'); %判断是否符合概率分布条件
end
if (abs(sum(P')-1)>10e-10)
error('Not a prob.vector,component do not add up to 1') %判断是否符合概率和为1
end
if (abs(sum(p)-1)>10e-10)
error('Not a prob.vector,component do not add up to 1') %判断是否符合概率和为1
end
% 1)排序
n=length(p);
x=1:n;
[p,x]=array(p,x);
% 2)计算代码组长度l
if (i<n)
for k=(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为编码返回的码字%
s0='很好！输入正确，编码结果如下：';
s1='Fano编码所得码字W:';
s2='Fanoຫໍສະໝຸດ Baidu码平均码字长度L:';
s3='Fano编码的编码效率q:';
disp(s0);
disp(s1),disp(B),disp(W);
disp(s2),disp(L);
disp(s3),disp(q);
附录E Huffman编码
n=length(current_P);
add(1)=current_P(1);
% 1)求概率的依次累加和
for i=2:n
add(i)=0;
add(i)=add(i-1)+current_P(i);
end
% 2)求概率和最接近的两小组
s=add(n);
for i=1:n
temp(i)=abs(s-2*add(i));
for i=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:';
disp(s3),disp(n);
附录C Shannon编码
%函数说明：%
% [p,x]=array(P,X)为按降序排序的函数%
% P为信源的概率矢量，X为概率元素的下标矢量%
% p为排序后返回的信源的概率矢量%
% x为排序后返回的概率元素的下标矢量%
%*******************************************%
附录Matlab源程序
附录A 信息熵
%函数说明：%
% H=entropy(P,r)为信息熵函数%
% P为信源的概率矢量, r为进制数%
% H为信息熵%
%******************************%
function H=entropy(P,r)
if (length(find(P<=0))~=0)
function [W,L,q]=fano(P)
%提示错误信息
if (length(find(P<=0))~=0)
error('Not a prob.vector,negative component'); %判断是否符合概率分布条件
end
if (abs(sum(P)-1)>10e-10)
error('Not a prob.vector,component do not add up to 1') %判断是否符合概率和为1
s5='输出符号数s:';
s6='迭代计算精度k:';
for i=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);
% next_index为返回的下次分组的信源的下标%
% code_num为返回的ASCII值%
%*********************************************************************%
function [next_P,code_num,next_index]=compare(current_P,current_index);
l=ceil(-log2(p));
% 3)计算累加概率P
P(1)=0;
n=length(p);
for i=2:n
P(i)=P(i-1)+p(i-1);
end
% 4)求得二进制代码组W
% a)将十进制数转为二进制数
for i=1:n
for j=1:l(i)
temp(i,j)=floor(P(i)*2);
% [next_P,next_index,code_num]=compare(current_P,current_index)%
%为比较函数，主要用于信源符号的分组%
% current_P为当前分组的信源的概率矢量%
% current_index为当前分组的信源的下标%
% next_P为返回的下次分组的信源的概率矢量%
error('Not a prob.vector,negative component'); %判断是否符合概率分布条件
end
if (abs(sum(P)-1)>10e-10)
error('Not a prob.vector,component do not add up to 1');
end
H=(sum(-P.*log2(P)))/(log2(r)+eps);
end
% 3)求信道容量C
C=sum(suma);
% 4)求下一次Pa,即Paa
L=zeros(1,r);
sumaa=0;
for i=1:r
for j=1:s
L(i)=L(i)+Pba(i,j)*log(Pab(j,i)+eps);
end
a(i)=exp( L(i));
end
sumaa=sum(a);
Huffman编码（1）
% huffman编码生成器%
%函数说明：%
% [W,L,q]=huffman(P)为huffman编码函数%
% P为信源的概率矢量，W为编码返回的码字%
% L为编码返回的平均码字长度，q为编码效率%
%*****************************************%
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;
q=H/L; %计算编码效率
%打印输出结果
for i=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='很好！输入正确，编码结果如下：';
H=entropy(P,2); %计算信源熵
q=H/L; %计算编码效率
%打印输出结果
for i=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);