一文带你看懂算术编码(C语言)

算术编码C语言

简介

算术编码是图像压缩的主要算法之一。 是一种无损数据压缩方法，也是一种熵编码的方法。和其它熵编码方法不同的地方在于，其他的熵编码方法通常是把输入的消息分割为符号，然后对每个符号进行编码，而算术编码是直接把整个输入的消息编码为一个数，一个满足(0.0 ≤ n < 1.0)的小数n(百度百科)(你伤害辽我，还一笑而过~)

原理

算术编码的基本原理是：根据信源的不同符号序列的概率，把[0，1]区间划分为互不重叠的子区间，子区间的宽度恰好是各符号序列的概率。这样信源发出的不同符号序列将与各子区间一一对应，每个子区间内的任意一个实数都可以用来表示对应的符号序列，这个二进制小数就是该符号序列所对应的码字，截取其长度与该序列的概率匹配，实现高效编码。
接下来我们看具体的步骤：
编码
1.假设对一段二进制a1,a2符号序列信源进行编码，先统计出a1,a2出现的概率。
2.计算出该序列的概率p=p(x1)…p(x2),然后根据下面的公式计算出截止长度$$L=\lceil lo{g}_2\frac{1}{p(S)}\rceil$$其中方框表示向上取整。
3.从左到右根据信源符号来截取对应的区间，如第一个符号为a1，则截取0到1区间的a1段，第二个符号为a2，则在上一步的基础上截取其中的a2段，即在0到a1区间取a2段，以此类推直到序列的最后一个符号。
4.在3计算出的序列区间内任取一个数作为码字，截取其二进制的小数点后前L位，即得信源序列的算数编码。

解码
根据算术编码的原理我们可以倒推出解码得步骤：
1.将编码转回十进制小数（主要由于编码时截取的操作使得有截止误差的存在，所以这里转回十进制最好加上一个0.5*2^(-L+1)的补偿量）
2.判断该小数落在哪个区间，如其落在第一个大区间内则可判断源信源序列的第一个符号为a1，又落在第一个大区间内的第二个小区间，则信源序列的第二个符号为a2，以此逐渐细化区间得到全部信源序列信源。

这里结合一个例题进行详细步骤说明：
例：设二进制无记忆信源S={0,1}，其p(0)=1/4,p(1)=3/4。对二元序列11110做算术编码解码。
解：
$p(S)=p(1{)}^6\ast p(0{)}^2=(\frac{3}{4}{)}^4\ast \frac{1}{4}$
$L=\lceil lo{g}_2\frac{1}{p(S)}\rceil =4$

因为x1=1，所以首先落在[$\frac{1}{4}$,1]区间，如图所示：

x2=1，所以落在$low=\frac{1}{4}+(1-\frac{1}{4})\ast \frac{1}{4}=\frac{7}{16}$,即[$\frac{7}{16}$,1]区间，如图所示。
x3=1,所以落在$low=\frac{7}{16}+(1-\frac{7}{16})\ast \frac{1}{4}=\frac{37}{64}$,即[$\frac{37}{64}$,1]区间。

x4=1,所以落在$low=\frac{37}{64}+(1-\frac{37}{64})\ast \frac{1}{4}=\frac{175}{256}$,即[$\frac{175}{256}$,1]区间。

x5=0,所以落在$high=\frac{175}{256}+(1-\frac{175}{256})\ast \frac{1}{4}=\frac{781}{1024}$,即[$\frac{175}{256}$,$\frac{781}{1024}$]区间。
$0.5\ast (\frac{175}{256}+\frac{781}{1024})$=0.7231445=0.10111…
取小数点后前4位得算术编码结果为1011。

解码：
$p(S)=1\ast 2^{-1}+0\ast 2^{-2}+1\ast 2^{-3}+1\ast 2^{-4}+0.5\ast 2^{-5}=0.701325$(最后一项为了补偿截止误差)
因为p(S)>1/4=0.25,
即在[0.25,1]，所以第一个符号为1，此时low=0.25,high=1.

又p(S)>low+(high-low)*0.25=0.4375,
即在[0.4375,1]区间,所以第二个符号为1,low=0.4375,high=1.

又p(S)>low+(high-low)*0.25=0.578125,
即在[0.578125,1]区间,所以第三个符号为1,low=0.578125,high=1.

又p(S)>low+(high-low)*0.25=0.68359375,
即在[0.68359375,1]区间,所以第四个符号为1,low=0.68359375,high=1.

又p(S)<low+(high-low)*0.25=0.7627…,
即在[0.68359375,0.7627]区间,所以第五个符号为0.

解码完成，得解码为11110，与源码一致。

从上面的过程中我们可以总结出一个公式，即low+(high-low)*p(0)，上一轮的low,high中的一个被该值取代而与另外一个构成新的区间。

C语言实现代码

代码编程思路与面的例题思路基本相同，这里不再过多赘述。在同目录下新建source.txt,compress.txt和decode.txt，即可直接运行。

注意：1.由于C语言的变量容量限制，因此当信源长度过长时需要对其进行分段编码
2.由于C语言中double变量精度的限制，所以分段编码的长度不能过长，太长了后面的码就会变成错码，亲测每段<55为宜。

# include <stdio.h>
# include <stdlib.h>
# include <time.h>
# include <math.h>
clock_t start,stop;
double duration;

