|
文 摘： 修正离散余弦变换（MDCT）在音视频信号编码中得到广泛地应用，其快速算法在实时编解码系统中尤为重要。论文给出了一种适用于数字信号处理器（DSP）实现的修正离散余弦反变换（IMDCT）快速算法—用M/2点时间抽取（decimation in time, DIT）分裂基FFT实现2M点的IMDCT。算法是基于蝶形运算组成，在DSP中可以获得很高的运算效率。该算法的蝶形运算结构同样适用于正向MDCT。在由定点DSP实现的活动图像专家组(MPEG)音频层III解码器中，与MPEG音频压缩标准ISO/IEC 11172-3中给出的IMDCT运算量相比较，该文提出的IMDCT快速算法节省了2/3的运算时间和1/2的存储空间。
关键词： 数字信号处理器（DSP）； 修正离散余弦变换（MDCT）； MPEG； 音频压缩
中图分类号： TN 911； TN 912 文献标识码： A
文章编号： 1000-0054(2000)03-0099-05
One of the Fast DSP-based
inverse MDCT algorithm
DOU Weibei, LIU Ruoheng, WANG Jianxin, DONG Zaiwang
(Department of Electronic Engineering， Tsinghua University, Beijing 100084, China)
Abstract： Modified discrete cosine transform (MDCT) is widely used in audio compression, such as the MPEG audio layer III codec. For the MDCT algorithm on a fixed-point digital signal processor (DSP), a butterfly operation of inverse MDCT was proposed to reduce the computational time and memory requirements. The 2M-length inverse MDCT algorithm was derived in the form of a butterfly computation with M/2 points DFT and M points windowing. This fast algorithm operating in ADSP-2181 has at least 66% less instructions and 50% less memory space then the inverse MDCT described in ISO/IEC 11172-3.
Key words： digital signal processor (DSP); inverse modified discrete cosine transform; MPEG; audio compression
离散余弦变换（DCT）常被认为是对语音与图像信号进行变换的好方法，其变换特性接近KLT(karhunen-loève transform)。1986年由Prencen和Bradly提出的一种改进的离散余弦变换（MDCT）［1］，利用50%的样点重叠和时域混叠消除（TDAC）滤波器组，在不降低变换编码性能的情况下，有效地克服了DCT块处理运算中的边缘效应（块边缘噪声）［2］。在相同编码效率的情况下， MDCT的性能优于DCT，目前在语音、宽带音频及图像信号的变换编码中，都普遍采用MDCT。
随着大规模集成电路技术的发展，高速度、低功耗、大内存、低成本的数字信号处理器DSP大量面世，用DSP实现的音频、视频信号的实时编解码系统也应运而生。用DSP实现实时系统时，由于DSP的性能约束，需要尽可能地减小运算时间开销和存储空间的占用量，并且要很好地解决二者之间的矛盾。在MDCT的计算中，由于每个样点需要变换两次，因此，直接计算MDCT所需要的计算时间较长，占用的存储空间也很大。在用DSP实现的实时编解码系统中， MDCT的快速算法是非常重要的。文［3］提出将MDCT转化为标准DCT，再利用快速余弦变换FCT进行计算。但该算法主要适用于变换长度为2m的MDCT。文［2］提出利用余弦函数的对称性，将2M点的MDCT计算分成数据移位、预处理、 M点快速Fourier变换和后处理来实现快速MDCT计算。文［4］提出了适用于并行VLSI实现的递归MDCT算法。本文从定点DSP芯片的特点出发，设计了蝶形反向MDCT（Batterfly-Inverse MDCT, B-IMDCT）算法，减小了MDCT运算中的存储空间和运算时间开销。算法采用了复序列的FFT技术，使块处理长度减小了一半。针对MDCT具有不同窗长的特点，算法中采用时间抽取（DIT）分裂基FFT(Radix-Splitted FFT)，以适应任意长度的IMDCT，进一步提高运算效率。与文［5］的算法相比，将B-IMDCT算法应用于DSP实现的MPEG音频层III解码器中，使IMDCT的运算量减小了2/3，存储空间节省了一半。本文设计的蝶形运算具有完全的对称性，同样适用于正向MDCT。
1 IMDCT的简化
MDCT的表达式为：
正变换
 (1.a)
