实验四——DPCM编码

一.实验原理

    DPCM是差分预测编码调制的缩写，是比较典型的系统。在DPCM系统中，需要注意的是预测器的输入是已经解码以后的样本。之所以不用原始值来做预测，是因为在解码端无法得到原始样本，只能用存在误差的解码之后的数据，因此DPCM编码器中实际内嵌了一个解码器。
    在一个 DPCM系统中，有两个因素需要设计：预测器 和量化。理想情况下系统中，预测器和量化采取联合优化的方式。而在实际的系统中采用的是次优方式：分别进行线性预测器和量化器的优化设计。

    信源的相邻符号通常是相关的，预测编码就是利用信源相邻符合之间的相关性来进行压缩，根据某一模型利用以往的样本值对新样本进行预测，然后将样本的实际值与其预测值相减得到一个误差值，最后对这一误差值进行编码，如果模型足够好，且样本序列在时间上相关性较强，则误差信号的幅度将远远小于原始信号，从而得到较大的数据压缩。

二.实验框图

虚线方框中的内容就为编码器中内嵌的解码器。

三.DPCM+量化使信号产生的三种失真

（1）粒状噪声（Granular noise）：平坦区域的随机噪声

对于小信号的量化区间过于大了，减小小信号的量化区间可以解决此类的问题，一般采用的是非均匀量化或者是压缩扩张法，即所谓的“小信号细量化，大信号粗量化”。

（2）Edge busyness：抖动的边缘

平坦区向斜率过载区过渡产生的信号失真。
（3）斜率过载（Slope overload）：边缘模糊

信号的量化区间过于小了，需要增大信号的量化区间。

四.代码分析

main主函数

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
#include "dpcm.h"
#include <windows.h>


#define u_int8_t	unsigned __int8
#define u_int		unsigned __int32
#define u_int32_t	unsigned __int32

#pragma pack(1)
	BITMAPFILEHEADER File_header;
	BITMAPINFOHEADER Info_header;
	long frameWidth ;			/* --width=<uint> */
	long frameHeight ;		/* --height=<uint> */

unsigned char rectify(unsigned char x)
	{
		if (x<0) x=0;
		if (x>255) x=255;
		return x;
	}

