如何用反余弦求pi c语言,反余弦变换 IDCT 8X8 的快速算法和SIMD优化（二,C程序优化）...

最新推荐文章于 2021-09-22 21:03:37 发布

weixin_39791152

最新推荐文章于 2021-09-22 21:03:37 发布

阅读量494

点赞数

文章标签： IDCT 快速算法 SIMD优化 DCT系数快速通道

反余弦变换 IDCT 8X8

的快速算法和SIMD优化(二,C程序优化)

8X8 IDCT 程序

前面我们学习了Feig等人对快速IDCT算法的论述。从算法的角度来说，2D-IDCT显然是最快的。但在在实际应用中，二维IDCT由于实现复杂，不常使用，常见的是一维IDCT。这里举出了几种IDCT的实现程序。包括笔者实现的AAN的几个变种，常见的CHENG_WANG，还有著名的JREV，以及几个MMX

IDCT。为了便于测试，我们设定了IDCT函数的输入与输出两个地址参数。实际使用中往往是IN-PLACE

，输入输出地址是同一个。这对性能的影响很微弱，并且修改成一个地址的工作也很容易。

AAN已经是很快的了。但是我们当然想追求更具效率的实现。出发点是，DCT系数中有很多的点数值为零，我们具体分析一下如何利用这个特性。

对于一维IDCT

AAN，如果存在二个以上系数不为零，计算量差别很小，因此我们只考虑一到二个系数非零的情况。我们姑且把这种特殊情况称之为快速通道。

为方便起见，我们称为零系数较多的情况为高负载，反之为低负载。

1)整个8x8矩阵，只有一个系数f(x，y)不为零。

此时一维IDCT： S8(u) =

Cu/2 * ∑ f(x)cos ((2x+1)*pi*u / 16 )，

退化为：S8(u) = Cu/2 * f(x)cos ((2x+1)*pi*u

/ 16 ),这说明每个一维IDCT通过8次乘法即可完成。根据x位置不同，这个表格是8 x 8 的。

a) f(x，y)在第0行第0列，两个一维IDCT可以直接合并，S8(u) S8(v) = f(0，0)/ 8;只需要64次移位。

b) x > 0,y = 0,

f(x，y)在第0行第x列。首先在0行作8次乘法，每个系数乘一个2√2；然后做8个列变换，S8(u) = f(0)/

8；8次移位，8次乘法。

c) y > 0,x = 0,

f(x，y)在第y行第0列。首先在0列作8次乘法，每个系数乘一个2√2；然后做8个行变换，S8(u) = f(0)/

8；也是8次移位，8次乘法。

d) x,y > 0,

f(x，y)在第y行第x列。首先在y行作8次乘法；然后做8个列变换，每个列变换也作8次乘法。一共是72次乘法。

当x不在0行0列的时候，还有另外的办法。

S8(u) S8(v) =

Cu/2 * Cv/2 *∑ ∑f(x,y)cos ((2x+1)*pi*u

/ 16 ) * cos ((2y+1)*pi*v / 16 )

退化为：

S8(u) S8(v) =

Cu/2 * Cv/2 *f(x,y)cos ((2x+1)*pi*u /

16 ) * cos ((2y+1)*pi*v / 16 )

这相当于把前面8 x 8 的表展开，使用一个64 x

64的超级大表，直接作64个乘法来搞定。具体可以参看JREV

IDCT的实现文件。但是过大的表格总是让人心有余悸的，不用也罢。使用两场一维IDCT，表格很小，速度更快。

我们姑且把f(0，0)不为零的情况称为快速通道0，其他的情况称为快速通道1。

2) 整个8x8矩阵，多于一个的系数不为零，但都在一行或者一列。

都在一行的时候，首先在该行作标准一维行变换，然后做8个退化的一维列变换。

都在一列的时候，首先在该列作标准一维列变换，然后做8个退化的一维行变换。

3) 整个8x8矩阵，多于一个的系数不为零，但都在二行或者二列。

当有2个系数不为零的时候，每个1D-IDCT需要做16次乘法，但是因为没有其他的操作，还是快一些的。

都在二行的时候，首先作2行标准一维行变换，然后做8个退化的一维列变换。

都在二列的时候，首先作2列标准一维列变换，然后做8个退化的一维行变换。

我们把2，3这种情况称为快速通道2。

单个行或者列，只有一个系数不为零。在该行或者该列作退化的一维变换。DC系数非零称为快速通道3，AC系数非零称为通道4。由于在1D-IDCT

里面，进行条件判断和程序跳转的代价很大，2个系数不为零的情况，如果不采用一些办法降低程序跳转的性能损失，２点的情况不会提高速度。

快速通道的总结：

快速通道代号

存在条件

8x8矩阵，只有f(0,0) ！= 0

8x8矩阵，只有f(x,y) != 0, x > 0, y >

8x8矩阵，所有f(x,y) != 0 都在一行或一列，或者二行或二列上

1x8或者8x1矩阵，只有f(0) ！= 0

1x8或者8x1矩阵，f(x) ！= 0(ｘ > 0)。

或者f(x１) ！= 0　和　f(x２) ！= 0。

我们根据以上分析来实现IDCT程序。限于篇幅，这里只给出AAN_EX的代码，其余代码在测试工程中都可以找到。

首先来看看我们需要获取的DCT系数位置信息。

1) 总的非零系数的数量。

2) 当只有一个非零系数的时候，该系数的位置。

3) 总的非零行的数量。非零行的位置。

4) 总的非零列的数量。非零列的位置。

5) 每行非零系数的数量和位置。

看起来记录这些数据很繁琐。但是要看到，所有这些信息都可以查表快速获取。

// LxIDCT.h: interface for the LxIDCT