void codes(float p,int L,int n)
{
	int i,j,len,N,temp; 
	double ps,PS,high,low,m;
	char source[L/n+1];    //source用来存放伯努利源 
	char code[100];        //code用来存放编码 
	start = clock();
	FILE *sc;
	FILE *file;
	file = fopen("compress.txt","w");   
	sc = fopen("source.txt","r");  
	if(file==NULL|sc==NULL)
	{
		printf("打开文件错误\n");
		exit(0);
	}
	for(i=0;i<n;i++)
	{
		//算出截止长度N 
		fgets(source,100,sc);
    //	printf("source:%s\n",source);  	
		ps=1;
		for(j=0;j<L/n;j++)
		{
			if(source[j]=='0')
				ps = ps*p;
			else
				ps = ps*(1-p);
		}
    	PS = log(1.0/ps)/log(2);
    	N = (int)PS;
    	if(PS-N>1e-6)
    		N = N+1;
   	//	printf("N的值:%d\n",N);
   	
   		//算出区间的两端端点值 
    	low=0;
		high=1;
		for(j=0;j<L/n;j++)
		{
			if(source[j]=='0')
			high=(high-low)*p+low;
			else
			low = (high-low)*p+low;			
		}
	
		m = (low+high)*0.5;
		//printf("\nm的值为 %lf",m);
		
		//将m转为二进制取前N位 
		len=0;
		while(m)
		{
			temp = (int)(m*2);
			if(temp)
				code[len]='1';
			else
				code[len]='0';
			fputc(code[len],file);
			len=len+1;
			m = 2*m-temp;	
			if(len==N) 
				break;       
		}
		fputc(10,file);

		/*
		printf("\n编码结果为:");
		for(j=0;j<len[i];j++)
		{
			printf("%c",code[i][j]);
		}
		printf("\n");
		*/
	}
	fclose(sc);
	fclose(file);
	stop = clock();
	duration=  (double)(stop-start)/CLK_TCK;
	printf("\n编码共计耗时：%lf\n",duration);
	
}

void decodes(float p,int L,int n)
{
	int i,j,f,len,sum;
	double P,low,high;

	char code[100];        //code用来存放编码 
	char decode[L/n+1];    //decode用来存放解码 
	start = clock();
	FILE *compress;
	FILE *file;
	compress = fopen("compress.txt","r");  
	file = fopen("decode.txt","w");   
	if(file==NULL|compress==NULL)
	{
		printf("打开文件错误\n");
		exit(0);
	}

	for(i=0;i<n;i++)
	{
		fgets(code,100,compress);   //获取第一行数据 
		for(f=0;f<100;f++)
		{
			if(code[f]==10)        //获取第一行长度  10为换行符 asc码
				break;
		}
		len=f;	
		//printf("len:%d\n",len);
		
		//将编码转回十进制小数 
		P=0;
		for(j=0;j<len;j++)
		{
			P = P+(code[j]-48)*pow(2,-(j+1));		
		} 
		P = P+pow(2,-(len+1))*0.5;
	//	printf(" P  = %lf\n",P);
	
		//判断转出的十进制小数在哪个区间进行解码并写入decode文件 
		low=0;
		high=1;
		for(j=0;j<L/n;j++)
		{
			if((P>low)&(P<((high-low)*p+low)))
			{
				decode[j]='0';
				high = (high-low)*p+low;
			}
			else
			{
				decode[j]='1';
				low = (high-low)*p+low;
			}
			fputc(decode[j],file);		
		}
		fputc(10,file);
	/*
		printf("解码结果为:");
		for(j=0;j<L/n;j++)
		{
			printf("%c",decode[i][j]);	
		}
		printf("\n");
*/
		sum = sum+len;	
	} 
	fclose(file);
	fclose(compress);
	
	stop = clock();
	duration=  (double)(stop-start)/CLK_TCK;
	printf("\n解码共计耗时：%lf\n",duration);
	printf("\n压缩率为:%f \n",sum/(L*1.0));
}

int main()
{

	float p=0;	
	int i,j,k,L=0,n=1;
	
	while(1)
	{
		printf("输入概率p:");
		scanf("%f",&p);
		printf("输入长度L：");
		scanf("%d",&L);              
		printf("输入分组数n：");
		scanf("%d",&n);
		if (L%n>1e-6|p>=1|p<=0)
			{
				printf("输入错误，要求0<p<1,L为n的整数倍，请重新输入\n");
				continue;				
			}
		else
	 		break;
	}

	
       
	char source[L/n+1];    //source用来存放伯努利源 
	
	FILE *file;
	file = fopen("source.txt","w");   
	if(file==NULL)
	{
		printf("打开文件错误\n");
		exit(0);
	}
	
	//生成伯努利源 
	for(i=0;i<n;i++)
	{
		for(j=0;j<L/n;j++)
		{
			k = rand()%100;
			if(k>100*p-1){
				source[j]='1';
				}
			else
			{
				source[j]='0';			
			}
			fputc(source[j],file);		
		}
		fputc(10,file);
	}
	fclose(file);
		
	//编码 
	codes(p,L,n);	
	
	//解码
	decodes(p,L,n);
	
	return 0;
}

参考资料:信息论与编码（曹雪虹著）

第一次认认真真写博客，有很多不足欢迎评论指正，另外还请多多点赞支持啊啊啊！