int BMP2YUV(char *url,char *url_out1,char *url_out2,int num)
{
	
	FILE *bmpFile=NULL;
	FILE *pFile=NULL;//存储重建图像值
	FILE *eFile=NULL;//存储差值
	bool flip = TRUE;/* --flip */
	unsigned int i;
	unsigned char* bmpBuf = NULL;
	unsigned char* pBuf = NULL;//存储重建图像值
	//int eBuf;
	unsigned char* eBuf_p = NULL;//存储差值
	unsigned char* yBuf = NULL;//仅对亮度进行预测编码
	unsigned char* uBuf = NULL;
	unsigned char* vBuf = NULL;
	u_int32_t videoFramesWritten = 0;
	bmpFile=fopen(url,"rb+");
	pFile=fopen(url_out1,"wb+");
	eFile=fopen(url_out2,"wb+");
	//判断文件头是否读取正确
	  if(fread(&File_header,sizeof(BITMAPFILEHEADER),1,bmpFile) != 1)
	{
		printf("read file header error!");
		exit(0);
	}
	  //判断是否为bmp文件
	if (File_header.bfType != 0x4D42)
	{
		printf("Not bmp file!");
		exit(0);
	}
	else
	{	printf("this is a bmp\n");
	}
	//判断信息头是否读取正确
if(fread(&Info_header,sizeof(BITMAPINFOHEADER),1,bmpFile) != 1)
	{	
		printf("read info header error!");
		exit(0);
	}
	//	end read header

	//判断长宽高是否读取正确，并赋值
	if (((Info_header.biWidth/8*Info_header.biBitCount)%4) == 0)
		frameWidth= Info_header.biWidth;
	else
		frameWidth = (Info_header.biWidth*Info_header.biBitCount+31)/32*4;
	if ((Info_header.biHeight%2) == 0)
		frameHeight = Info_header.biHeight;
	else
		frameHeight = Info_header.biHeight + 1;
	printf("The width is %d\n",frameWidth);
	printf("The height is %d\n",frameHeight);
	/* open the BMP file */
	if (bmpFile == NULL)
	{
		printf("cannot find bmp file\n");
		exit(1);
	}
	else
	{
		printf("The input bmp file is open\n"/*, bmpFileName*/);
	}

	/* open the pfile */
	if (pFile == NULL)
	{
		printf("cannot find p file\n");
		exit(1);
	}
	else
	{
		printf("The output pfile is open \n"/*, yuvFileName*/);
	}
	/* open the efile */
	if (eFile == NULL)
	{
		printf("cannot find e file\n");
		exit(1);
	}
	else
	{
		printf("The output efile is open \n");
	}
	
	bmpBuf = (unsigned char*)malloc(frameWidth * frameHeight * 3);
	memset(bmpBuf,0,frameHeight*frameWidth*3);
	/* get the output buffers for a frame */
	yBuf = (unsigned char*)malloc(frameWidth * frameHeight);
	//memset(yBuf,0,frameHeight*frameWidth);
	pBuf = (unsigned char*)malloc(frameWidth * frameHeight);
	eBuf_p =(unsigned char*)malloc(frameWidth * frameHeight);
	uBuf = (unsigned char*)malloc((frameWidth * frameHeight)/4 );
	vBuf = (unsigned char*)malloc((frameWidth * frameHeight) /4);
	/* get an input buffer for a frame */
	if (bmpBuf == NULL || yBuf == NULL || pBuf == NULL || eBuf_p == NULL)
	{
		printf("no enought memory\n");
		exit(1);
	}
	//
	void ReadRGB(FILE * pFile,BITMAPFILEHEADER & file_h,BITMAPINFOHEADER & info_h,unsigned char * rgbData);
		ReadRGB(bmpFile, File_header, Info_header,bmpBuf);
	

	//while (fread(bmpBuf, 1, frameWidth * frameHeight*3 , bmpFile)) 
	//{
		if(RGB2YUV(frameWidth, frameHeight, bmpBuf , yBuf, uBuf, vBuf, flip))
		{
			printf("error");
			return 0;
		}
		int bits=8;
		int scale=512/(1<<bits);
		printf("%d\n",scale);
		int eBuf ;//存储差值
		//一阶
		//for( int j=0;j<frameHeight;j++)//j：行
		//{
		//	eBuf=yBuf[j*frameWidth]-128;
		//	eBuf_p[j*frameWidth]=(eBuf+255)/2;
		//	yBuf[j*frameWidth]=128+(eBuf_p[j*frameWidth]-128)*2;
		//	pBuf[frameWidth*j] = (eBuf_p[frameWidth*j] - 128)*2 + 128;
		//	for(i=1+j*frameWidth;i<frameWidth*(j+1);i++)
		//	{//i:列
		//		eBuf=yBuf[i]-yBuf[i-1];
		//		eBuf_p[i]=(eBuf+255)/2;//预测误差的量化值
		//		yBuf[i] = (eBuf_p[i] - 128)*2 + yBuf[i-1];
		//		pBuf[i] = (eBuf_p[i] - 128)*2 + pBuf[i-1];	 
		//	}//
		//}
		
		//4bit量化
		//for( int j=0;j<frameHeight;j++)//j：行
		//{
		//	eBuf=yBuf[j*frameWidth]-128;
		//	eBuf_p[j*frameWidth]=(eBuf+255)/2;
		//	yBuf[j*frameWidth]=128+(eBuf_p[j*frameWidth]-127)*2;
		//	pBuf[frameWidth*j] = (eBuf_p[frameWidth*j] - 127)*2 + 128;
		//	for(i=1+j*frameWidth;i<frameWidth*(j+1);i++)
		//	{//i:列
		//		eBuf=yBuf[i]-yBuf[i-1];
		//		eBuf_p[i]=(eBuf+255)/2;//预测误差的量化值
		//		yBuf[i] =(eBuf_p[i] - 127)*2+ yBuf[i-1];
		//		pBuf[i] = (eBuf_p[i] - 127)*2 + pBuf[i-1];	 
		//	}//
		//}
		//二阶2
		for( int j=1;j<frameHeight;j++)//j：行
		{
			eBuf=yBuf[j*frameWidth]-128;
			eBuf_p[j*frameWidth]=(eBuf+255)/2;
			yBuf[j*frameWidth]=128+(eBuf_p[j*frameWidth]-255/2)*2;
			pBuf[frameWidth*j] = rectify((eBuf_p[frameWidth*j] - 255/2)*2 + 128);
			for(i=1;i<frameWidth;i++)
			{//i:列
				eBuf=yBuf[i+j*frameWidth]-0.5*yBuf[i+j*frameWidth-1]-0.25*yBuf[i+(j-1)*frameWidth]-0.25*yBuf[i+1+(j-1)*frameWidth];
				eBuf_p[i+j*frameWidth]=(eBuf+255)/2;//预测误差的量化值
				yBuf[i+j*frameWidth] = rectify((eBuf_p[i+j*frameWidth] - 255/2)*2 + 0.5*yBuf[i+j*frameWidth-1]+0.25*yBuf[i+(j-1)*frameWidth]+0.25*yBuf[i+1+(j-1)*frameWidth]);
				pBuf[i+j*frameWidth] = rectify((eBuf_p[i+j*frameWidth] - 255/2)*2 + 0.5*pBuf[i+j*frameWidth-1]+0.25*pBuf[i+(j-1)*frameWidth]+0.25*pBuf[i+1+(j-1)*frameWidth]);	 
			}//
		}
		//一阶0.95
		//for( int j=0;j<frameHeight;j++)//j：行
		//{
		//	eBuf=yBuf[j*frameWidth]-128;
		//	eBuf_p[j*frameWidth]=(eBuf+255)/2;
		//	yBuf[j*frameWidth]=128+(eBuf_p[j*frameWidth]-255/2)*2;
		//	pBuf[frameWidth*j] = (eBuf_p[frameWidth*j] - 255 /2)*2 + 128;
		//	for(i=1+j*frameWidth;i<frameWidth*(j+1);i++)
		//	{//i:列
		//		eBuf=yBuf[i]-0.95*yBuf[i-1];
		//		eBuf_p[i]=(eBuf+255)/2;//预测误差的量化值
		//		yBuf[i] = (eBuf_p[i] - 127)*2 + 0.95*yBuf[i-1];
		//		pBuf[i] = (eBuf_p[i] - 127)*2 + 0.95*pBuf[i-1];	 
		//	}//
		//}
		
		for (i = 0; i < frameWidth*frameHeight; i++)
		{
			if (pBuf[i] < 16) pBuf[i] = 16;     
			if (pBuf[i] > 235) pBuf[i] = 235;
		}

		for (i = 0; i < frameWidth*frameHeight/4; i++)
		{
			if (uBuf[i] < 16) uBuf[i] = 16;
			if (uBuf[i] > 240) uBuf[i] = 240;

			if (vBuf[i] < 16) vBuf[i] = 16;
			if (vBuf[i] > 240) vBuf[i] = 240;
		}
		fwrite(pBuf, 1, frameWidth * frameHeight, pFile);
		fwrite(uBuf, 1, (frameWidth * frameHeight)/4 , pFile);
		fwrite(vBuf, 1, (frameWidth * frameHeight)/4 , pFile);

		fwrite(eBuf_p, 1, frameWidth * frameHeight, eFile);
		fwrite(uBuf, 1, (frameWidth * frameHeight)/4 , eFile);
		fwrite(vBuf, 1, (frameWidth * frameHeight)/4 , eFile);


		printf("\r...%d", ++videoFramesWritten);
	//}
	printf("\n%u %ux%u video frames written\n", 
		videoFramesWritten, frameWidth, frameHeight);
	
	free(bmpBuf);
	free(yBuf);
	free(uBuf);
	free(vBuf);
	/* cleanup */
	fclose(bmpFile);
	fclose(pFile);
	fclose(eFile);
	return 0;
}
void ReadRGB(FILE * pFile,BITMAPFILEHEADER & file_h,BITMAPINFOHEADER & info_h,unsigned char * rgbData)
{
	bool MakePalette(FILE * pFile,BITMAPFILEHEADER &file_h,BITMAPINFOHEADER & info_h,RGBQUAD *pRGB_out);
	unsigned long Loop,iLoop,jLoop,width,height,w,h;
	unsigned char mask, *Index_Data,* Data;

	if (((info_h.biWidth/8*info_h.biBitCount)%4) == 0)
		w = info_h.biWidth;
	else
		w = (info_h.biWidth*info_h.biBitCount+31)/32*4;
	if ((info_h.biHeight%2) == 0)
		h = info_h.biHeight;
	else
		h = info_h.biHeight + 1;

	width = w/8*info_h.biBitCount;
	height = h;

	Index_Data = (unsigned char *)malloc(height*width);
	Data = (unsigned char *)malloc(height*width);

	fseek(pFile,file_h.bfOffBits,0);
	if(fread(Index_Data,height*width,1,pFile) != 1)
	{
		printf("read file error!\n\n");
		exit(0);
	}
	for ( iLoop = 0;iLoop < height;iLoop ++)
		for (jLoop = 0;jLoop < width;jLoop++)
		{
 			Data[iLoop*width+jLoop] = Index_Data[(height-iLoop-1)*width+jLoop];
		}

	switch(info_h.biBitCount)
	{
	case 24:
		memcpy(rgbData,Data,height*width);
		if(Index_Data)
			free(Index_Data);
		if(Data)
			free(Data);
		return;
	case 16:
		if(info_h.biCompression == BI_RGB)
		{
			
			for (Loop = 0;Loop < height * width;Loop +=2)
			{
				*rgbData = (Data[Loop]&0x1F)<<3;
				*(rgbData + 1) = ((Data[Loop]&0xE0)>>2) + ((Data[Loop+1]&0x03)<<6);
				*(rgbData + 2) = (Data[Loop+1]&0x7C)<<1;
				rgbData +=3;
			}
		}
		if(Index_Data)
			free(Index_Data);
		if(Data)
			free(Data);
		return;
		default:
		RGBQUAD *pRGB = (RGBQUAD *)malloc(sizeof(RGBQUAD)*(unsigned long)pow(float(2),info_h.biBitCount));
		int temp = sizeof(pRGB);
		if(!MakePalette(pFile,file_h,info_h,pRGB))
			printf("No palette!\n\n");

		for (Loop=0;Loop<height*width;Loop++)
		{
		
			switch(info_h.biBitCount)
			{
			case 1:
				mask = 0x80;
				break;
			case 2:
				mask = 0xC0;
				break;
			case 4:
				mask = 0xF0;
				break;
			case 8:
				mask = 0xFF;
			}

			int shiftCnt = 1;

			while (mask)
			{
				unsigned char index = mask == 0xFF ? Data[Loop] : ((Data[Loop] & mask)>>(8 - shiftCnt * info_h.biBitCount));
				* rgbData = pRGB[index].rgbBlue;
				* (rgbData+1) = pRGB[index].rgbGreen;
				* (rgbData+2) = pRGB[index].rgbRed;

				if(info_h.biBitCount == 8)
					mask =0;
				else
					mask >>= info_h.biBitCount;
				rgbData+=3;
				shiftCnt ++;
			}
		}
		if(Index_Data)
			free(Index_Data);
		if(Data)
			free(Data);
// 		if(pRGB)
// 			free(pRGB);
	}

}