#define PI 3.14159265359

#define

IDCT_TAB_RC 12

#define

IDCT_TAB_RC_DELTA 4

#define

IDCT_TAB_RC0 (IDCT_TAB_RC - IDCT_TAB_RC_DELTA)

#define

IDCT_TAB_ROUND0 (1<

#define

IDCT_TAB_RC1 (IDCT_TAB_RC + IDCT_TAB_RC_DELTA + 1)

#define

IDCT_TAB_ROUND1 (1<

#define

IDCT_TAB_DC_RC0 (IDCT_TAB_RC_DELTA -

1) // 4;

#define

IDCT_TAB_DC_RC1 (IDCT_TAB_RC_DELTA +

2) // 6;

#define

IDCT_TAB_DC_ROUND1 (1<

#define SET_NZ_MASK(A) (A)

#define SET_NZ_INF(A) (A&0x00ffff00)

typedef struct LxNzInf{

BYTE byColMask; // 1 << (j&7);

BYTE byBlkNz; // byBlkNz ++;

BYTE byBlkNzPos; // byBlkNzPos = j;

BYTE byRowMask; // 1 <<

(j>>3);

}LxNzInf;

typedef struct LxIdctInf{

DWORD dwNzInf; // += SET_NZ_INF

DWORD dwColMask; // |= SET_NZ_MASK

BYTE byRowPos[8]; // byRowPos[j>>3] |= msk[j]

; // 1 <<

(j&7);

}LxIdctInf;

__declspec (align(16)) const int

IDCT_TAB[8][8] =

{

{2048, 2048, 2048, 2048, 2048, 2048, 2048, 2048, },

{2841, 2408, 1609, 565,

-565, -1609, -2408, -2841, },

{2676, 1108, -1108, -2676, -2676, -1108, 1108, 2676, },

{2408, -565, -2841, -1609, 1609, 2841, 565,

-2408, },

{2048, -2048, -2048, 2048, 2048, -2048, -2048, 2048, },

{1609, -2841, 565,

2408, -2408, -565, 2841, -1609, },

{1108, -2676, 2676, -1108, -1108, 2676, -2676, 1108, },

{565, -1609, 2408, -2841, 2841, -2408, 1609, -565, },

};

extern DWORD g_dwBlockInf[64];

void Initialize_Fast_IDCT();

void idct_aan_bridge_ex(short*

block,short* dst, LxIdctInf* pos );

void idct_sparse_dc(short* block,short*

dst, LxIdctInf* pos);

void idct_sparse_ac(short* block,short*

dst, LxIdctInf* pos);

// end of LxIDCT.h .

// start of LxIDCT.cpp: implementation

of the LxIDCT .

static short iclip[1024];

short *iclp;

BYTE g_RowNz[256];

DWORD g_RowPos[256];

DWORD g_dwColMsk[8];

DWORD g_dwBlockInf[64];

#define

GET_COL_MASK(A) ( A & 0xff )

#define

GET_ROW_MASK(A) ( A

>> 24 )

#define

GET_NZ_NUM(A) ( (A >> 8) &

0xff )

#define

GET_NZ_POS(A) ( A >> 16 )

void Initialize_Fast_IDCT()

{

int i;

//填写裁减表

iclp = iclip+512;

for (i= -512; i<512;

i++) {

iclp[i] = (i

-256 : ((i>255) ? 255 : i);

}

#if 0

for( i = 0; i < 8; i ++

) {

double scale = (i == 0) ? sqrt(0.125) :

0.5;

for( int j = 0; j < 8; j

++ ) {

double s = scale * cos((PI/8.0)*i*(j +

0.5));

s = s * sqrt(2);

double delt = s >= 0 ?

0.5 : -0.5;

IDCT_TAB[i][j] = (short)(

(1<

}

#endif

for( i = 0; i < 256; i

++ ) {

int nz = 0;

g_RowPos[i] = 0;

for( int j = 0; j < 8; j

++ ) {

if( i &

(1<

g_RowPos[i] |=

nz ++;

}

g_RowNz[i] = nz;

}

for( i = 0; i < 8; i ++

) {

// g_RowPos[1<

g_dwColMsk[i] =

(1<

}

// DCT系数表。

for( i = 0; i < 8; i ++

) {

for( int j = 0; j < 8; j

++ ) {

LxNzInf inf;

inf.byColMask =

inf.byBlkNz = 1;

inf.byBlkNzPos = i*8+j;

inf.byRowMask =

g_dwBlockInf[i*8+j] =

*(DWORD*)&inf;

}

void idct_sparse_dc(short* block,short*

dst, LxIdctInf* pos)

{

int val32;

int val =

block[0];

val = ( val + 4)

>> 3;

val32 = (val

<< 16) | ( val &

0xffff );

int* b32 = (int*)dst;

for( int i = 0; i < 32;

i +=8 ){

b32[i+0] = val32;

b32[i+1] = val32;

b32[i+2] = val32;

b32[i+3] = val32;

b32[i+4] = val32;

b32[i+5] = val32;

b32[i+6] = val32;

b32[i+7] = val32;

}

void idct_sparse_ac( short* block,

short* dst, LxIdctInf* pos )

{

int tmp[8];

int i;

DWORD const dwRowMask =

GET_ROW_MASK(pos->dwColMask);

DWORD const dwColMask =

GET_COL_MASK(pos->dwColMask);

int hnc = g_RowPos[dwColMask];

int vnc = g_RowPos[dwRowMask];

int np = hnc*8 + vnc ;

if( 0 == np ) {

idct_sparse_dc(block,dst,pos);

return;

}

int ac = block[ np ];

int nc = hnc ? hnc : vnc;

tmp[0] = ( ( ac * IDCT_TAB[nc][0] +