fdct_sse2_skal.asm

最新的视频开发包哦

;/****************************************************************************
; *
; *  XVID MPEG-4 VIDEO CODEC
; *  - SSE2 forward discrete cosine transform -
; *
; *  Copyright(C) 2003 Pascal Massimino <[e-mail]>
; *
; *  This program is free software; you can redistribute it and/or modify it
; *  under the terms of the GNU General Public License as published by
; *  the Free Software Foundation; either version 2 of the License, or
; *  (at your option) any later version.
; *
; *  This program is distributed in the hope that it will be useful,
; *  but WITHOUT ANY WARRANTY; without even the implied warranty of
; *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
; *  GNU General Public License for more details.
; *
; *  You should have received a copy of the GNU General Public License
; *  along with this program; if not, write to the Free Software
; *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
; *
; * $Id$
; *
; ***************************************************************************/

BITS 32

%macro cglobal 1
  %ifdef PREFIX
    global _%1
    %define %1 _%1
  %else
    global %1
  %endif
%endmacro

;-----------------------------------------------------------------------------
;
;                          -=FDCT=-
;
; Vertical pass is an implementation of the scheme:
;  Loeffler C., Ligtenberg A., and Moschytz C.S.:
;  Practical Fast 1D DCT Algorithm with Eleven Multiplications,
;  Proc. ICASSP 1989, 988-991.
;
; Horizontal pass is a double 4x4 vector/matrix multiplication,
; (see also Intel's Application Note 922:
;  http://developer.intel.com/vtune/cbts/strmsimd/922down.htm
;  Copyright (C) 1999 Intel Corporation)
;  
; Notes:
;  * tan(3pi/16) is greater than 0.5, and would use the
;    sign bit when turned into 16b fixed-point precision. So,
;    we use the trick: x*tan3 = x*(tan3-1)+x
; 
;  * There's only one SSE-specific instruction (pshufw).
;
;  * There's still 1 or 2 ticks to save in fLLM_PASS, but
;    I prefer having a readable code, instead of a tightly
;    scheduled one...
;
;  * Quantization stage (as well as pre-transposition for the
;    idct way back) can be included in the fTab* constants
;    (with induced loss of precision, somehow)
;
;  * Some more details at: http://skal.planet-d.net/coding/dct.html
;
;
;//////////////////////////////////////////////////////////////////////
;
;  == Mean square errors ==
;   0.000 0.001 0.001 0.002 0.000 0.002 0.001 0.000    [0.001]
;   0.035 0.029 0.032 0.032 0.031 0.032 0.034 0.035    [0.032]
;   0.026 0.028 0.027 0.027 0.025 0.028 0.028 0.025    [0.027]
;   0.037 0.032 0.031 0.030 0.028 0.029 0.026 0.031    [0.030]
;   0.000 0.001 0.001 0.002 0.000 0.002 0.001 0.001    [0.001]
;   0.025 0.024 0.022 0.022 0.022 0.022 0.023 0.023    [0.023]
;   0.026 0.028 0.025 0.028 0.030 0.025 0.026 0.027    [0.027]
;   0.021 0.020 0.020 0.022 0.020 0.022 0.017 0.019    [0.020]
;  
;  == Abs Mean errors ==
;   0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000    [0.000]
;   0.020 0.001 0.003 0.003 0.000 0.004 0.002 0.003    [0.002]
;   0.000 0.001 0.001 0.001 0.001 0.004 0.000 0.000    [0.000]
;   0.027 0.001 0.000 0.002 0.002 0.002 0.001 0.000    [0.003]
;   0.000 0.000 0.000 0.000 0.000 0.001 0.000 0.001    [-0.000]
;   0.001 0.003 0.001 0.001 0.002 0.001 0.000 0.000    [-0.000]
;   0.000 0.002 0.002 0.001 0.001 0.002 0.001 0.000    [-0.000]
;   0.000 0.002 0.001 0.002 0.001 0.002 0.001 0.001    [-0.000]
;
;  =========================
;  Peak error:   1.0000
;  Peak MSE:     0.0365
;  Overall MSE:  0.0201
;  Peak ME:      0.0265
;  Overall ME:   0.0006
;
;-----------------------------------------------------------------------------
;
;                          -=IDCT=-
;
; A little slower than fdct, because the final stages (butterflies and
; descaling) require some unpairable shifting and packing, all on
; the same CPU unit.
;
;   THIS IDCT IS NOT IEEE-COMPLIANT: IT WILL FAIL THE [-300,300]
;   INPUT RANGE TEST (because of overflow). But the [-256,255] one
;   is OK, and I'm fine with it (for now;)
;
;  == Mean square errors ==
;   0.007 0.006 0.005 0.007 0.006 0.007 0.005 0.007    [0.006]
;   0.006 0.008 0.007 0.007 0.007 0.008 0.008 0.008    [0.007]
;   0.008 0.008 0.008 0.008 0.007 0.009 0.010 0.007    [0.008]
;   0.007 0.007 0.006 0.007 0.008 0.007 0.006 0.008    [0.007]
;   0.007 0.006 0.006 0.006 0.006 0.005 0.006 0.006    [0.006]
;   0.008 0.007 0.006 0.008 0.007 0.008 0.009 0.009    [0.008]
;   0.008 0.006 0.010 0.008 0.008 0.008 0.007 0.007    [0.008]
;   0.007 0.006 0.006 0.007 0.007 0.006 0.006 0.007    [0.006]
;  
;  == Abs Mean errors ==
;   0.001 0.000 0.000 0.001 0.001 0.000 0.000 0.000    [0.000]
;   0.000 0.002 0.002 0.000 0.001 0.001 0.000 0.002    [0.000]
;   0.001 0.002 0.001 0.001 0.001 0.001 0.000 0.001    [-0.001]
;   0.000 0.002 0.000 0.000 0.001 0.000 0.000 0.001    [-0.000]
;   0.000 0.001 0.001 0.001 0.000 0.001 0.000 0.001    [0.000]
;   0.000 0.001 0.001 0.001 0.001 0.000 0.001 0.000    [0.000]
;   0.001 0.001 0.002 0.001 0.001 0.002 0.001 0.001    [0.001]
;   0.000 0.000 0.001 0.000 0.000 0.000 0.000 0.000    [0.000]
;
;  =========================
;
;  Peak error:   1.0000
;  Peak MSE:     0.0096
;  Overall MSE:  0.0070
;  Peak ME:      0.0024
;  Overall ME:   0.0001
;
;-----------------------------------------------------------------------------

