dct8c.i

H.264完整的C语言代码和DCT的代码

// dct.i      (based on dct/dct.i)
// Ujval Kapasi
// 1/22/97
// 3/28/97
// 7/22/97
//
// 8x8 DCT (for JPEG and MPEG)
// 
// From Pennebaker/Mitchell, pg. 50-52.  See also Arai, Agui, Nakajima.
// This algorithm is based on the 16-pt DFT.  Basically, the 8-pt DCT can
//   be calculated by scaling the real parts of the output of the 16-pt DFT.
//
// This code performs two DCTs every iteration of the loop. Thus the input
// data must have one 8x8 block in the upper 16 bits of every word, and another
// 8x8 block in the lower 16 bits of every word.

kernel dct(istream<half2> consts,
           istream<half2> datain,
           ostream<int>   out)
{
  // DCT constants

  // Stored in 2.14 format
  // COS_2             = 0x2d412d41;    // cos(2*pi/8) || cos(2*pi/8);
  // COS_3             = 0x187e187e;    // cos(3*pi/8) || cos(3*pi/8);
  // COS_1_plus_COS_3  = 0x539f539f;    // cos(pi/8) + cos(3*pi/8) || same
  // COS_1_minus_COS_3 = 0x22a322a3;    // cos(pi/8) - cos(3*pi/8) || same

  half2 COS_2, COS_3, COS_1_plus_COS_3, COS_1_minus_COS_3;

  consts >> COS_2 >> COS_3 >> COS_1_plus_COS_3 >> COS_1_minus_COS_3;


  // Stored in 2.14 format
  // K0 = 0x16a116a1           // 0.25 * sqrt(2)       || 0.25 * sqrt(2);
  // K1 = 0x10501050           // 0.25 * sec(pi/16)    || 0.25 * sec(pi/16);
  // K2 = 0x11511151           // 0.25 * sec(2*pi/16)  || 0.25 * sec(2*pi/16);
  // K3 = 0x133e133e           // 0.25 * sec(3*pi/16)  || 0.25 * sec(3*pi/16);
  // K4 = 0x16a116a1           // 0.25 * sec(4*pi/16)  || 0.25 * sec(4*pi/16);
  // K5 = 0x1ccd1ccd           // 0.25 * sec(5*pi/16)  || 0.25 * sec(5*pi/16);
  // K6 = 0x29cf29cf           // 0.25 * sec(6*pi/16)  || 0.25 * sec(6*pi/16);
  // K7 = 0x52035203           // 0.25 * sec(7*pi/16)  || 0.25 * sec(7*pi/16);

  half2 K0, K1, K2, K3, K4, K5, K6, K7;

  consts >> K0 >> K1 >> K2 >> K3 >> K4 >> K5 >> K6 >> K7;


  // SP arrays  (not really persistent)

  array<half2> buf1(8); // intermediate dct output.  ie, do rows then
  array<half2> buf2(8); //   store here.  Then index into this
                        //   differently to get the columns

  // Comm permutations used to transpose the block

  uc<int> perm_a = 0x07654321;
  uc<int> perm_b = 0x10765432;
  uc<int> perm_c = 0x21076543;
  uc<int> perm_d = 0x32107654;
  uc<int> perm_e = 0x43210765;
  uc<int> perm_f = 0x54321076;
  uc<int> perm_g = 0x65432107;

  int src_idx = 0;
  int idx0 = cid();
  int idx1 = (idx0 - 1) & 7;
  int idx2 = (idx0 - 2) & 7;
  int idx3 = (idx0 - 3) & 7;
  int idx4 = (idx0 - 4) & 7;
  int idx5 = (idx0 - 5) & 7;
  int idx6 = (idx0 - 6) & 7;
  int idx7 = (idx0 - 7) & 7;

  loop_stream(datain) pipeline(1) {
    half2 a0, a1, a2, a3, a4, a5, a6, a7;

    datain >> a0 >> a1 >> a2 >> a3 >> a4 >> a5 >> a6 >> a7;

    // do the 1d dct
    half2 s16, s07, s25, s34, s1625, s0734;

    s07 = a0 + a7;
    s16 = a1 + a6;
    s25 = a2 + a5;
    s34 = a3 + a4;
    s1625 = s16 + s25;
    s0734 = s07 + s34;
    // 12 OPS (count double because we are using half2's)

    half2 d16, d07, d25, d34, d1625, d0734;

    d07 = a0 - a7;
    d16 = a1 - a6;
    d25 = a2 - a5;
    d34 = a3 - a4;
    d1625 = s16 - s25;
    d0734 = s07 - s34;
    // 12 OPS

    half2 sd16d07, sd25d34;

    sd16d07 = d07 + d16;
    sd25d34 = d25 + d34;
    // 4 OPS

    half2 m1_over_2, m2, m5, m6, m7, m8, m9;

    // All results in 16.0
    m1_over_2 = s0734 + s1625;
    m2 = s0734 - s1625;
    m5 = hi(COS_2 * shift(d1625 + d0734, 2));
    m6 = hi(COS_2 * shift(d25 + d16, 2));
    m7 = hi(COS_3 * shift(sd16d07 - sd25d34, 2));
    m8 = hi((COS_1_plus_COS_3) * shift(sd16d07, 2));
    m9 = hi((COS_1_minus_COS_3) * shift(sd25d34, 2));
    // 30 OPS

    half2 s5, s6, s7, s8;

    s5 = d07 + m6;
    s6 = d07 - m6;
    s7 = m8 - m7;
    s8 = m9 - m7;
    // 8 OPS

    // All results in 16.0
    buf1[0] = hi(K0 * shift(m1_over_2, 2));
    buf1[1] = hi(K1 * shift(s5 + s7, 2));
    buf1[2] = hi(K2 * shift(d0734 + m5, 2));
    buf1[3] = hi(K3 * shift(s6 - s8, 2));
    buf1[4] = hi(K4 * shift(m2, 2));
    buf1[5] = hi(K5 * shift(s6 + s8, 2));
    buf1[6] = hi(K6 * shift(d0734 - m5, 2));
    buf1[7] = hi(K7 * shift(s5 - s7, 2));
    // 44 OPS

    // Do comm stuff to transpose the matrix to do rows now
    
    buf2[idx0] = buf1[idx0];
    buf2[idx7] = commucperm(perm_a, buf1[idx1]);
    buf2[idx6] = commucperm(perm_b, buf1[idx2]);
    buf2[idx5] = commucperm(perm_c, buf1[idx3]);
    buf2[idx4] = commucperm(perm_d, buf1[idx4]);
    buf2[idx3] = commucperm(perm_e, buf1[idx5]);
    buf2[idx2] = commucperm(perm_f, buf1[idx6]);
    buf2[idx1] = commucperm(perm_g, buf1[idx7]);
    // 0 OPS
  
    // get a's from scratchpad -- In 16.0 format
    a0 = buf2[0];
    a1 = buf2[1];
    a2 = buf2[2];
    a3 = buf2[3];
    a4 = buf2[4];
    a5 = buf2[5];
    a6 = buf2[6];
    a7 = buf2[7];

    s07 = a0 + a7;
    s16 = a1 + a6;
    s25 = a2 + a5;
    s34 = a3 + a4;

    s1625 = s16 + s25;
    s0734 = s07 + s34;
    // 12 OPS

    d07 = a0 - a7;
    d16 = a1 - a6;
    d25 = a2 - a5;
    d34 = a3 - a4;
    d1625 = s16 - s25;
    d0734 = s07 - s34;
    // 12 OPS

    sd16d07 = d07 + d16;
    sd25d34 = d25 + d34;
    // 4 OPS

    // All results in 16.0
    m1_over_2 = s0734 + s1625;
    m2 = s0734 - s1625;
    m5 = hi(COS_2 * shift(d1625 + d0734, 2));
    m6 = hi(COS_2 * shift(d25 + d16, 2));
    m7 = hi(COS_3 * shift(sd16d07 - sd25d34, 2));
    m8 = hi((COS_1_plus_COS_3) * shift(sd16d07, 2));
    m9 = hi((COS_1_minus_COS_3) * shift(sd25d34, 2));
    // 30 OPS

    s5 = d07 + m6;
    s6 = d07 - m6;
    s7 = m8 - m7;
    s8 = m9 - m7;
    // 8 OPS

    half2 d0, d1, d2, d3, d4, d5, d6, d7;

    // All results in 16.0
    d0 = hi(K0 * shift(m1_over_2, 2));
    d1 = hi(K1 * shift(s5 + s7, 2));
    d2 = hi(K2 * shift(d0734 + m5, 2));
    d3 = hi(K3 * shift(s6 - s8, 2));
    d4 = hi(K4 * shift(m2, 2));
    d5 = hi(K5 * shift(s6 + s8, 2));
    d6 = hi(K6 * shift(d0734 - m5, 2));
    d7 = hi(K7 * shift(s5 - s7, 2));
    // 44 OPS

    // TOTAL :  220 OPS per lop iter per cluster
    //          1760 OPS TOTAL for 8 clusters

    ///////////////////////////////////
    //      -->  110 OPS per BLOCK  //
    //           880   OPS total     //
    ///////////////////////////////////

    out << d0 << d1 << d2 << d3 << d4 << d5 << d6 << d7;
  }
}