bool MakePalette(FILE * pFile,BITMAPFILEHEADER &file_h,BITMAPINFOHEADER & info_h,RGBQUAD *pRGB_out)
{
	if ((file_h.bfOffBits - sizeof(BITMAPFILEHEADER) - info_h.biSize) == sizeof(RGBQUAD)*pow(float(2),info_h.biBitCount))
	{
		fseek(pFile,sizeof(BITMAPFILEHEADER)+info_h.biSize,0);
		fread(pRGB_out,sizeof(RGBQUAD),(unsigned int)pow(float(2),info_h.biBitCount),pFile);
		return true;
	}
	else
		return false;
// 			free(pRGB);
	}

int main(int argc, char **argv)
{	
	int count=1;
	u_int num;
	FILE *bmpFile=NULL;
	FILE *pFile=NULL;
	FILE *eFile=NULL;
	char *bmpFileName=NULL;
	char *pFileName=NULL;
	char *eFileName=NULL;
	unsigned char* pBuf = NULL;
	unsigned char* eBuf = NULL;
	bmpFileName=argv[1];
	pFileName=argv[2];
	eFileName=argv[3];
	num=1;
	bmpFile=fopen(bmpFileName,"rb+");
	pFile=fopen(pFileName,"wb+");
	eFile=fopen(eFileName,"wb+");

	BMP2YUV(bmpFileName,pFileName,eFileName,1);
	//pBuf = (unsigned char*)malloc((frameWidth * frameHeight)*3/2);
	return 0;
}

