快速傅氏变换离小波变换.docVIP

10 快速傅氏变换和离散小波变换 长期以来，快速傅氏变换(Fast Fourier Transform)和离散小波变换(Discrete Wavelet Transform)在数字信号处理、石油勘探、地震预报、医学断层诊断、编码理论、量子物理及概率论等领域中都得到了广泛的应用。各种快速傅氏变换(FFT)和离散小波变换(DWT)算法不断出现，成为数值代数方面最活跃的一个研究领域，而其意义远远超过了算法研究的范围，进而为诸多科技领域的研究打开了一个崭新的局面。本章分别对FFT和DWT的基本算法作了简单介绍，若需在此方面做进一步研究，可参考文献[2]。 快速傅里叶变换FFT 离散傅里叶变换是20世纪60年代是计算复杂性研究的主要里程碑之一，1965年Cooley和Tukey所研究的计算离散傅里叶变换(Discrete Fourier Test)的快速傅氏变换(FFT)将计算量从О(n2)下降至О(nlogn)，推进了FFT更深层、更广法的研究与应用。FFT算法有很多版本，但大体上可分为两类：迭代法和递归法，本节仅讨论迭代法，递归法可参见文献[1]、[2]。 串行FFT迭代算法 假定a[0],a[1], …,a[n-1] 为一个有限长的输入序列，b[0], b[1], …,b[n-1]为离散傅里叶变换的结果序列，则有：，其中 W，实际上，上式可写成矩阵W和向量a的乘积形式： 一般的n阶矩阵和n维向量相乘，计算时间复杂度为n2，但由于W是一种特殊矩阵，故可以降低计算量。FFT的计算流图如 REF _Re\h \* MERGEFORMAT 图 22.1所示，其串行算法如下： 算法22.1 单处理器上的FFT迭代算法 输入：a=(a0,a1, …,an-1) 输出：b=(b0,b1, …,bn-1) Begin (1)for k=0 to n-1 do ck=ak end for (2)for h=logn-1 downto 0 do (2.1) p=2h (2.2) q=n/p (2.3) z=wq/2 (2.4) for k=0 to n-1 do if (k mod p=k mod2p) then (i)ck = ck + ck +p (ii)ck +p=( ck - ck +p)z k modp /* ck不用(i)计算的新值 */ end if end for end for (3)for k=1 to n-1 do br(k) = ck /* r(k)为k的位反 */ end for End 图 STYLEREF 1 \s 22. SEQ 图 \* ARABIC \s 1 1 n=4时的FFT蝶式变换图 显然, FFT算法的计算复杂度为O(nlogn)。 并行FFT算法 设P为处理器的个数，一种并行FFT实现时是将n维向量a划分成p个连续的m维子向量，这里，第i个子向量中下标为i×m, …, (i+1)×m-1，其元素被分配至第i号处理器（i=0,1, …, p-1）。由 REF _Re\h 图 22.1可以看到，FFT算法由logn步构成，依次以2logn-1、2logn-2、…、2、1为下标跨度做蝶式计算，我们称下标跨度为2h的计算为第h步（h=logn-1, logn-2, …, 1, 0）。并行计算可分两阶段执行：第一阶段，第logn-1步至第logm步，由于下标跨度h≥ m，各处理器之间需要通信；第二阶段，第logm-1步至第0步各处理器之间不需要通信。具体并行算法框架描述如下： 算法22.2 FFT并行算法 输入：a=(a0,a1, …,an-1) 输出：b=(b0,b1, …,bn-1) Begin 对所有处理器my_rank(my_rank=0,…, p-1)同时执行如下的算法: (1)for h=logp-1 downto 0 do /* 第一阶段，第logn-1步至第logm步各处理器之间需要通信*/ (1.1) t=2i, ,l=2(i+logm) ,q=n/l , z=wq/2 , j= j+1 ,v=2j /*开始j=0*/ (1.2)if ((my_rank mod t)=(my_rank mod 2t)) then /*本处理器的数据作为变换的前项数据*/ (i) tt= my_rank+p/v (ii)接收由tt 号处理器发来的数据块，并将接收的数据块存于 c[my_rank*m+n/v]到c[my_rank*m+n/v+m]之中 (iii)for k=0 to m-1 do c[k]=c[k]+c[k+n/v] c[k+n/v]=( c[k]- c[k+n/v])*z(my_rank*m+k