信息论实验课上的作业，其中包括【香农编码、哈弗曼编码、算术编码和费诺编码】虽然有很大一部分都是网上借鉴别人的，但是自己琢磨了几次也算是看懂了原理，动手写了一些自己理解的东西，索性就记录一下。
1.香农编码
编码步骤：
将信源发出的q个消息符号按其概率大小的次序排列；
计算出各个符号的累加概率。
按下式计算第i个消息对应的码字长度。
将累加概率（十进制小数）变换成二进制小数，根据码长Li取小数点后Li个二进制符号作为第i个消息的码字。
代码：
【shannon.m】
%香农编码
%input：信源分布p
%output:二进制编码code
function shannon(p)
n = length(p);%信源符号数
p = sort(p,'descend');%降序排列
F = zeros(1,n);%累加概率序列
for i=2:n
F(i)=F(i-1)+p(i-1);
end
l = zeros(1,n);%码长序列
for i=1:n
l(i)=ceil(-log2(p(i)));%向上取整
end
disp('最终编码结果：');
for i=1:n
f = max(l(i));
w=trans(F(i),f); %调用二进制转换编码函数
for j=1:l(i) %这个循环是指取前多少个作为码字
code(j)=w(j);
end
disp([num2str(p(i)),'->',num2str(code)]); %这句是作为输出语句
end【trans.m】
%十进制转换为二进制，应用于香农编码
%原理：乘二取整
%输入：需要编码的概率F，f为码长
%输出：编码后的码字
function w = trans( F,f )
for i=1:f
temp=F.*2;
if(temp<1)
w(i)=0;
F=temp;
else
F=temp-1;
w(i)=1;
end
end运行结果：
>> p=[0.20 0.19 0.18 0.17 0.15 0.10 0.01]
p =
0.2000 0.1900 0.1800 0.1700 0.1500 0.1000 0.0100
>> shannon(p)
最终编码结果：
0.2->0 0 0
0.19->0 0 1
0.18->0 1 1
0.17->1 0 0
0.15->1 0 1
0.1->1 1 1 0
0.01->1 1 1 1 1 1 02.哈弗曼编码（二元）
代码：
【myHuffman.m】
clc;
clear;
%输入格式：[p1 p2 p3 … pn]
%A = [0.1 0.18 0.4 0.05 0.06 0.1 0.07 0.04];
A = input('输入信源概率分布：\n')
n = length(A);
%----先判断输入格式是否正确
for i = 1:n
if A(i)<0
fprintf('概率小于0，请重新输入。\n');
A = input('输入信源概率分布：\n')
end
end
A = sort(A,'descend'); %降序排列
T = A;%设置一个当前编码序列
[m,n] = size(A);%行列数
B=zeros(n,n-1);%生成一个空的编码表
%-----编码表的第一列
for i=1:n
B(i,1)=T(i);
end
sum=B(i,1)+B(i-1,1);%最后两个概率相加
T(n-1)=sum;
T(n)=0;%把刚所求最后两个数之和放进当前编码序列，最后一个设置为0
T = sort(T,'descend');%然后重新进行排序
t = n-1;%当前序列数目
%-----编码表的其他列
for j=2:n-1 %从第二列开始
for i=1:t
B(i,j)=T(i);
end
K=find(T==sum); %找到当前序列中sum的位置，记在K中
B(n,j)=K(end);%从第二列开始，编码表中的最后一位元素是sum在当前序列中的位置
sum=(B(t-1,j)+B(t,j));%继续把最后两个元素相加记为sum
T(t-1)=sum;
T(t)=0;%把刚所求最后两个数之和放进当前编码序列，最后一个设置为0
T = sort(T,'descend'); %然后重新进行排序
t=t-1;%当前序列数目
end
%编码表
disp('编码表为：')
B
%-----开始编码
%先设置最后一列的两个元素编码为0和1
%这里写下划线是因为matlab会显示不出来0，所以只能加上下划线让它显示出来
%因为有下划线，所以后面计算码长的时候会除2
END1=sym('[_0,_1]');
END=END1;
t=3;
d=1;
for j=n-2:-1:1 %从倒数第二列开始依次对各列元素编码
for i=1:t-2
if i>1 & B(i,j)==B(i-1,j)
d=d+1;
else
d=1;
end
%把相加后的那个数设置为-1方便找到
B(B(n,j+1),j+1)=-1;
temp=B(:,j+1);
x=find(temp==B(i,j));
%给没有变过的值编码
END(i)=END1(x(d));
end
y=B(n,j+1);
%给变过的值编码，上面的补0，下面的补1
END(t-1)=[char(END1(y)),'_0'];
END(t)=[char(END1(y)),'_1'];
t=t+1;
END1=END;
end
%编码结果
disp('编码结果为：')
END
%-----计算编码效率p
for i=1:n
[a,b]=size(char(END(i)));
M(i)=b/2; %这里除以2是因为有下划线占了一半的码字长度
end
%计算平均码长
L1 = 0;
for i=1:n
L1 = L1 + M(i).*A(i);
end
%计算信息熵
H1=log2(A);
H=-A*(H1');
p=H/L1;
disp('编码效率为：')
p备注：= = 这个有个很大的bug就是0在matlab里居然显示不出来！我也又气又无奈……只能用下划线来让它表示出来了。
运行结果：
输入信源概率分布：
[1/6 1/4 5/12 1/6]
A =
0.1667 0.2500 0.4167 0.1667
编码表为：
B =
0.4167 0.4167 0.5833
0.2500 0.3333 0.4167
0.1667 0.2500 0
0.1667 2.0000 1.0000
编码结果为：
END =
[ _1, _0_1, _0_0_0, _0_0_1]
编码效率为：
p =
0.98503.算术编码
备注：哈哈哈哈这个是我自己写的！！（= = 虽然很简单但本菜鸟还是有自豪感的……）
代码：
【suanshu.m】
%算术编码程序
%输入：信源u、信源分布p、需要进行编码的序列s
%输入备注：输入格式：a=[0 1] p = [1/4 3/4] s = [1 0 1 1 1]
%输出：该序列的算术编码S
%输出备注：以序列所在区间的左端点值作为所求序列的码字
%输入备注：如果信源符号是字符型，输入方式为：a = 'ab'; s = 'bbbababb';
%其实符号是数字可以写成数组也可以写成字符类型，matlab都支持运算，如：a = '01' s = '1101'
%序列比较长的时候写成字符类型比较方便
function [S] = suanshu(a,p,s)
%设置初始值,序列为空，左端点为0，右端点为1，长度为1
Slow(1)=0;
Shigh(1)=1;
range=Shigh(1)-Slow(1);
%计算信源符号的数目和所需要编码的序列长度
n = length(a);
m = length(s);
%计算信源a中每个元素对应的左右端点low high
low(1) = 0;
for i = 1:n
if i == 1
high(i) = p(i);
else
low(i) = p(i-1) + low(i-1);
high(i) = low(i) + p(i);
end
end
%开始遍历所要编码序列中的每个元素
for i = 1:m
%寻找与当前元素相等的信源符号，把索引放入num中
num = find(a==s(i));
%开始计算添加新符号后的序列左右端点newSlow和newShigh
Slow(i+1) = Slow(i) + range*low(num);
Shigh(i+1) = Slow(i) + range*high(num);
range = Shigh(i+1) - Slow(i+1);
end
%最后以序列所在区间的左端点值作为所求序列的码字
%即输出最后计算得到的Slow(m+1)
S = Slow(m+1);运行结果：
>> a = '01'
a =
01
>> p = [1/8 7/8]
p =
0.1250 0.8750
>> s = '11111110111110'
s =
11111110111110
>> suanshu(a,p,s)
ans =
0.63124.费诺编码
备注：这个是朋友给我的= =，只看懂了一点，至于各种递归调用我已经晕了……
代码：
【f1.m】
function x=f1(i,j,p,r)
global x;
x=char(x);
if(j<=i)
return;
else
q=0;
%求累加概率
for t=i:j
q=p(t)+q;y(t)=q;
end
%累加概率值与该总概率值减该累加概率值之差取绝对值
for t=i:j
v(t)=abs(y(t)-(q-y(t)));
end
%取最小的一个值作为分界点
for t=i:j
if(v(t)==min(v))
%赋值
for k=i:t
x(k,r)='0';
end
for k=(t+1):j
x(k,r)='1';
end
d=t;
%递归
f1(i,d,p,r+1);
f2(d+1,j,p,r+1);
f1(d+1,j,p,r+1);
f2(i,d,p,r+1);
else
end
end
end
return;【f2.m】
function x=f2(i,j,p,r)
global x;
x=char(x);
if(j<=i)
return;
else
q=0;
%求累加概率值
for t=i:j
q=p(t)+q;y(t-i+1)=q;
end
%%累加概率值与该总概率值减该累加概率值之差取绝对值
for t=1:j-(i-1)
v(t)=abs(y(t)-(q-y(t)));
end
%取最小的一个值作为分界点
for t=1:j-(i-1)
if(v(t)==min(v))
d=t+i-1;
%赋值
for k=i:d
x(k,r)='0';
end
for k=(d+1):j
x(k,r)='1';
end
%递归
f2(d+1,j,p,r+1);
f1(i,d,p,r+1);
f2(i,d,p,r+1);
f1(d+1,j,p,r+1);
else
end
end
end
return;【myfano.m】
fprintf('请输入信源概率');
p=input('p=');
N = length(p);
s=0;l=0;H=0;
%计算信息熵
for i=1:N
H=H+(- p(i)*log2(p(i)));
end
tic;
for i=1:N-1 %按概率分布大小对信源排序
for j=i+1:N
if p(i)<p(j)
m=p(j);p(j)=p(i);p(i)=m;
end
end
end
x=f1(1,N,p,1);
%计算平均码长
for i=1:N
L(i)=length(find(x(i,:)));
l=l+p(i)*L(i);
end
%编码效率
n=H/l;
fprintf('编码效率:');disp(n)
fprintf('编码结果:\n');
disp(x)
fprintf('平均码长:K=');disp(l)运行结果：
>> myfano
请输入信源概率p=[0.2 0.19 0.18 0.17 0.15 0.1 0.01]
编码效率: 0.9521
编码结果:
00
010
011
10
110
1110
1111
平均码长:K= 2.7400本博客所有文章除特别声明外，均采用 CC BY-SA 3.0协议 。转载请注明出处！