DCMP.cpp

int RGB2YUV (int x_dim, int y_dim, void *bmp, void *y_out, void *u_out, void *v_out, int flip)
{
	static int init_done = 0;

	long i, j, size;
	unsigned char *r, *g, *b;
	unsigned char *y, *u, *v;
	unsigned char *pu1, *pu2, *pv1, *pv2, *psu, *psv;
	unsigned char *y_buffer, *u_buffer, *v_buffer;
	unsigned char *sub_u_buf, *sub_v_buf;

	if (init_done == 0)
	{
		InitLookupTable();//程序一开始，便运行一次亮度计算公式的组件
		init_done = 1;
	}

	// check to see if x_dim and y_dim are divisible by 2，查看是否可以被2整除
	if ((x_dim % 2) || (y_dim % 2)) return 1;//如果不能被2 整除，则return1
	size = x_dim * y_dim;

	// allocate memory
	y_buffer = (unsigned char *)y_out;
	sub_u_buf = (unsigned char *)u_out;
	sub_v_buf = (unsigned char *)v_out;
	u_buffer = (unsigned char *)malloc(size * sizeof(unsigned char));//分配指定字节的内存空间
	v_buffer = (unsigned char *)malloc(size * sizeof(unsigned char));
	if (!(u_buffer && v_buffer))
	{
		if (u_buffer) free(u_buffer);//如果其中任何一个没有分配内存空间，就释放
		if (v_buffer) free(v_buffer);
		return 2;
	}

	b = (unsigned char *)bmp;//b指向bmp文件的首地址的指针
	y = y_buffer;
	u = u_buffer;
	v = v_buffer;

	// convert RGB to YUV
	if (!flip) {
		for (j = 0; j < y_dim; j ++)
		{
			y = y_buffer + (y_dim - j - 1) * x_dim;
			u = u_buffer + (y_dim - j - 1) * x_dim;
			v = v_buffer + (y_dim - j - 1) * x_dim;

			for (i = 0; i < x_dim; i ++) {
				g = b + 1;
				r = b + 2;
				*y = (unsigned char)(  RGBYUV02990[*r] + RGBYUV05870[*g] + RGBYUV01140[*b]);
				*u = (unsigned char)(- RGBYUV01684[*r] - RGBYUV03316[*g] + (*b)/2          + 128);
				*v = (unsigned char)(  (*r)/2          - RGBYUV04187[*g] - RGBYUV00813[*b] + 128);
				b += 3;
				y ++;
				u ++;
				v ++;
			}
		}//倒序
	} else {
		for (i = 0; i < size; i++)
		{
			g = b + 1;
			r = b + 2;
			*y = (unsigned char)(  RGBYUV02990[*r] + RGBYUV05870[*g] + RGBYUV01140[*b]);//强制类型转换会产生失真
			*u = (unsigned char)(- RGBYUV01684[*r] - RGBYUV03316[*g] + (*b)/2          + 128);//调用函数按照亮度计算公式计算yuv
			*v = (unsigned char)(  (*r)/2          - RGBYUV04187[*g] - RGBYUV00813[*b] + 128);
			b += 3;
			y ++;
			u ++;
			v ++;
		}
	}

 原始图像                                                                         误差图像                                                                       重建图像

第一幅图片的原始概率分布与误差概率分布

和第二幅图的原始概率分布与误差分布（所有图片的误差概率分布形状近似，相对集中）