信源编码的Matlab实现：费诺编码、香农编码、哈夫曼编码、算术编码、LZ编码、连续信源编码量化 M02014167方正，P02114173陆康超

Infor_Theory_AHU

已于 2023-07-18 16:51:18 修改

阅读量1.1w

点赞数 36

文章标签： matlab 开发语言信息论编码算法人工智能机器学习

于 2023-07-18 15:19:16 首次发布

本文链接：https://blog.youkuaiyun.com/lily_ahu/article/details/131787889

版权

文章目录

一、概要

信源编码是为了提高通信系统的有效性，通过压缩信息的冗余度来实现。与之对应的信道编码，则是增加信息冗余度来实现，从而提高信息传输的安全性和信源与信道的适配性。信源分为离散信源和连续信源，要先判断信源的类型，然后采取相应的编码方法，但编码途径都是解除序列符号之间的相关性与概率均匀化。本文介绍几种基本的离散信源编码：费诺编码、香农编码、哈夫曼编码、算术编码及LZ-78编码并给出Matlab实现代码，最后简要介绍连续信源编码的量化方法。

二、离散信源编码

1.1 费诺编码

费诺码即‘’香农-范诺编码‘’（Shannon–Fano coding）是一种基于一组符号集及其出现的或然率（估量或测量所得）构建前缀码的技术，它是一种无失真信源编码方法，由Ralph Fano 于1949年在香农的"A Mathematical Theory of Communication" (1948)文章基础上提出。费诺码在编码效率上，它并不能与霍夫曼编码一样实现编码（code word）长度的最低期望；然而，与霍夫曼编码不同的是，它确保了所有的编码长度在一个理想的理论范围之内。

费诺码的构建过程是递归的，它根据符号的概率分布来确定每个符号所对应的编码。具体步骤如下：

将概率按从大到小的顺序排列，不失一般性，令 $p(a_{1})\geqslant p(a_{2})\geqslant p(a_{n})$ ；
按编码进制数将概率分组，使每组概率和尽可能接近或相等。如编二进制码就分成两组，编m进制码就分成m组；
给每组分配一位码元；
将每一分组再按同样原则划分，重复步骤2和3，直至概率不再可分为止。

费诺码适用于符号概率已知的情况，且在概率分布不平衡时表现较好，可以达到较高的压缩效率。然而，费诺码并不是最优的编码方法，它不能达到信息熵的下界，即香农编码的压缩效率。在实际应用中，霍夫曼编码等其他编码方法通常比费诺码更常用，因为它们可以更接近信息熵的极限。

程序

迭代函数

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

编码函数

