mpeg 变换DCT

背景

uncorrelated image has more sharp intensity vairations-->

has more high frequency content.
has its energy spread out
low spatial content

energy compactiion<-->correlation

定义

1维N点FDCT和IDCT:
 C(u) =\sum_{x=0}^{N-1}  a(u) f(x) cos(P(2x+1)u/(2N))
 f(x) =\sum_{u=0}^{N-1}  a(u) C(u) cos(P(2x+1)u/(2N))

2维N点FDCT和IDCT:
 C(u,v) =\sum_{x=0}^{N-1}\sum_{x=0}^{N-1}  
           a(u)a(v) f(x,y) cos(P(2x+1)u/(2N)) cos(P(2y+1)v/(2N))
 f(x,y) =\sum_{u=0}^{N-1}\sum_{v=0}^{N-1}  
           a(u)a(v) C(u,v) cos(P(2x+1)u/(2N)) cos(P(2y+1)v/(2N))
x,y,u,v=0..N-1

可以看出，只要记住1维N点FDCT，其他3个类似。
另外，2维N点FDCT可以理解为
1. 先行变换，对y求和:f(x,y) --> C(x,v)，
然后列变换，对x求和:C(x,v) --> C(u,v)

2. 矩阵形式 C = AfA'
矩阵A的i,j项为 a(i)cos(P(2j+1)i/(2N))
a(i)=\sqrt(1/N) 若i=0，否则\sqrt(2/N)

注意到这个矩阵是正交的: A'A = E，因此有
 f = A'CA

dctref.c 就是按两个矩阵形式计算的，需要注意的是
ff_ref_fdct()结果大了8倍，目的是与其他优化版本直接比较。
其中有个细节
block[i + j] = floor(tmp + 0.499999999999);
为什么不是0.5呢？

3.图形化显示
我见到的文档图形都是灰色的，如下面的[1]，打算写个ppm, 暖色调表示正值，冷色调表示负值。不过转念一想，Y分量自然是黑白的，没必要。

实现代码
-dct -idct的参数见 options_tables.h。默认是mmx和sse2,
纯C整数版本的有int, fastint, int, simple, 分别对应这些文件:
jfdctint_template.c -- 使用LLM算法[2]，12个乘法，慢但是准确
jfdctfst.c -- 使用AAN算法[4]，5个乘法，post-scale
jrevdct.c -- 使用LLM算法，优化稀疏结构(5个if嵌套)。
simple_idct_template.c -- mmx的C版本。

测试
(环境:core i3 ubuntu 14.04 64bit 2GB memeory)
make libavcodec/dct-test && cd libavcodec
$ ./dct-test  -t | grep kdct | sort -nuk 3
DCT REF-DBL: 447.2 kdct/s
DCT FAAN: 2781.0 kdct/s
DCT IJG-LLM-INT: 5146.8 kdct/s
DCT IJG-AAN-INT: 5902.2 kdct/s
DCT MMX: 9386.2 kdct/s
DCT MMXEXT: 11020.4 kdct/s
DCT SSE2: 17417.6 kdct/s
$ ./dct-test -i -t | grep kdct | sort -nuk 3
IDCT REF-DBL: 1235.3 kdct/s
IDCT FAANI: 1863.7 kdct/s
IDCT PR-C: 2714.8 kdct/s
IDCT PR-SSE2: 4737.0 kdct/s
IDCT SIMPLE-C: 6712.0 kdct/s
IDCT INT: 7281.5 kdct/s
IDCT XVID-MMX: 9729.4 kdct/s
IDCT SIMPLE-MMX: 10373.5 kdct/s
IDCT XVID-MMXEXT: 10913.4 kdct/s
IDCT XVID-SSE2: 25000.6 kdct/s

IMDCT变换
基于人耳的带通特性,人们开始使用FFT和PQMF进行带通分析,后来逐步通过MDCT进行取代.主要问
题是带通分析不可避免的带入频带交叠误差.这也是频带处理技术引入的主要误差.技术发展是从早期的大运算量
的FFT和PQMF到小运算量且容易进行频带交叠误差处理的MDCT技术.更多见[3]。
变换技术主要过程：加窗,变换,去交叠处理
音频解码场景: r -i in.mp3 out.mp2
#0  ff_imdct36_blocks_fixed at libavcodec/mpegaudiodsp_template.c:385
#1  compute_imdct at libavcodec/mpegaudiodec_template.c:1243
#2  mp_decode_layer3  at libavcodec/mpegaudiodec_template.c:1540
#3  mp_decode_frame  at libavcodec/mpegaudiodec_template.c:1572
#4  decode_frame at libavcodec/mpegaudiodec_template.c:1684
#5  avcodec_decode_audio4 at libavcodec/utils.c:2432
#6  try_decode_frame   at libavformat/utils.c:2759

#7  avformat_find_stream_info at libavformat/utils.c:3376

熵解码
huffman_decode
反量化
exponents_from_scale_factors
反变换
compute_imdct

杂项
用高精度计算器来验证常数FIX_xxx是否正确:
echo "0.298631336*2^13 + 0.5" | bc #FIX_0_298631336
pi=$(echo "scale=10; 4*a(1)" | bc -l)
echo "sqrt(2)*c($pi*6/16)" | bc -l #FIX_0_541196100

Ref

[1] https://en.wikipedia.org/wiki/Discrete_cosine_transform

[2] C.Loeffler, A.Ligtenberg, G.Moschytz Practical fast 1-D DCT Alogirthms with 11 Multiplications. 1989 IEEE

[3] http://people.xiph.org/~xiphmont/demo/daala/demo1.shtml.

[4] http://vsr.informatik.tu-chemnitz.de/~jan/MPEG/HTML/IDCT.html