;=============================================================================
; Read only data
;=============================================================================

%ifdef FORMAT_COFF
SECTION .rodata
%else
SECTION .rodata align=16
%endif

ALIGN 16
tan1:    times 8 dw 0x32ec    ; tan( pi/16)
tan2:    times 8 dw 0x6a0a    ; tan(2pi/16)  (=sqrt(2)-1)
tan3:    times 8 dw 0xab0e    ; tan(3pi/16)-1
sqrt2:   times 8 dw 0x5a82    ; 0.5/sqrt(2)

;-----------------------------------------------------------------------------
; Inverse DCT tables
;-----------------------------------------------------------------------------

ALIGN 16
iTab1:
  dw 0x4000, 0x539f, 0x4000, 0x22a3
  dw 0x4000, 0xdd5d, 0x4000, 0xac61
  dw 0x4000, 0x22a3, 0xc000, 0xac61
  dw 0xc000, 0x539f, 0x4000, 0xdd5d
  dw 0x58c5, 0x4b42, 0x4b42, 0xee58
  dw 0x3249, 0xa73b, 0x11a8, 0xcdb7
  dw 0x3249, 0x11a8, 0xa73b, 0xcdb7
  dw 0x11a8, 0x4b42, 0x4b42, 0xa73b

iTab2:
  dw 0x58c5, 0x73fc, 0x58c5, 0x300b
  dw 0x58c5, 0xcff5, 0x58c5, 0x8c04
  dw 0x58c5, 0x300b, 0xa73b, 0x8c04
  dw 0xa73b, 0x73fc, 0x58c5, 0xcff5
  dw 0x7b21, 0x6862, 0x6862, 0xe782
  dw 0x45bf, 0x84df, 0x187e, 0xba41
  dw 0x45bf, 0x187e, 0x84df, 0xba41
  dw 0x187e, 0x6862, 0x6862, 0x84df