function[W,L,q]=fano(P)
 %1)排序
 n=length(P);
 x=1:n;
 % 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
     end
     l(i)=length(find(abs(W(i,:))~=0));%得到各码宇的长度
 end
 L=sum(P.*l);%计算平均码字长度
 H=sum(P.*(-log2(P)));%计算信源熵
 q=H/L;  %计算编码效率
  
 %打印输出结果
 for i=1:n
     B{i}=x(i);
 end
 [m,n]=size(W);
 TEMP=blanks(m);
 W=[W,TEMP',TEMP',TEMP'];
 [m,n]=size(W);
 W=reshape(W',1,m*n);
  
 fprintf('信源符号出现的概率为:\n');
 disp(P);
 fprintf('Fano编码所得码字W：\n');
 disp(B),disp(W);
 fprintf('Fano编码平均码字长度L：');
 disp(L);
 fprintf('信源的信息熵H：');
 disp(H);
 fprintf('Fano编码的编码效率q：');
 disp(q);

样例

clc；clear；
s=[1 2 3 4];
p=[0.10 0.32 0.18 0.10 0.07 0.06 0.15 0.12];%给定信源
fano(p);

信源符号出现的概率为:
  列 1 至 4

    0.1000    0.3200    0.1800    0.1000

  列 5 至 8

    0.0700    0.0600    0.1500    0.1200

Fano编码所得码字W：
  列 1 至 5

    {[1]}    {[2]}    {[3]}    {[4]}    {[5]}

  列 6 至 8

    {[6]}    {[7]}    {[8]}

000    001    01     100    1010   1011   110    111    
Fano编码平均码字长度L：    3.2500

信源的信息熵H：    2.9254

Fano编码的编码效率q：    0.9001

1.2 香农编码

香农（Shannon）在1948年首次提出“有扰信道编码定理”，讲述了怎样在有干扰的信道中进行可靠通信，该定理成为了纠错码理论的基石。纠错码理论的发展从此拉开序幕，日渐得到了相关研究者的重视，也对通信的发展产生了重要的影响。

在香农提出信道编码定理后，编码学者们就开始致力于构造逼近信道容量且实用的好码，一开始是对简单的线性分组码进行改造，得到了汉明码、格雷码、循环码和BCH码，但由于译码复杂度限制了码的长度，这些线性分组码的短码性能与香农极限的距离仍较远。后来，卷积码和级联码的出现为寻找近香农极限码提供了新的思路，如1993年提出的Turbo码。Turbo码是一种级联型卷积码，它在编码结构和译码算法上进行了巧妙的设计，具有突出的纠错性能，性能接近香农极限，且编译码的计算复杂度不高。

香农编码流程和费诺码相似，具体步骤为：

将信源符号按概率递减的方式进行排列 $p(a_{1})\geqslant p(a_{2})\geqslant p(a_{n})$ ；
按码长计算公式 $l=-log(p(a_{i}))$ 计算出每个信源符号的码长；
计算至每个符号的累加概率 $P_{i}=\sum_{1}^{i-1}p(a_{i})$ ；
取累加概率对应二进制数的小数点后 $l$ 位构成该信源符号的二进制码字

程序

clear all
clc
% 用户输入符号概率
p = input('请输入离散信源概率分布，例如[0.5,0.5]：\n');
N = length(p);
L = ceil(-log2(p));% 获得码长向量，元素表示每个符号所对应的码长

% 获得累加概率P及对应码字
[p_SortDescend,reflect] = sort(p,'descend');% 将概率从大到小进行排序
%注：reflect所表示的映射关系很重要
P = zeros(1,N);     % 初始化累加概率
CODE = strings(1,N); % 初始化对应码字(字符串形式)
for i=1:N           % i表示排序后第几个符号
    code = zeros(1,L(reflect(i)));% 初始化对应码字（数组形式）
    if i==1         % 定义第一个编码为0
        P(1)=0;
        CODE(reflect(i)) = num2str(code);
    else
        P(i) = sum(p_SortDescend(1,1:i-1)); % 获得累加概率
    end
    % 下面计算香农码（计算累加概率的二进制数，并取前Li位）
    p_count = P(i)*2;       % p_count用于逐步的计算累加概率的二进制数
    for m=1:L(reflect(i))   % m表示这个符号里第几个码字
       if p_count >= 1
            code(m) = 1;
            p_count = p_count-1;
       else
            code(m) = 0;
       end
       p_count = p_count*2;
    end
    % 将香农码赋值给对应的符号
    CODE(reflect(i)) = num2str(code);
end

H = sum(-p.*log2(p));   % 计算信源信息熵
L_ave = sum(L.*p);      % 计算平均码长
yita = H/L_ave;         % 计算编码效率

% 展示输出码字、平均码长和编码效率
fprintf('\n运行结果:\n');
disp(['信号符号: ',num2str(1:N)]);
disp(['对应概率: ',num2str(p)]);
fprintf('对应码字:');disp(CODE);
disp(['平均码长:',num2str(L_ave)]);
disp(['编码效率:',num2str(yita)]);

样例

请输入离散信源概率分布，例如[0.5,0.5]：
[0.10 0.32 0.18 0.10 0.07 0.06 0.15 0.12]

运行结果:
信号符号: 1  2  3  4  5  6  7  8
对应概率: 0.1        0.32        0.18         0.1        0.07        0.06        0.15        0.12
对应码字:    "1  1  0  0"    "0  0"    "0  1  0"    "1  1  0  1"    "1  1  1  1"    "1  1  1  1  1"    "1  0  0"    "1  0  1  0"

平均码长:3.49
编码效率:0.83823

1.3 哈夫曼编码

哈夫曼编码(Huffman Coding)，是可变字长编码(VLC)的一种。Huffman于1952年提出一种编码方法，该方法完全依据字符出现概率来构造异字头的平均长度最短的码字，有时称之为最佳编码，一般就叫做Huffman编码（有时也称为霍夫曼编码）。

哈夫曼编码的具体步骤如下：

将信源符号按概率递减的方式进行排列 $p(a_{1})\geqslant p(a_{2})\geqslant p(a_{n})$ ；
从队尾给m个最小的概率进行m元标号如二进制下将两个概率最小的符号标为0和1，并相加合并再次排序，重复此步骤直至概率和为1；
画出由概率１处到每个信源符号的路径，顺序记下沿路径的标号，所得就是该符号的霍夫曼码字。

程序

function huffman(P)

%输入：P信源概率矩阵
%返回h编码矩阵，l平均码长
N=length(P);
Q=P;
Index=zeros(N-1,N);  %初始化Index  
for i=1:N-1  
   [Q,L]=sort(Q);  %将P中的元素按升序排序后，元素放到Q中，对应的索引值存到L中
   Index(i,:)=[L(1:N-i+1),zeros(1,i-1)];
   G(i,:)=Q;%缩减信源得到的最终矩阵
   %Index为N-1行、N列矩阵，用来记录每行最小两概率叠加后概率排列次序,元素不足的地方补0
   %参考doc sort
   Q=[Q(1)+Q(2),Q(3:N),1]; %将Q中概率最小的两个元素合并，元素不足的地方补1
end
%根据以上建立的Index矩阵，进行回溯，获取信源编码
for i=1:N-1
    Char(i,:)=blanks(N*N);%初始化一个由空格符组成的字符矩阵N*N,用于存放编码
end
%从码树的树根向树叶回溯，即从G矩阵的最后一行按与Index中的索引位置的对应关系向其第一行进行编码
Char(N-1,N)='0';%G中的N-1行即最后一行第一个元素赋为0，存到Char中N-1行的N列位置
Char(N-1,2*N)='1';%G中的N-1行即最后一行第二个元素赋为1，存到Char中N-1行的2*N列位置
%以下从G的倒数第二行开始向前编码
for i=2:N-1  
   Char(N-i,1:N-1)=Char(N-i+1,N*(find(Index(N-i+1,:)==1)) -(N-2):N*(find(Index(N-i+1,:)==1))); 
   %将Index后一行中索引为1的编码码字填入到当前行的第一个编码位置
   Char(N-i,N)='0'; %然后在当前行的第一个编码位置末尾填入'0' 
   Char(N-i,N+1:2*N-1)=Char(N-i,1:N-1); %将G后一行中索引为1的编码码字填入到当前行的第二个编码位置 
   Char(N-i,2*N)='1';  %然后在当前行的第二个编码位置末尾填入'1'
   for j=1:i-1  
 %内循环作用：将Index后一行中索引不为1处的编码按照左右顺序填入当前行的
 %第3个位置开始的地方,最后计算到Index的首行为止
      Char(N-i,(j+1)*N+1:(j+2)*N)=Char(N-i+1,N*(find(Index(N-i+1,:)==j+1)-1)+1:N*find(Index(N-i+1,:)==j+1));  
 end  
end  

 %Char中第一行的编码结果就是所需的Huffman 编码输出，通过Index中第一行索引将编
 % 码对应到相应概率的信源符号上。
   for i=1:N  
      h(i,1:N)=Char(1,N*(find(Index(1,:)==i)-1)+1:find(Index(1,:)==i)*N);
      %根据Index第一行索引将Char中第一行编码值还原为输入概率矩阵中的顺序填入Result
      ll(i)=length(find(abs(h(i,:))~=32)); 
   end 
  l=sum(P.*ll);
H=0; 
for i=1:length(P)
    H=H-P(i)*log(P(i))/log(2);
end
eff=H/l;%编码效率
disp('信息源的平均信息量:');
disp(H);
disp('霍夫曼编码的平均码长:');
disp(l);
disp('所对应的编码:');
disp(h);
disp('编码效率:'),
disp(eff)

样例

 p=[0.10 0.32 0.18 0.10 0.07 0.06 0.15 0.12];
huffman(p)
信息源的平均信息量:
    2.9254

霍夫曼编码的平均码长:
    3.1100

所对应的编码:
     000
      10
     111
     001
    0111
    0110
     110
     010
编码效率:
    0.9407

1.4 算术编码

和哈夫曼编码一样，算数编码是熵编码的一种，也是基于数据中字符出现的概率，给不同字符以不同的编码。相较下哈夫曼编码所划分出来的子区间并不是严格按照概率的大小等比例划分的。哈夫曼编码可以看作是对算数编码的一种近似，它并不是完美地呈现原始数据中字符的概率分布。也正是因为这一点微小的偏差，使得哈夫曼编码的压缩率通常比算数编码略低一些。或者说，算数编码能更逼近香农给出的理论熵值。

简单来说，算数编码做了这样一件事情：

统计需要编码的序列里面所有的字符和出现的次数。
将区间 [0,1) 连续划分成多个子区间，每个子区间代表一个上述字符，区间的大小正比于这个字符在文中出现的概率 p。概率越大，则区间越大。所有的子区间加起来正好是 [0,1)。
编码从一个初始区间 [0,1) 开始，设置: $low=0,hight=1$
不断读入原始数据的字符，找到这个字符所在的区间，比如 [ L, H )，更新:

$low=low+(hight-low)*L$

$hight=low+(hight-low)*H$

5.最后将得到的区间 [low, high)中任意一个小数以二进制形式输出即得到编码的数据。

程序

function acode = suanshubianma(symbol, ps, inseq)
% 输入参数：
%   symbol: 符号数组，包含所有可能的符号
%   ps: 符号对应的概率数组
%   inseq: 待编码的输入符号序列
% 输出参数：
%   acode: 编码后的值
    % 计算符号的高区间边界
    high_range = [];
    for k = 1:length(ps)
        high_range = [high_range, sum(ps(1:k))];
    end

    % 计算符号的低区间边界
    low_range = [0, high_range(1:length(ps)-1)];

    % 将输入符号序列转换为符号对应的索引
    sbidx = zeros(size(inseq));
    for i = 1:length(inseq)
        sbidx(i) = find(symbol == inseq(i));
    end

    % 初始化编码区间的上下边界
    low = 0;
    high = 1;

    % 进行算术编码
    for i = 1:length(inseq)
        range = high - low;
        high = low + range * high_range(sbidx(i));
        low = low + range * low_range(sbidx(i));
    end

    % 编码结果为编码区间的下边界
    acode = low;
end

样例

symbol=['abcd'];
ps=[0.5 0.25 0.125 0.125];
inseq=('abda');
codeword=suanshubianma(symbol,ps,inseq)

codeword =

     3.593750000000000e-01

1.5 LZ78编码

1977年两位以色列研究者J.Ziv和A.Lempel独辟蹊径，完全脱离Huffman及算术编码的设计思路，创造出一系列比Huffman编码更有效，比算术编码更快捷的通用压缩算法。将这些算法统称为LZ系列算法。Ziv和Lempel于1977年提出了LZ77算法，1978年他们提出改进算法，即LZ78算法。1984年T.A.Welch提出了LZ77算法一个变种，即LZW算法。

LZ78词典编码实现文本数据压缩的基本原理是，在编码过程中一边编码一边自适应建立词典,为新字符串(即短语)建立索引。将当前短语用已出现的与之相同的短语索引号代替，从而压缩冗余，若不存在匹配串，则新建索引。

LZ系列算法的杰出性能使它在数据压缩领域获得广泛应用，UNIX系统最先出现了使用LZW算法的compress程序，该程序很快成为了UNIX世界的压缩标准。目前LZ77、LZ78、LZW算法以及他们的各种变体几乎垄断了整个通用数据压缩领域。PKZip、WinRAR、gzip等压缩工具以及ZIP、GIF、PNG等文件格式都是LZ系列算法的受益者。

程序

字典生成函数

%author:TimDing
%生成字典
function [D,m]=make_dictionary(origin,U,N)
%dictionary为中间变量，其第一行为字典的短语，
%第二行为段号，第三行为短语的长度，第四行为码在U中的编号，
%第五行为段号的二进制表示，第六行为码的二进制表示。
ll=length(origin);%源文件的长度
if (ll==1)%当只有一个字符时的处理
    D(1,1)={U};
    D(2,1)={'0'};
    m=1;
    return;
end
flag=1;%用于标记当前字符串是否与字典中的字符相等，1为相等，0为不相等
string='';%表示当前字符串
dictionary(1,1)={char(origin(1,1))};%设置第一个的字符
m=2;%字典长度

[U,N]=Num(origin);%求得信源的字符序列

N_bit=ceil(log2(N));%表示码编号的位数

for i=2:1:ll   
   string=char([string,origin(1,i)]);%将当前字符与下一个字符组合在一起   
    for j=1:1:m-1%进行字符串匹配
        if(strcmp(dictionary(1,j),string))
         flag=1;
         dictionary(2,m)={j};%记下段号
         break;
        else          
            flag=0;
        end
    end
        if(flag==0)
            k=length(string);
        string=reshape(string,1,k);%将string变为行向量
        dictionary(1,m)={string};
        dictionary(3,m)={k};%记下短语的长度
        string='';
        m=m+1;       
        end
     
    
end

if(~isempty(string))%处理最后面的特殊字符
    dictionary(1,m)={string};
    dictionary(3,m)={length(string)};
else
    m=m-1;
end
dictionary(3,1)={1};%补上第一个的长度


para_bit=ceil(log2(m-1));%向上取整得到段号所需的码长

for i=1:1:m
   if(isempty(dictionary{2,i}))
       dictionary(2,i)={0};%置零
   end
end

for i=1:1:m
    k=length(dictionary{1,i});
    for j=1:1:N
    if(dictionary{1,i}(k)==U(j))
        dictionary(4,i)={j};%得到在原文件中的次序
    end
    end
end




for i=2:1:m
 dictionary(5,i)={dec2bin(dictionary{2,i}(1),para_bit)};%用二进制表示段号

  dictionary(6,i)={dec2bin(dictionary{4,i}(1)-1,N_bit)};%用二进制表示码编号  
   
end
dictionary(5,1)={dec2bin(0,para_bit)};%补上第一个的段号
dictionary(6,1)={dec2bin(0,N_bit)};%用二进制表示码编号 
for i=1:1:m
    D(1,i)=dictionary(1,i);%最终的字典
    D(2,i)={[dictionary{5,i} dictionary{6,i}]};%最终的编码
end
end

样例

input='abbaabbabbabab ';
[dictionary,m]=make_dictionary(input);
code=encode(dictionary,m);
disp('编码为：');
disp(code);

编码为：
00000000010100000101011011010001010

二、连续信源编码

根据限失真信源编码定理，由于连续信源在时间和取值上都是连续的，因此要使得信源在时间和取值上离散，转变为离散信源，就可以采用离散编码定理进行编码。首先使其在时间上离散，采取的方法是抽样，抽取离散点，这些离散点必须可以代表整个信号。然后再经过均匀量化或非均匀量化使信号在取值上也离散，成为离散的数字信号。在值域上选取有限个量化值中的一个来代替信号值即量化，量化肯定带来误差。量化后转化为离散信源编码，而离散信源编码属无失真编码，因此连续信源编码的误差来自量化过程,下面简要介绍两类量化方法。

2.1均匀量化

均匀量化又称线性量化，在整个量化范围内的量化间隔都是相等的，其中当信号的量化间隔为 $\Delta$ 时，码长为

$k=log_{2}\tfrac{1}{\Delta }+1$

均匀量化分为平量化和升量化，主要以有无量化值来区分。均匀量化编码可以大致分为两个步骤，即极性判断和信号绝对值量化。其一般编码过程为：判断信号值的极性，确定极性码。将信号绝对值与量化码各位权值组合的逐次比较，确定量化码。把极性码与量化码组合起来，得到均匀量化码。

2.2非均匀量化

非均匀量化即量化范围内的量化间隔不完全相等。以1３折线Ａ律为例，了解非均匀编码的主要概念和方法。其非均匀量化编码过程为：判断信号值的极性，确定极性码；中位搜素段落码起始量化值，确定段落码；计算信号绝对值与所确定段落的起始量化值之差，然后逐次比较其与段内码各位权值组合的，确定段内码。组合起来即得到1３折线Ａ律非线性量化编码。由于各个段落的宽度不同，每个段落的内段内码各位的权值也不同。