反变换
 (n)=z(n)+z(n+M), n=0,1,…，M-1， (1.b)
 (1.c)
其中： 变换核
 (1.d)
2M称为MDCT的变换窗长，在不同的应用中变换窗长度不一样；
h(n)为窗函数，满足如下完全重构（PR）条件［6］:
 (2)
通常，定义
 (3)
由式（1.d）时移3M/2，可得
定义：
则：
bnk=-a(n-3M/2)k， (5.a)
 (5.b)
式（1.b）是IMDCT的最后一步重叠相加计算，由式（1.c）和（2），（1.b）可以写成：
 (6)
将式（5.b）代入（6），得到：］
 (7)
并根据式(4.a)的特性得到：
则式(7)变成：
 (8)
式(8)是蝶形计算式，该算式的蝶形结构如图1所示。上述推导说明： 2M点的IMDCT可利用M点的蝶形算式计算得到。
|
图1 IMDCT的加窗蝶形结构
2 IMDCT的FFT实现
IMDCT的化简过程将式(1.c)计算2M点的y(n)转化成了计算M点的g(n)，这样，式（4.b）变成
 (9)
式(9)是第4类DCT(DCT-IV)表达式［7］，可利用式(9)的变换核式(4.a)与DFT核
 (10)
的关系，将式（9）转换成Cooley-Tukey FFT实现。
对于一个给定的输入实序列［X0,X1,…，XM-1］，可定义两个新的向量：
 (11)
则式（9）可以写成
 (12)
其中：
 (13)
那么利用式（10），可得到关系式：
 (14)
将式(14)代入(12)，得到：
 (15)
