matlab的编码大全

附录Matlab源程序

附录A 信息熵

% 函数说明：%

% H=entropy(P,r) 为信息熵函数%

% P为信源的概率矢量, r为进制数%

% H为信息熵%

%****************************** %

function H=entropy(P,r)

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');

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('Not a prob.vector,negative component'); % 判断是否符合概率分布条件end

294

信息论基础及应用

if (abs(sum(P')-1)>10e-10)

error('Not a prob.vector,component do not add up to 1') % 判断是否符合概率和为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;

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

% (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));

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);

for i=1:r

Paa(i)=a(i)/(sumaa+eps);

end

% 5)求下一次C,即CC

附录Matlab源程序295

CC=log2(sumaa);

Pa=Paa;

end

% 打印输出结果

s0='很好！输入正确，迭代结果如下：';

s1='最佳输入概率分布Pa:';

s2='信道容量C:';

s3='迭代次数n:';

s4='输入符号数r:';

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);

disp(s3),disp(n);

附录C Shannon编码

% 函数说明：%

% [p,x]=array(P,X) 为按降序排序的函数%

% P为信源的概率矢量，X为概率元素的下标矢量%

% p为排序后返回的信源的概率矢量%

% x为排序后返回的概率元素的下标矢量%

%*******************************************%

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

MAX=P(:,j);

max=P(1,j);

maxN=j;

end

296

信息论基础及应用

end

if (maxN>1)

if (i

for k=(maxN-1):-1:i

P(:, k+1)=P(:,k);

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('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

% 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);

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)

附录Matlab源程序297

temp(i,j)=floor(P(i)*2);

P(i)=P(i)*2-temp(i,j);

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

L=sum(p.*l); % 计算平均码字长度

H=entropy(p,2); % 计算信源熵

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='很好！输入正确，编码结果如下：';

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编码

% 函数说明：%

% [next_P,next_index,code_num]=compare(current_P,current_index) %

% 为比较函数，主要用于信源符号的分组%

% current_P为当前分组的信源的概率矢量%

298

信息论基础及应用

% 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)求概率的依次累加和

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));

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('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

附录Matlab源程序299

% 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);

current_index=next_index;

current_P=next_P;

W(i,j)=code_num;

j=j+1;

if (length(current_P)==1)

break;

end

l(i)=length(find(abs(W(i,:)) ~=0)); % 得到各码字的长度

end

L=sum(P.*l); % 计算平均码字长度

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);

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);

300

信息论基础及应用

附录E Huffman编码

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('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

n=length(P); % 计算输入元素个数

p=P;

mark=zeros(n-1,n); % mark为n-1行、n列矩阵，用来记录每行最小两概率叠加后概率排列次序% 1) 确定概率大小值的排列，得到mark矩阵。

for i=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）。

for i=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

for i=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';

for j=1:i-1

table(n-i,(j+1)*n+1:(j+2)*n)=table(n-i+1,...

附录Matlab源程序301

n*(find(mark(n-i+1,:)==j+1)-1)+1:n*find(mark(n-i+1,:)==j+1));

% 按mark的记录依次赋值

end

% 4)得到编码后的码字

for i=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; % 计算编码效率

% 打印输出结果

for i=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('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

302

信息论基础及应用

n=length(P); % 计算输入元素个数

p=P;

mark=zeros(n-1,n); % mark为n-1行、n列矩阵，用来记录每行最小两概率叠加后概率排列次序% 1) 确定概率大小值的排列，得到mark矩阵。

t=1;

for i=1:n-1

[p,num]=sort(p); % 对输入元素排序并纪录

if (i~=1)

if (count~=0)

k=max(a(t,:));

for s=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; % 用于计数

for j=2:n-i

if (p(1)==p(j)) % 判断p中是否有与求和后的新项相等的项

count=count+1;

a(t,count)=j;

end

% 2) 生成一个n-1行、n1(n×n)列矩阵table，每行可看做n个段，

