快速小波变换在DSP中的实现方法.docx

快速小波变换在DSP中的实现方法2008-03-25 09:21姜新华范征宇（上海交通大学信息与控制工程系上海，200240）摘要小波分析是分析非稳定信号的一种非常有效的方法。Mallat快速算法使得小波分析的广泛应用成为现实。在实时信号分析中，离散小波变换在DSP上的有效应用受到了特别的关注。本文简单介绍了FWT，详细阐述了在DSP上FWT的周期性扩展的实现。其中，特别介绍了DSP的循环寻址。最后，给出了针对TI公司的TMS320C3X系列DSP的相应的汇编代码。关键词：小波变换；快速算法；数字信号处理器引言小波变换是从信号理论领域发展形成的。小波变换涵盖许多领域，它在这些领域得到了（或将得到）广泛的应用。作为信号处理的经典工具，傅里叶变换的缺点是它没有时域局部化性质。它仅仅给出信号的频率成分表达，但一般不能从变换后的信号中得出频率与时域局部有关的信息，其原因就在于分解函数的周期性（三角函数）。这些函数可在频域定位，但在时域则不行。因此，傅立叶变换是表示稳定信号的理想工具。小波变换的分解函数，即所谓的“小波”在时域和频域都具有局部性。因此，小波变换是表示包含不连续或不稳定成分的信号的理想工具。在实际的信号处理应用中，常采用的是小波变换的离散形式，也即所谓的离散小波变换（DWT）。1988年，Mallat在Burt和Adelson的图像分解和重构的塔形算法启示下，基于多分辨分析框架，给出了信号与图像（或函数）分解为不同频率通道的算法及其重构算法Mallat算法1，大大提高了小波变换的计算速度，其重要性相当于傅里叶变换中的FFT算法。从而使小波变换具有了明显的工程应用价值。通用的微处理器在运算速度上难以适应信号实时处理的要求。DSP处理器就是为了适应这种需求而开发的。DSP处理器中集成有高速的乘法器硬件，能快速地进行大量数据的乘法和加法运算。小波变换和DSP处理器的结合使得小波变换具有很强的可实现性。1快速小波变换（FWT）基于多分辨分析的理论，Mallat给出了快速小波变换（FWT）的算法Mallat算法。对于一个多分辨分析VjjZ，和分别为相应的尺度函数和小波函数，函数f属于由尺度函数生成的MRA的基本子空间V0。它可以表述成5对上述表达式的几点说明如下：（1）hk是多分辨分析的尺度系数，gk是多分辨分析的小波系数，两者的关系是gn（1）nh2k1n，kZ；（2）dmn和cmn分别是函数f在第m层尺度上的小波分解系数和尺度逼近系数；（3）实际应用中，函数f不是连续函数，而是所要分析的信号的离散采样序列fk。当采样率相当高时，信号的采样就非常近似于展开系数c0。这样，计算离散小波变换可以描述成递归计算式（1）和式（2）两个和式。因此，分析算法有如下形式：输入c0（输入序列）M（分解层次）计算for m1，M2有限序列的FWT以上算法推导时，假定序列是无限的。然而，实际应用时序列总是有限的。为了满足多分辨分析的要求，可以用合适的数来扩展有限序列，常采用的有零插值法和周期性扩展法。零插值法是将序列支撑以外的元素都设为零。如果序列不是以零开头或结尾，这就在序列的边界处人为地产生了不连续性，从而使得在dm相应位置出现较大的值。另外，零插值法中，输入序列与系数序列的支撑的长度不是相互独立的，这增加了应用时的计算量。周期性扩展法是把序列当作一个周期信号进行扩展，序列支撑的长度就是周期的长度。周期性扩展克服了零插值法的不足，而且这种周期性能沿着尺度传递。另外，可以利用DSP处理器的循环寻址功能实现序列的周期性扩展，从而更加提高了FWT的计算速度。因此，本文采用周期性扩展法在DSP处理器上实现FWT。对式（1，2）变换下标可得 输入序列定义为x（n），尺度展开系数为z（k），而系数为f（l）（包括尺度系数h（l）和小波系数g（l）。当算法应用于DSPs上时，这些序列的下标必须从0开始。系数序列的下标的范围为：0，lupper，lupper表示下标的最大值 如果输入序列x（n）的长度N2u，uZ，则z（k）的长度沿尺度不断地被二等分；因此，下标k的范围是：0，N/2-1。 这样，就可以通过Mallat算法计算某个尺度值：将系数序列放至输入序列的开头，各个系数与相应的输入值相乘再累加。最后所得就是尺度展开序列上的第一个值。然后，系数序列向前跳过输入序列的两个位置（基于2的向下抽样），再计算各个系数与相应的输入值的乘积的和，这样就得到第二个值。依次类推，直到系数序列跳出输入序列的范围。此时为了使每个系数都与输入序列中的值相乘，把输入序列的开始部分的一些值加到序列的结尾。整个计算完成后，得到的序列就是下一个尺度计算的输入序列。算法的原理如图1所示。图中长度为4的系数序列f（l）（l0，3）平移过长度为N输入序列x（n）（n0，N1；N2u，uZ）。当k0，计算尺度“m1”上的第一个值 xm1（0）xm1（0）f（0）xm（0）f（1）xm（1） f（2）xm（2）f（3）xm（3）系数序列然后平移至“k1”，计算第二个值xm1（1）xm1（1）f（0）xm（2）f（1）xm（3） f（2）xm（4）f（3）xm（5）如此反复，直到移至“kN21”。这时，只有头两个系数可以与输入序列相乘。为此，将xm（0）和xm（1）加到序列的结尾处。这样，尺度“m1”的最后一个值xm1（N21）就可以由下式计算xm+1（N/21）=f（0）xm(N-2)+ f（1）xm（N1） f（2）xm（0） f（3）xm（1）然后，把xm1（n）（n0，N21）作为输入值计算下一个尺度“m2”。3循环寻址大多数的DSPs的指令系统中有“循环寻址”（Circular addressing）的指令。周期性扩展可以通过循环寻址很容易地实现。循环寻址是间接编址的特殊情况，它用到一个循环缓存区。这种循环缓存区通常用于卷积、相关、数字滤波等算法中。它可以看成一个滑动的窗口，窗口内是最新的待处理的数据，如果有新的数据输入，就覆盖掉原有的数据，从而形成循环缓存。 数据x（n）的开头由一个指针定位，x（n）之前的数据从这个指针开始按顺时针方向连续写入。当得到一个新数据后，其将被放置到x（n）。在有关的算法计算中（比如FIR滤波器的应用），通过指针平移来连续定位每一个数据，直到最后一个输出数据的得出。这时，指针又回到了x（n）。然后，指针再逆时针方向移动一个位置，指向x（nN）。因为这个值是“最旧的”，不再有用，下一个数据就写到x（nN）的位置。周而复始，直至整个计算完成。在文2中讨论了循环缓存区在数字滤波器中的应用。所有的微处理器（包括DSP处理器）都提供线性存储空间，而循环寻址使得线性存储空间可以用作循环缓存空间。本文中，开始地址和循环缓存区的大小必须预先设定。指针指向开始地址，并且每次循环执行完成后随即增加“1”以指向下一个较高的存储位置，直到指针指向缓存尾部。循环缓存的这种结构很好地满足了输入信号的周期性扩展的要求，进而可用于实现FWT算法。上面提到为了在序列的边界上，每个系数都有输入值与之相乘，将序列开始部分的一些值添加到序列的末尾。如果输入序列应用循环缓存区，那么在到达数据缓冲末尾时，指针又一次跳到开始部分。也就是说，添加过程是自动完成的。循环寻址的原理在文3中阐述。 4小波变换在DSP中的应用对DSP编写代码主要有两种方法：（1）用高级语言编写，常用的如C语言；（2）用汇编语言编写。两者有各自的优缺点，见表1。大多数的DSP程序在一小段的代码上往往占去极大多数的执行时间。根据经验法则，90的执行时间将被10的代码占据。为了充分利用高级语言和汇编语言的优点，最好的方法是编写一个高效的汇编语言的子程序，在C语言程序中象普通函数一样调用它4。本文以TMS320C3X系列为例来阐述小波变换在DSP中的应用。Texas Instruments的TMS320C3X系列在相应的用户手册上有详细的说明。汇编语言子程序的主要部分是计算一个尺度上的尺度展开值，以下是针对TMS320C3X DSP的子程序的主要部分： LDIXARRAY，AR0LDI ZARRAY，AR1LDI WLETTAB，AR2LDI 0，IR0LDI DX，BKLDI DZ，RCSUBI 1，RCRPTB WTBLOCKLDF 00，R0LDF 00，R2PUSH STPUSH RSPUSH REPUSH RCRPTS LUPPERMPYF AR0，AR2，R0ADDF R0，R2ADDF R0，R2POP RCPOP REPOP RSPOP STSTF R2，AR1ADDI 2，IR0LDI XARRAY，AR0NOP AR0（IR0）WTBLOCKLDIWLETTAB，AR2尺度计算是在两个嵌套的循环中完成的。在内循环中，每次通过系数与输入值相乘累加求得一个尺度值（见式（6）。外循环保证每次内循环执行完毕后系数序列能平移到输入序列的正确位置（基于2向下抽样）。外循环运行的次数对应于该尺度下所要计算的尺度值的个数。比如输入序列的长度是N（偶数），那么尺度值个数就是N2。输入值、尺度值和小波系数都由指针来定位。这些指针必须在程序的开头初始化 LDIXARRAY，AR0LDIZARRAY，AR1LDIWLETTAB，AR2指针AR0指向输入序列，AR1指向要计算的尺度值的数据缓存区，AR2指向所用的小波系数表。外循环通过循环WTBLOCK代码块来计算某个尺度所有的尺度值（RPTB指令，Repeat block of code）。决定循环执行次数的循环次数寄存器RC必须赋予所要分析的尺度的长度（DZ），最小是1。LDIDZ，RCSUBI 1，RC内循环通过指令RPTS（Repeat single instruction）来计算其中一个尺度值。RPTSsrc指令循环执行下一条指令src1次。在程序中，循环次数与所用小波系数序列长度相等（RPTSLUPPER）。尺度的任何一个值正是通过循环执行下一条指令计算得出：把输入序列与小波系数相乘，再累加，这些乘法和加法由指令的平行操作结构来实现。在这种结构中，两个指令是同时执行的。这种结构具有高度的并发性。平行操作的标识符是“”，放在作平行操作的两条指令的第二条指令前。关于平行操作指令的详细说明可参考文3。相应的代码是RPTSLUPPERMPYF AR0，AR2，R0ADDF R0，R2 ADDF R0，R2在执行指令RPTS前，寄存器R0和R2必须置零（LDF00，R0和LDF00，R2）。指针AR0指向的输入值与指针AR2指向的小波系数相乘，结果存入寄存器R0。此外，R0中的值加到寄存器R2中。紧接着，指针AR0和AR2各增加1以指向下一个更高的存储位置。指针AR0是按循环寻址递增的（以“”标识），这样，当到达数据缓存区的末尾时，AR0在下一次操作中将重新指向缓存的开头。为了保证循环缓存区的正常工作，其大小必须在第一次执行外循环前定义。为此，寄存器BK被赋予输入序列的长度值（DZ）LDIDX，BK如果平行操作指令循环lupper次，寄存器R0中值反复地被加到寄存器R2中，最终在R2中得到该尺度下的一个尺度值。TMS320C3X执行指令采用流水线技术3。一条指令的执行经过四个独立的步骤：取指、译码、访问数据、执行。流水线技术是将各指令的执行时间重叠起来，使得每一条指令的最终执行时间是在单个指令周期内完成的，大大提高了速度。所有的寄存器都在执行段开始时读出，结束时写入。也就是说在上述例程第一次执行平行操作指令时，在把输入值与小波系数的乘积写入寄存器R0前，R0中的值（0）先被读出并加到R2中。第一次执行完成后，R2中的值仍是0，并不是第一次乘法操作的结果。因此，RPTS指令执行完成后，R0中的值必须再加入R2中一次。第二条ADDF指令累加了最后一次的乘积。紧接着，计算得出的尺度值（R2中的内容）存储到由指针AR1指向的存储位置中。AR1指向存储尺度值的存储空间。AR1在执行下面的指令后增加1STFR2，AR1在循环块WTBLOCK的结尾，指针AR0重新指向输入序列的开头，其内容增加IR0中的值LDIXARRAY，AR0NOPAR0（IR0）寄存器IR0在第一次执行循环模块WTBLOCK前初始化为0，在模块的结尾其值增加2。第一次执行时，AR0指向X（0），第二次指向X（2），第三次指向X（4），等。这样就实现了Mallat算法中的基于2的向下抽样。 指向小波系数的指针AR2在每次执行新的循环前必须指向小波的第一个系数LDIWLETTAB，AR2当执行嵌套循环操作时（在指令RPTB中执行指令RPTS），必须注意这样一个事实：指令RPTB和指令RPTS都利用了寄存器RS（Repeat start address register），RE（Repeat end address register）和RC（Repeat count register）。另外，两者都影响到寄存器RM（Repeat mode），它决定处理器是否工作在循环模式下RM0标准指令获取模式RM1循环模式因此，有必要在执行RPTS指令前将寄存器RS，RE，RC和ST的内容存储起来（指令PUSH），并在RPTS指令后将它们恢复（指令POP）。这样就避免了循环指令间的相互影响。有一点必须指出：寄存器存储和恢复的顺序对程序的正确运行非常重要。状态寄存器ST必须在寄存器RC，RE和RC之后恢复。指令RPTS完成后，循环计数寄存器RC的内容为0，表示内循环执行完成。此时，RM位也是0。POPST指令从堆栈中恢复执行RPTS指令前的状态寄存器的内容。如果此时寄存器RC0，POP