文档介绍：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的计算流图如图 ,其串行算法如下:
单处理器上的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
() p=2 h
() q=n/p
() z=wq/2
() 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
图 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)。由图 ,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步各处理器之间不需要通信。具体并行算法框架描述如下:
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步各处理器之间需要通信*/
() t=2i, ,l=2(i+logm) ,q=n/l , z=wq/2 , j= j+1 ,v=2 j /*开始j=0*/
()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) mod l
end for
(iv)将存于c[my_rank*m+n/v]到c[my_rank*m+n/v+m]之中的数据