iTab3:
  dw 0x539f, 0x6d41, 0x539f, 0x2d41
  dw 0x539f, 0xd2bf, 0x539f, 0x92bf
  dw 0x539f, 0x2d41, 0xac61, 0x92bf
  dw 0xac61, 0x6d41, 0x539f, 0xd2bf
  dw 0x73fc, 0x6254, 0x6254, 0xe8ee
  dw 0x41b3, 0x8c04, 0x1712, 0xbe4d
  dw 0x41b3, 0x1712, 0x8c04, 0xbe4d
  dw 0x1712, 0x6254, 0x6254, 0x8c04

iTab4:
  dw 0x4b42, 0x6254, 0x4b42, 0x28ba
  dw 0x4b42, 0xd746, 0x4b42, 0x9dac
  dw 0x4b42, 0x28ba, 0xb4be, 0x9dac
  dw 0xb4be, 0x6254, 0x4b42, 0xd746
  dw 0x6862, 0x587e, 0x587e, 0xeb3d
  dw 0x3b21, 0x979e, 0x14c3, 0xc4df
  dw 0x3b21, 0x14c3, 0x979e, 0xc4df
  dw 0x14c3, 0x587e, 0x587e, 0x979e

  ; the original rounding trick is by
  ; Michel Lespinasse (hi Walken!) <[e-mail]>

ALIGN 16
Idct_Rnd0: dd  65535, 65535, 65535, 65535
Idct_Rnd1: dd   3612,  3612,  3612,  3612
Idct_Rnd2: dd   2271,  2271,  2271,  2271
Idct_Rnd3: dd   1203,  1203,  1203,  1203
Idct_Rnd4: dd   1023,  1023,  1023,  1023
Idct_Rnd5: dd    102,   102,   102,   102
Idct_Rnd6: dd    398,   398,   398,   398
Idct_Rnd7: dd    469,   469,   469,   469

Idct_Sparse_Rnd0: times 4 dw  (65535>>11)
Idct_Sparse_Rnd1: times 4 dw  ( 3612>>11)
Idct_Sparse_Rnd2: times 4 dw  ( 2271>>11)
  ; other rounders are zero...

;-----------------------------------------------------------------------------
; Forward DCT tables
;-----------------------------------------------------------------------------

ALIGN 16
fTab1:
  dw 0x4000, 0x4000, 0x58c5, 0x4b42,
  dw 0xdd5d, 0xac61, 0xa73b, 0xcdb7,
  dw 0x4000, 0x4000, 0x3249, 0x11a8,
  dw 0x539f, 0x22a3, 0x4b42, 0xee58,
  dw 0x4000, 0xc000, 0x3249, 0xa73b,
  dw 0x539f, 0xdd5d, 0x4b42, 0xa73b,
  dw 0xc000, 0x4000, 0x11a8, 0x4b42,
  dw 0x22a3, 0xac61, 0x11a8, 0xcdb7

fTab2:
  dw 0x58c5, 0x58c5, 0x7b21, 0x6862,
  dw 0xcff5, 0x8c04, 0x84df, 0xba41,
  dw 0x58c5, 0x58c5, 0x45bf, 0x187e,
  dw 0x73fc, 0x300b, 0x6862, 0xe782,
  dw 0x58c5, 0xa73b, 0x45bf, 0x84df,
  dw 0x73fc, 0xcff5, 0x6862, 0x84df,
  dw 0xa73b, 0x58c5, 0x187e, 0x6862,
  dw 0x300b, 0x8c04, 0x187e, 0xba41

fTab3:
  dw 0x539f, 0x539f, 0x73fc, 0x6254,
  dw 0xd2bf, 0x92bf, 0x8c04, 0xbe4d,
  dw 0x539f, 0x539f, 0x41b3, 0x1712,
  dw 0x6d41, 0x2d41, 0x6254, 0xe8ee,
  dw 0x539f, 0xac61, 0x41b3, 0x8c04,
  dw 0x6d41, 0xd2bf, 0x6254, 0x8c04,
  dw 0xac61, 0x539f, 0x1712, 0x6254,
  dw 0x2d41, 0x92bf, 0x1712, 0xbe4d