% 每段长为n，记录一个码字（每个码字的长度不会超过n）。

for i=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

for i=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';

for j=1:i-1

table(n-i,(j+1)*n+1:(j+2)*n)=table(n-i+1,...

附录Matlab源程序303

n*(find(mark(n-i+1,:)==j+1)-1)+1:n*find(mark(n-i+1,:)==j+1));

% 按mark的记录依次赋值

end

% 4)得到编码后的码字

for i=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; % 计算编码效率

% 打印输出结果

for i=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：允许的最大失真度%

304

信息论基础及应用

% Pb：信源的输出概率矩阵，D：允许的失真度，R：信息率%

% r：输入信源数，s：输出信源数%

%*********************************************************************%

function [Pba,Rmin,Dmax,Smax]=RateDF(Pa,d,S)

% 提示错误信息

[r,s]=size(d);

if (length(find(Pa<=0)) ~=0)

error('Not a prob.vector, shoud be positive component!'); % 判断是否符合概率分布条件end

if (abs(sum(Pa)-1)>10e-10)

error('Not a prob.vector,component do not add up to 1!') % 判断是否符合概率和为1 end

if (r~=length(Pa))

error('The parameters do not match!'); % 判断参数是否一致

end

% 第一步

pba=[];

RS=[]; % R(S)函数初始化

DS=[]; % D(S)函数初始化

m=1; % m为S循环的次数

while (1) % 外层循环，对S 的循环

Pba(1:r,1:s,1)=1/s*ones(r,s); % 求信道正向转移矩阵Pba

% 第二步

for j=1:s

Pb(j,1)=0;

for i=1:r

Pb(j,1)=Pb(j,1)+Pa(i)*Pba(i,j,1);

% 求信源的输出概率矩阵,即Pb(j,1)

end

for i=1:r

temp(i)=0;

for j=1:s

temp(i)=temp(i)+Pb(j,1)*exp(S(m)*d(i,j));

% temp为临时项,求Pba(i,j,2)时表达式的分母

end

for i=1:r

for j=1:s

Pba(i,j,2)=(Pb(j,1)*exp(S(m)*d(i,j)))/temp(i);

end

D(1)=0;

for i=1:r

附录Matlab源程序305

for j=1:s

D(1)=D(1)+Pa(i)*Pba(i,j,1)*d(i,j); % 求D(1)

end

R(1)=0;

for i=1:r

for j=1:s

if (Pba(i,j,1) ~=0)

R(1)=R(1)+Pa(i)*Pba(i,j,1)*log2(Pba(i,j,1)/Pb(j,1)); % 求R(1)

end

n=2; % n为内层循环次数

while (1) % 内层循环，对精度的循环

% 第三步

for j=1:s

Pb(j,n)=0;

for i=1:r

Pb(j,n)=Pb(j,n)+Pa(i)*Pba(i,j,n); % 求输出的信源概率分布

end

for i=1:r

temp(i)=0;

for j=1:s

temp(i)=temp(i)+Pb(j,n)*exp(S(m)*d(i,j));

% temp为临时项,求Pba(i,j,n+1)时表达式的分母

end

for i=1:r

for j=1:s

if (temp(i) ~=0)

Pba(i,j,n+1)=(Pb(j,n)*exp(S(m)*d(i,j)))/temp(i);

% 求Pba(i,j,n+1)

end

% 第四步

D(n)=0;

for i=1:r

for j=1:s

D(n)=D(n)+Pa(i)*Pba(i,j,n)*d(i,j); % 求D(n)

end

R(n)=0;

306

信息论基础及应用

for i=1:r

for j=1:s

if (Pba(i,j,n) ~=0)

R(n)=R(n)+Pa(i)*Pba(i,j,n)*log2(Pba(i,j,n)/Pb(j,n));

% 求R(n)

end

% 判断差别是否在允许的精度范围之内