式(15)可以被看成 (Ek+jOk)W4k+18M的M/2点的DFT。至此，由式(8)和(15)，将一个长度为M的输入实序列组合成一个长度为M/2的复序列，从而用M/2点的复数DFT实现了2M点的IMDCT。
3 IMDCT蝶形快速算法
通过上面的推导，可以得出结论，IMDCT可由DFT和蝶形加窗运算完成。而DFT很容易被设计成同位蝶形运算FFT。在 DSP中，DFT的蝶形计算过程可分为级、组和基本蝶形三重循环。由三点构成的基本蝶形DFT运算被称之为基3 FFT（Radix-3 FFT），图2是由DIT基3 FFT的基本蝶形运算和叠形加窗运算级联组成的36点IMDCT计算过程。它表明，不论输入是实序列还是复序列， FFT的运算总是复数运算过程。因此，将一个长度为M的输入实序列组合成一个长度为M/2的复序列，并没有增加运算复杂度，而是在充分利用FFT的运算模块的同时，又将块处理长度减小了一半，从而大大提高了IMDCT的运算速度。
如图2所示，一个输入长度为M的实序列， {X(k), k=0,1,…,M-1} ，其窗长为2M的IMDCT可由适于同位运算的IMDCT蝶形快速算法计算：
第一步： 输入序列重新排序：
X(k)←X(2k)， |
|
图2 36点IMDCT蝶形运算图
X(k+M/2)←X(M-1-2k)，
k＝0，1，…，M/2-1.
第二眇:预处理,乘旋转因子
第三步： 对复序列作M/2点DFT, n=0,1,…,M/2-1.
g(n)←Re{DFT{［X(0),X(1)，…，X(M-1)］,
M/2}},
g(n+M/2)←Im{DFT{［X(0),X(1)，…，
X(M-1)］, M/2}}.
第四步： 后处理，对复序列乘旋转因子，重排序为
g(n)←Re{g(n)exp(-jnπ/M)},
g(M-2n-1)←-Im{g(n)exp(-jnπ/M)}.
第五步： 进行蝶形加窗运算：
(n)←g(n+M/2)h(n)-g(M/2-1-n).
h(M-1-n),
(M/2-1-n)←-g(M/2-1-n)h(n)-
g(n+M/2)h(M-1-n).
这就是由蝶形FFT和蝶形加窗运算两部分级联构成的蝶形IMDCT快速算法。由于窗函数h(n)的特性（见式(2)），该算法完全适用于正向MDCT。
4 应用举例
将上述算法应用于用单片定点DSP实现MPEG音频层III解码。文［5］中IMDCT窗长为36和12，去混叠运算模块输出每子带18个样值，对应32个子带中各子带的窗长进行相应的IMDCT。利用上述蝶形IMDCT算法，对于长窗2M=36，需要作9点的FFT和18点的加窗蝶形运算； 而对于短窗只需作3点的FFT和6点的加窗运算。如果选用基3 FFT，无论长窗还是短窗都可以共用3点的基本蝶形运算模块。短窗时只需一级基本蝶形运算，而长窗时需要两级循环，第一级有三组各一个基本蝶形，第二级只有一组共3个基本蝶形运算。采用Analog Device Inc.(ADI)定点DSP—ADSP2181实现这一运算过程，应用块浮点技术以提高定点运算的精度，其计算所需的指令周期如表1所示。其中用作对比的计算量，是在同样的DSP应用技术下，由文［5］给出的标准计算过程所需的计算量。
表1 蝶形IMDCT运算量比较 |
窗长
(样点数) |
蝶形IMDCT
运算量
B（指令周期数） |
标准IMDCT［5］
运算量
S（指令周期数） |
运算量比较
B.S－1／% |
| 12 |
152 |
173 |
87.5 |
| 36 |
525 |
1552 |
33.8 |
由于ADSP-21xx芯片所特有的并行处理能力与乘加运算指令已经简化了运算量，所以在短窗时，蝶形IMDCT算法的优势不太明显，只比对比算法快12.5% 。随着窗长度的增加，蝶形IMDCT算法的优势就愈加明显。窗长为36时，比对比算法快了66.2%.
蝶形IMDCT不仅在运算量上明显优于对比算法，而且，由于计算长度减少和同位运算的应用，IMDCT运算中的存储空间也相应的减少了1/2，在空间资源的节省方面也有明显的优势。
5 结束语
近来，MDCT被广泛应用于变换编码中，由于其运算量大，对于一个实用的编码系统而言，MDCT的快速算法显得尤为重要。本文根据定点DSP的特点，利用时间抽取(DIT)分裂基快速Fourier变换，设计了适用于定点DSP的蝶形IMDCT（Inverse MDCT）快速算法，该算法完全适用于正向MDCT。在由单片定点ADSP2181实现的MPEG音频层III解码器中，蝶形IMDCT快速算法的应用节省了2/3的运算量。该算法适用于任意长度的MDCT，随着长度的增加，运算量的节省优势越明显。
基金项目： 国家科委高技术发展研究委托项目
作者简介： 窦维蓓(1963-), 女(汉), 重庆， 副教授
窦维蓓(清华大学 电子工程系， 北京 100084)
刘若珩(清华大学 电子工程系， 北京 100084)
王建昕(清华大学 电子工程系， 北京 100084)
董在望(清华大学 电子工程系， 北京 100084)
［参考文献］
［1］Princen J, Bradley A. Analysis/synthesis filter bank design based on time domain aliasing cancellation［J］. IEEE Transaction Acoustic, Speech, Signal Processing, 1986, 34(5): 11531161.
［2］Masahiro Iwadare, Akihiko Sugiyama, Fumie Hazu, et al. A 128kb/s Hi-Fi audio CODEC based on adaptive transform coding with adaptive block size MDCT［J］. IEEE Journal on Selected Areas in Communications, 1992, 10(1): 138144.
［3］Wang Jianxin, Dong Zaiwang. A fast algorithm for modified discrete cosine transform［A］. 1996 International Conference on Communication Technology Proceedings(ICCT'96)［C］, Beijing:PHEI & IEEE PRESS(Publishing House of Electronics Industy & Institute of Electrical and Electronics Engineers), 1996. 445448.
［4］Chiang Hwang-cheng, Liu Jie-Cherng. Regressive Implementations for the Forward and Inverse MDCT in MPEG Audio Coding［J］. IEEE Signal Processing Letters, 1996, 3(4): 116118.
［5］ISO/IEC 11172-3. Information technology—coding of moving pictures and associated audio for digital storage media at up to about 1.5Mbits/s—part3 Audio［S］.
［6］Malvar Henrique S. Signal Processing with Lapped Transforms［M］, Norwood, MA: Artech House, 1992. ［7］Wang Z. Fast algorithms for the discrete W transform and for the discrete Fourier transform［J］. IEEE Transaction Acoustic, Speech, Signal Processing, 1984, 32(4): 803816. |
|