fTab4:
  dw 0x4b42, 0x4b42, 0x6862, 0x587e,
  dw 0xd746, 0x9dac, 0x979e, 0xc4df,
  dw 0x4b42, 0x4b42, 0x3b21, 0x14c3,
  dw 0x6254, 0x28ba, 0x587e, 0xeb3d,
  dw 0x4b42, 0xb4be, 0x3b21, 0x979e,
  dw 0x6254, 0xd746, 0x587e, 0x979e,
  dw 0xb4be, 0x4b42, 0x14c3, 0x587e,
  dw 0x28ba, 0x9dac, 0x14c3, 0xc4df


ALIGN 16
Fdct_Rnd0: dw  6,8,8,8, 6,8,8,8
Fdct_Rnd1: dw  8,8,8,8, 8,8,8,8
Fdct_Rnd2: dw 10,8,8,8, 8,8,8,8
Rounder1:  dw  1,1,1,1, 1,1,1,1

;=============================================================================
; Code
;=============================================================================

SECTION .text

cglobal idct_sse2_skal
cglobal idct_sse2_sparse_skal
cglobal fdct_sse2_skal

;-----------------------------------------------------------------------------
; Helper macro iMTX_MULT
;-----------------------------------------------------------------------------

%macro iMTX_MULT 4   ; %1=src, %2 = Table to use, %3=rounder, %4=Shift

  movdqa  xmm0, [ecx+%1*16]     ; xmm0 = [01234567]

  pshuflw xmm0, xmm0, 11011000b ; [0213]
  pshufhw xmm0, xmm0, 11011000b ; [02134657]
  pshufd  xmm4, xmm0, 00000000b ; [02020202]
  pshufd  xmm5, xmm0, 10101010b ; [46464646]
  pshufd  xmm6, xmm0, 01010101b ; [13131313]
  pshufd  xmm7, xmm0, 11111111b ; [57575757]

  pmaddwd xmm4, [%2+ 0]   ; dot [M00,M01][M04,M05][M08,M09][M12,M13]
  pmaddwd xmm5, [%2+16]   ; dot [M02,M03][M06,M07][M10,M11][M14,M15]
  pmaddwd xmm6, [%2+32]   ; dot [M16,M17][M20,M21][M24,M25][M28,M29]
  pmaddwd xmm7, [%2+48]   ; dot [M18,M19][M22,M23][M26,M27][M30,M31]
  paddd   xmm4, [%3]      ; Round

  paddd   xmm6, xmm7      ; [b0|b1|b2|b3]
  paddd   xmm4, xmm5      ; [a0|a1|a2|a3]

  movdqa  xmm7, xmm6
  paddd   xmm6, xmm4      ; mm6=a+b
  psubd   xmm4, xmm7      ; mm4=a-b
  psrad   xmm6, %4        ; => out [0123]
  psrad   xmm4, %4        ; => out [7654]

  packssdw xmm6, xmm4     ; [01237654]
  pshufhw xmm6, xmm6, 00011011b ; [01234567]

  movdqa  [ecx+%1*16], xmm6

%endmacro

;-----------------------------------------------------------------------------
; Helper macro iLLM_PASS
;-----------------------------------------------------------------------------

%macro iLLM_PASS 1  ; %1: src/dst

  movdqa xmm0, [tan3]    ; t3-1
  movdqa xmm3, [%1+16*3] ; x3
  movdqa xmm1, xmm0       ; t3-1
  movdqa xmm5, [%1+16*5] ; x5

  movdqa xmm4, [tan1]    ; t1
  movdqa xmm6, [%1+16*1] ; x1
  movdqa xmm7, [%1+16*7] ; x7
  movdqa xmm2, xmm4       ; t1

  pmulhw xmm0, xmm3       ; x3*(t3-1)
  pmulhw xmm1, xmm5       ; x5*(t3-1)
  paddsw xmm0, xmm3       ; x3*t3
  paddsw xmm1, xmm5       ; x5*t3
  psubsw xmm0, xmm5       ; x3*t3-x5 = tm35
  paddsw xmm1, xmm3       ; x3+x5*t3 = tp35

  pmulhw xmm4, xmm7       ; x7*t1
  pmulhw xmm2, xmm6       ; x1*t1
  paddsw xmm4, xmm6       ; x1+t1*x7 = tp17
  psubsw xmm2, xmm7       ; x1*t1-x7 = tm17