if (abs(R(n)-R(n-1))<=10^(-7)) % R(n)精度判断

if (abs(D(n)-D(n-1))<=10^(-7)) % D(n)精度判断

break;

end

n=n+1; % 内层循环次数加1

end % 内层循环结束

% 第五步

S(m+1)=S(m)+0.5;

% 第六步

if (abs(R(n)<10^(-7)))

if (m<=10)

disp('此时S的值为:'),disp(S(m));

error('初始拉氏乘子S取得大了，请取小些!');

% 判断拉氏乘子S的初始值是否合适

else

[k,l,q]=size(pba);

Pba=pba(:,:,q);

Rmin=min(RS);

Dmax=max(DS);

Smax=S(m-1);

break;

end

pba=[Pba(:,:,:)];

RS=[RS R(n)];

DS=[DS D(n)];

m=m+1; % 外层循环次数加1

end % 外层循环结束

% 打印输出结果

s0='很好！输入正确，迭代结果如下:';

s1='最佳转移概率分布Pba:';

s2='最小信息率Rmin:';

s3='最大失真度Dmax:';

s4='最大拉氏乘子Smax:';

附录Matlab源程序307

disp(s0);

disp(s1),disp(Pba);

disp(s2),disp(Rmin);

disp(s3),disp(Dmax);

disp(s4),disp(Smax);

% 画出信息率失真函数R(D)

plot(DS,RS)

xlabel('允许的失真度D')

ylabel('信息率失真函数R(D)')

title('信息率失真函数R(D)的曲线图')

参考文献

1 C.E.Shannon.A Mathematical Theory of Communication.B.S.T.J.vol.27, July, and Oct.1948

2 C.E.Shannon.Coding Theorems for a Discrete Source with a Fidelity Criterion.IRE Nat.Conv.Rec.part 4, 1959

3 C.E.Shannon.Two way Communication channels.Pro.4th Berkeley Symp.on Math.Statisties and Probability,

vol.I.pp.611—644, 1961

4 N. https://www.360docs.net/doc/2110850606.html,rmation Theory and Coding.McGraw-Hall, Inc., 1963

5 https://www.360docs.net/doc/2110850606.html,rmation Theory and Reliable Communication. John Wley & Sons., 1968

6 T.Berger.Rate Distortion Theory.Prentice-Hall, Inc., 1971

7 https://www.360docs.net/doc/2110850606.html,rmation Theory With Application. McGraw-Hill，Inc.，1977

8 D.S.Jones.Elementary Information Theory.Clarendon Press.Oxford，1979

9 Robert J.McEliece.The Theory of Information and Coding(第2版). 电子工业出版社，2002

10 Thomas M.Cover, Joy A.Thomas. Elements of Information Theory. 清华大学出版社，2003

11 周炯槃.信息理论基础.人民邮电出版社，1983

12 傅祖芸.信息论基础.电子工业出版社，1989

13 周荫清.信息理论基础.北京航空航天大学出版社，1993

14 钟义信.信息科学原理.北京邮电大学出版社，1996

15 朱雪龙.应用信息论基础.清华大学出版社，2001

16 傅祖芸.信息论——基础理论与应用.电子工业出版社，2001

17 姜丹.信息论与编码.中国科学技术大学出版社，2001

18 曹雪虹等.信息论与编码.北京邮电大学出版社，2001

19 沈世镒等.信息论与编码理论.科学出版社，2002

20 周荫清.信息理论基础（修订版）. 北京航空航天大学出版社，2002

21 余成波.信息论与编码.重庆大学出版社，2002

22 陈运等.信息论与编码.电子工业出版社，2002

23 仇佩亮.信息论与编码.高等教育出版社，2003

24 叶中行.信息论基础.高等教育出版社，2003

25 唐朝京等.信息论与编码基础.国防科技大学出版社，2003

26 吕锋等.信息理论与编码.人民邮电出版社，2004