《快速傅氏变换和离散小波变换》由会员分享，可在线阅读，更多相关《快速傅氏变换和离散小波变换（13页珍藏版）》请在金锄头文库上搜索。
1、1. 10 快速傅氏变换和快速傅氏变换和离散小波变换离散小波变换长期以来，快速傅氏变换快速傅氏变换(Fast Fourier Transform)和离散小波变换离散小波变换(Discrete Wavelet Transform)在数字信号处理、石油勘探、地震预报、医学断层诊断、编码理论、量子物理及 概率论等领域中都得到了广泛的应用。各种快速傅氏变换(FFT)和离散小波变换(DWT)算法 不断出现，成为数值代数方面最活跃的一个研究领域，而其意义远远超过了算法研究的范 围，进而为诸多科技领域的研究打开了一个崭新的局面。本章分别对 FFT 和 DWT 的基本 算法作了简单介绍，若需在此方面做进一步研
2、究，可参考文献2。 1.1 快速傅里叶变换快速傅里叶变换 FFT离散傅里叶变换是 20 世纪 60 年代是计算复杂性研究的主要里程碑之一，1965 年 Cooley 和 Tukey 所研究的计算离散傅里叶变换离散傅里叶变换(Discrete Fourier Test)的快速傅氏变换(FFT) 将计算量从 (n2)下降至 (nlogn)，推进了 FFT 更深层、更广法的研究与应用。FFT 算法 有很多版本，但大体上可分为两类：迭代法和递归法，本节仅讨论迭代法，递归法可参见 文献1、2。1.1.1 串行串行 FFT 迭代算法迭代算法假定 a0,a1, ,an-1 为一个有限长的输入序列，b0, b
3、1, ,bn-1为离散傅里叶变换的结果序列，则有：，其中 W，实际上，上式可) 1,.,2 , 1 , 0(10nkWkakbnmkmnnine2 写成矩阵 W 和向量 a 的乘积形式：1210) 1)(1() 1(2) 1(0) 1(2420121000001210nnnnnnnnaaaawwwwwwwwwwwwwwwwbbbbM LMOMMMLLLM一般的 n 阶矩阵和 n 维向量相乘，计算时间复杂度为 n2，但由于 W 是一种特殊矩阵， 故可以降低计算量。FFT 的计算流图如图 22.1 所示，其串行算法如下： 算法算法 22.1 单处理器上的单处理器上的 FFT 迭代算法迭代算法 输入
4、：输入：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 t
5、o n-1 do br(k) = ck /* r(k)为 k 的位反 */ end forEnda0a1a2a3a4a5a6a7b0b4b2b6b1b5b3b7000000221230图图 1.1 n=4 时的时的 FFT 蝶式变换图蝶式变换图显然, FFT 算法的计算复杂度为 O(nlogn)。1.1.2 并行并行 FFT 算法算法设 P 为处理器的个数，一种并行 FFT 实现时是将 n 维向量 a 划分成 p 个连续的 m 维子向量，这里，第 i 个子向量中下标为 im, , (i+1)m-1，其元素被分配至第 ipnm/号处理器（i=0,1, , p-1） 。由图 1.1 可以看到，FF
6、T 算法由 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=
7、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 号处理器发来的数据块，并将接收的数据块存于 cmy_rank*m+n/v到 cmy_rank*m+n
8、/v+m之中 (iii)for k=0 to m-1 do ck=ck+ck+n/v ck+n/v=( ck- ck+n/v)*z(my_rank*m+k) mod l end for (iv)将存于 cmy_rank*m+n/v到 cmy_rank*m+n/v+m之中的数据块发送 到 tt 号处理器 else /*本处理器的数据作为变换的后项数据*/ (v)将存于之中的数据块发送到 my_rank-p/v 号处理器 (vi)接收由 my_rank-p/v 号处理器发来的数据块存于 c 中 end if end for (2)for i=logm-1 downto 0 do /*第二阶段，第
9、logm-1 步至第 0 步各处理器之间 不需要通信*/ (2.1) l=2i ,q=n/l , z=wq/2 (2.2)for k=0 to m-1 do if (k mod l)=(k mod 2l) then ck=ck+ck+l ck+l=( ck- ck+l)*z(my_rank*m+k) mod l end if end for end for End由于各处理器对其局部存储器中的 m 维子向量做变换计算，计算时间为；点对nmlog点通信次，每次通信量为 m，通信时间为，因此快速傅里叶变换的plog2pmttwslog)(2并行计算时间为。pmttnmTwsplog)(2logMP
10、I 源程序请参见章末附录。1.2 离散小波变换离散小波变换 DWT1.2.1 离散小波变换离散小波变换 DWT 及其串行算法及其串行算法先对一维小波变换作一简单介绍。设 f(x)为一维输入信号，记，这里与分别称为定标函)2(2)(2/kxxjjjk)2(2)(2/kxxjjjk)(x)(x数与子波函数，与为二个正交基函数的集合。记 P0f=f，在第级上的一)(xjk)(xjkj维离散小波变换离散小波变换 DWT(Discrete Wavelet Transform)通过正交投影 Pjf 与 Qjf 将 Pj-1f 分解为： kkjkj kjkj kjjjdcfQfPfP1其中：， ，这里，h(
11、n)与 101 2)(pnj nkj kcnhc 101 2)(pnj nkj kcngd) 12,.,1 , 0,.,2 , 1(jNkLjg(n)分别为低通与高通权系数，它们由基函数与来确定，p 为权系数)(xjk)(xjk的长度。为信号的输入数据，N 为输入信号的长度，L 为所需的级数。由上式可见，0 nC每级一维 DWT 与一维卷积计算很相似。所不同的是：在 DWT 中，输出数据下标增加 1 时， 权系数在输入数据的对应点下标增加 2，这称为“间隔取样” 。算法算法 22.3 一维离散小波变换串行算法一维离散小波变换串行算法 输入：输入：c0=d0(c00, c10, cN-10)h=
12、(h0, h1, hL-1)g=(g0, g1, gL-1) 输出：输出：cij , dij (i=0, 1, N/2j-1, j0) Begin (1)j=0, n=N(2)While (n1) do (2.1)for i=0 to n-1 do (2.1.1)cij+1=0, dij+1=0 (2.1.2)for k=0 to L-1 doj nikkj ij ijni) (kkj ij idgdd ch ccmod)2(11 mod211, end for end for (2.2)j=j+1, n=n/2end while End显然，算法 22.3 的时间复杂度为 O(N*L)。 在
13、实际应用中，很多情况下采用紧支集小波紧支集小波（Compactly Supported Wavelets） ，这时相 应的尺度系数和小波系数都是有限长度的，不失一般性设尺度系数只有有限个非零值：h1,hN，N 为偶数，同样取小波使其只有有限个非零值：g1,gN。为简单起见，设尺度系数与小波函数都是实数。对有限长度的输入数据序列：(其余点的值都Mnxcnn, 2 , 1,0L看成 0)，它的离散小波变换为:kn Znj nj khcc21 kn Znj nj kgcd21 1, 1 , 0JjL其中 J 为实际中要求分解的步数，最多不超过 log2M，其逆变换为knZkj kknZkj kjnh
14、chcc221 1 ,LJj 注意到尺度系数和输入系列都是有限长度的序列，上述和实际上都只有有限项。若完全按照上述公式计算，在经过 J 步分解后，所得到的 J+1 个序列和1, 1 , 0,Jjdj kLj kc的非零项的个数之和一般要大于 M，究竟这个项目增加到了多少？下面来分析一下上述计算过程。 j=0 时计算过程为knMnnkhxc211 knMnnkgxd2 11 不难看出，1 kc的非零值范围为：即有, 12, 0 , 1, 12 MNkLL个非零值。1 kd的非零值范围相同。继续往下分解时，非零项 21 212NMMNk出现的规律相似。分解多步后非零项的个数可能比输入序列的长度增加
15、较多。例如，若输入序列长度为 100，N=4，则1 kd有 51 项非零，2 kd有 27 项非零，3 kd有 15 项非零，4 kd有9 项非零，5 kd有 6 项非零，6 kd有 4 项非零，6 kc有 4 项非零。这样分解到 6 步后得到的序列的非零项个数的总和为 116，超过了输入序列的长度。在数据压缩等应用中，希望总的 长度基本不增加，这样可以提高压缩比、减少存储量并减少实现的难度。 可以采用稍微改变计算公式的方法，使输出序列的非零项总和基本上和输入序列的非 零项数相等，并且可以完全重构。这种方法也相当于把输入序列进行延长（增加非零项） ， 因而称为延拓法。 只需考虑一步分解的情形，下面考虑第一步分解(j=1)。将输入序列作延拓，若 M 为偶数，直接将其按 M 为周期延拓；若 M 为奇数，首先令01Mx。然后按 M+1 为周期延拓。作了这种延拓后再按前述公式计算，相应的变换矩阵已不再是 H 和 G，事实上这时的变换 矩阵类似于循环矩阵。例如，当 M=8，N=4 时矩阵 H 变为：2143432
