快速傅氏變換和離散小波變換

1、1. 10 快速傅氏變換和離散小波變換長期以來，快速傅氏變換(Fast Fourier Transform)和離散小波變換(Discrete Wavelet Transform)在數(shù)字信號處理、石油勘探、地震預(yù)報、醫(yī)學(xué)斷層診斷、編碼理論、量子物理及概率論等領(lǐng)域中都得到了廣泛的應(yīng)用。各種快速傅氏變換(FFT)和離散小波變換(DWT)算法不斷出現(xiàn)，成為數(shù)值代數(shù)方面最活躍的一個研究領(lǐng)域，而其意義遠(yuǎn)遠(yuǎn)超過了算法研究的范圍，進而為諸多科技領(lǐng)域的研究打開了一個嶄新的局面。本章分別對FFT和DWT的基本算法作了簡單介紹，若需在此方面做進一步研究，可參考文獻2。 1.1 快速傅里葉變換FFT離散傅里葉變換是2

2、0世紀(jì)60年代是計算復(fù)雜性研究的主要里程碑之一，1965年Cooley和Tukey所研究的計算離散傅里葉變換(Discrete Fourier Test)的快速傅氏變換(FFT)將計算量從(n2)下降至(nlogn)，推進了FFT更深層、更廣法的研究與應(yīng)用。FFT算法有很多版本，但大體上可分為兩類：迭代法和遞歸法，本節(jié)僅討論迭代法，遞歸法可參見文獻1、2。1.1.1 串行FFT迭代算法假定a0,a1, ,an-1 為一個有限長的輸入序列，b0, b1, ,bn-1為離散傅里葉變換的結(jié)果序列，則有：，其中 W，實際上，上式可寫成矩陣W和向量a的乘積形式：一般的n階矩陣和n維向量相乘，計算時間復(fù)雜

3、度為n2，但由于W是一種特殊矩陣，故可以降低計算量。FFT的計算流圖如圖 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=akend 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 doif (k mod p=k mod2p) then (i)ck = ck + ck +p(ii)ck +p=( ck - ck +p)z

4、 k modp /* ck不用(i)計算的新值 */end ifend forend for(3)for k=1 to n-1 do br(k) = ck /* r(k)為k的位反 */end forEnd圖 1.1 n=4時的FFT蝶式變換圖顯然, FFT算法的計算復(fù)雜度為O(nlogn)。1.1.2 并行FFT算法設(shè)P為處理器的個數(shù)，一種并行FFT實現(xiàn)時是將n維向量a劃分成p個連續(xù)的m維子向量，這里，第i個子向量中下標(biāo)為i×m, , (i+1)×m-1，其元素被分配至第i號處理器（i=0,1, , p-1）。由圖 1.1可以看到，F(xiàn)FT算法由logn步構(gòu)成，依次以2lo

5、gn-1、2logn-2、2、1為下標(biāo)跨度做蝶式計算，我們稱下標(biāo)跨度為2h的計算為第h步（h=logn-1, logn-2, , 1, 0）。并行計算可分兩階段執(zhí)行：第一階段，第logn-1步至第logm步，由于下標(biāo)跨度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)同時執(zhí)行如下的算法:(1)for h=logp-1 downto 0 do /* 第一階段，第

6、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 /*本處理器的數(shù)據(jù)作為變換的前項數(shù)據(jù)*/(i) tt= my_rank+p/v(ii)接收由tt 號處理器發(fā)來的數(shù)據(jù)塊，并將接收的數(shù)據(jù)塊存于cmy_rank*m+n/v到cmy_rank*m+n/v+m之中(iii)for k=0 to m-1 do ck=ck+ck+n/vck+n/v=( ck- ck+n/v)*z(my

7、_rank*m+k) mod lend for(iv)將存于cmy_rank*m+n/v到cmy_rank*m+n/v+m之中的數(shù)據(jù)塊發(fā)送到tt 號處理器else /*本處理器的數(shù)據(jù)作為變換的后項數(shù)據(jù)*/(v)將存于之中的數(shù)據(jù)塊發(fā)送到my_rank-p/v 號處理器(vi)接收由my_rank-p/v 號處理器發(fā)來的數(shù)據(jù)塊存于c中end if end for(2)for i=logm-1 downto 0 do /*第二階段，第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

8、l)=(k mod 2l) then ck=ck+ck+lck+l=( ck- ck+l)*z(my_rank*m+k) mod lend if end forend forEnd由于各處理器對其局部存儲器中的m維子向量做變換計算，計算時間為；點對點通信次，每次通信量為m，通信時間為，因此快速傅里葉變換的并行計算時間為。MPI源程序請參見章末附錄。1.2 離散小波變換DWT1.2.1 離散小波變換DWT及其串行算法先對一維小波變換作一簡單介紹。設(shè)f(x)為一維輸入信號，記，這里與分別稱為定標(biāo)函數(shù)與子波函數(shù)，與為二個正交基函數(shù)的集合。記P0f=f，在第級上的一維離散小波變換DWT(Discret

9、e Wavelet Transform)通過正交投影Pjf與Qjf將Pj-1f分解為： 其中：， ，這里，h(n)與g(n)分別為低通與高通權(quán)系數(shù)，它們由基函數(shù)與來確定，p為權(quán)系數(shù)的長度。為信號的輸入數(shù)據(jù)，N為輸入信號的長度，L為所需的級數(shù)。由上式可見，每級一維DWT與一維卷積計算很相似。所不同的是：在DWT中，輸出數(shù)據(jù)下標(biāo)增加1時，權(quán)系數(shù)在輸入數(shù)據(jù)的對應(yīng)點下標(biāo)增加2，這稱為“間隔取樣”。算法22.3 一維離散小波變換串行算法輸入：c0=d0(c00, c10, cN-10) h=(h0, h1, hL-1) g=(g0, g1, gL-1)輸出：cij , dij (i=0, 1, N/2j

10、-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 do end forend for(2.2)j=j+1, n=n/2 end whileEnd顯然，算法22.3的時間復(fù)雜度為O(N*L)。在實際應(yīng)用中，很多情況下采用緊支集小波（Compactly Supported Wavelets），這時相應(yīng)的尺度系數(shù)和小波系數(shù)都是有限長度的，不失一般性設(shè)尺度系數(shù)只有有限個非零值：h1,hN，N為偶數(shù)，同樣取小波使其只有有限個非零值：g1,gN。

11、為簡單起見，設(shè)尺度系數(shù)與小波函數(shù)都是實數(shù)。對有限長度的輸入數(shù)據(jù)序列：(其余點的值都看成0)，它的離散小波變換為:其中J為實際中要求分解的步數(shù)，最多不超過log2M，其逆變換為注意到尺度系數(shù)和輸入系列都是有限長度的序列，上述和實際上都只有有限項。若完全按照上述公式計算，在經(jīng)過J步分解后，所得到的J+1個序列和的非零項的個數(shù)之和一般要大于M，究竟這個項目增加到了多少？下面來分析一下上述計算過程。j=0時計算過程為 不難看出，的非零值范圍為：即有個非零值。的非零值范圍相同。繼續(xù)往下分解時，非零項出現(xiàn)的規(guī)律相似。分解多步后非零項的個數(shù)可能比輸入序列的長度增加較多。例如，若輸入序列長度為100，N=4，

12、則有51項非零，有27項非零，有15項非零，有9項非零，有6項非零，有4項非零，有4項非零。這樣分解到6步后得到的序列的非零項個數(shù)的總和為116，超過了輸入序列的長度。在數(shù)據(jù)壓縮等應(yīng)用中，希望總的長度基本不增加，這樣可以提高壓縮比、減少存儲量并減少實現(xiàn)的難度。可以采用稍微改變計算公式的方法，使輸出序列的非零項總和基本上和輸入序列的非零項數(shù)相等，并且可以完全重構(gòu)。這種方法也相當(dāng)于把輸入序列進行延長（增加非零項），因而稱為延拓法。只需考慮一步分解的情形，下面考慮第一步分解(j=1)。將輸入序列作延拓，若M為偶數(shù)，直接將其按M為周期延拓；若M為奇數(shù)，首先令。然后按M+1為周期延拓。作了這種延拓后再按

13、前述公式計算，相應(yīng)的變換矩陣已不再是H和G，事實上這時的變換矩陣類似于循環(huán)矩陣。例如，當(dāng)M=8，N=4時矩陣H變?yōu)椋寒?dāng)M=7，N=4時矩陣H變?yōu)椋簭纳鲜龅木仃嚤硎究梢钥闯?，兩種情況下的矩陣內(nèi)都有完全相同的行，這說明作了重復(fù)計算，因而從矩陣中去掉重復(fù)的那一行不會減少任何信息量，也就是說，這時我們可以對矩陣進行截短（即去掉一行），使得所得計算結(jié)果仍然可以完全恢復(fù)原輸入信號。當(dāng)M=8，N=4時截短后的矩陣為：當(dāng)M=7，N=4時截短后的矩陣為：這時的矩陣都只有行。分解過程成為：向量C1 和D1都只有個元素。重構(gòu)過程為：可以完全重構(gòu)。矩陣H，G有等式H*H+G*G=I一般情況下，按上述方式保留矩陣的行，

14、可以完全恢復(fù)原信號。這種方法的優(yōu)點是最后的序列的非0元素的個數(shù)基本上和輸入序列的非0元素個數(shù)相同，特別是若輸入序列長度為2的冪，則完全相同，而且可以完全重構(gòu)輸入信號。其代價是得到的變換系數(shù)Dj中的一些元素已不再是輸入序列的離散小波變換系數(shù)，對某些應(yīng)用可能是不適合的，但在數(shù)據(jù)壓縮等應(yīng)用領(lǐng)域，這種方法是可行的。1.2.2 離散小波變換并行算法下設(shè)輸入序列長度N=2t，不失一般性設(shè)尺度系數(shù)只有有限個非零值：h0，hL-1，L為偶數(shù)，同樣取小波使其只有有限個非零值：g0，gL-1。為簡單起見，我們采用的延拓方法計算。即將有限尺度的序列按周期N延長，使他成為無限長度的序列。這時變換公式也稱為周期小波變換

15、。變換公式為：其中<n+2k>表示n+2k對于模N/2j的最小非負(fù)剩余。注意這時和是周期為N/2j的周期序列。其逆變換為從變換公式中可以看出，計算輸出點和，需要輸入序列在n=2k附近的值（一般而言，L遠(yuǎn)遠(yuǎn)小于輸入序列的長度）。設(shè)處理器臺數(shù)為p，將輸入數(shù)據(jù)按塊分配給p臺處理器，處理器i得到數(shù)據(jù)，讓處理器i負(fù)責(zé)和的計算，則不難看出，處理器i 基本上只要用到局部數(shù)據(jù)，只有L/2個點的計算要用到處理器i+1中的數(shù)據(jù)，這時做一步并行數(shù)據(jù)發(fā)送：將處理器i+1中前L-1個數(shù)據(jù)發(fā)送給處理器i，則各處理器的計算不再需要數(shù)據(jù)交換，關(guān)于本算法其它描述可參見文獻1。算法22.4 離散小波變換并行算法輸入：

16、hi(i=0, L-1), gi(i=0, L-1), ci0(i=0, N-1) 輸出：cik (i=0, N/2k-1,k>0)Begin對所有處理器my_rank(my_rank=0, p-1)同時執(zhí)行如下的算法:(1)j=0;(2)while (j<r) do(2.1)將數(shù)據(jù)按塊分配給p臺處理器(2.2)將處理器i+1中前L-1個數(shù)據(jù)發(fā)送給處理器i(2.3)處理器i負(fù)責(zé)和的計算(2.4)j=j+1end whileEnd這里每一步分解后數(shù)據(jù)和已經(jīng)是按塊存儲在P臺處理器上，因此算法第一步中的數(shù)據(jù)分配除了j=0時需要數(shù)據(jù)傳送外，其余各步不需要數(shù)據(jù)傳送（數(shù)據(jù)已經(jīng)到位）。因此，按L

17、ogP模型，算法的總的通信時間為：2(Lmax(o,g)+l)，遠(yuǎn)小于計算時間O(N)。MPI源程序請參見所附光盤。1.3 小結(jié) 本章主要討論一維FFT和DWT的簡單串、并行算法，二維FFT和DWT在光學(xué)、地震以及圖象信號處理等方面起著重要的作用。限于篇幅，此處不再予以討論。同樣，F(xiàn)FT和DWT的并行算法的更為詳盡描述可參見文獻2和文獻3,專門介紹快速傅氏變換和卷積算法的著作可參見4。另外，二維小波變換的并行計算及相關(guān)算法可參見文獻5，LogP模型可參考文獻3。參考文獻 1. 王能超 著數(shù)值算法設(shè)計華中理工大學(xué)出版社,1988.92. 陳國良 編著并行計算結(jié)構(gòu)·算法·編程高

18、等教育出版社,1999.10 3. 陳國良 編著并行算法設(shè)計與分析（修訂版）高等教育出版社，2002.114. Nussbaumer H J. Fast Fourier Transform and Convolution Algorithms.2nded. Springer- Verlag,19825. 陳崚二維正交子波變換的VLSI并行計算電子學(xué)報,1995,23(02):95-97 附錄 FFT并行算法的MPI源程序1. 源程序fft.c#include <stdio.h>#include <stdlib.h>#include <complex.h>#i

19、nclude <math.h>#include "mpi.h"#define MAX_PROCESSOR_NUM 12#define MAX_N 50#define PI 3.#define EPS 10E-8#define V_TAG 99#define P_TAG 100#define Q_TAG 101#define R_TAG 102#define S_TAG 103#define S_TAG2 104void evaluate(complex<double>* f, int beginPos, int endPos, const compl

20、ex<double>* x, complex<double> *y, int leftPos,int rightPos, int totalLength);void shuffle(complex<double>* f, int beginPos,int endPos);void print(const complex<double>* f, int fLength);void readDoubleComplex(FILE *f,complex<double> &z); int main(int argc, char *arg

21、v)complex<double> pMAX_N, qMAX_N, s2*MAX_N, r2*MAX_N;complex<double> w2*MAX_N;complex<double> temp;int variableNum;MPI_Status status;int rank, size;int i, j, k, n;int wLength;int everageLength;int moreLength;int startPos;int stopPos;FILE *fin;MPI_Init(&argc, &argv); MPI_Get

22、_rank(MPI_COMM_WORLD, &rank); MPI_Get_size(MPI_COMM_WORLD, &size);if(rank = 0) fin = fopen("dataIn.txt", "r"); if (fin = NULL) puts("Not find input data file"); puts("Please create a file "dataIn.txt""); puts("<example for dataIn.txt&

23、gt; "); puts("2"); puts("1.0 2"); puts("2.0 -1"); exit(-1);readDoubleComplex(fin, variableNum); if (variableNum < 1)|(variableNum > MAX_N) puts("variableNum out of range!"); exit(-1); for(i = 0; i < variableNum; i +) readDoubleComplex(fin, pi);fo

24、r(i = 0; i < variableNum; i +) readDoubleComplex(fin, qi); fclose(fin); puts("Read from data file "dataIn.txt""); printf("p(t) = "); print(p, variableNum); printf("q(t) = "); print(q, variableNum); for(i = 1; i < size; i +) MPI_Send(&variableNum,1,MP

25、I_INT,i, V_TAG,MPI_COMM_WORLD); MPI_Send(p,variableNum, MPI_DOUBLE_COMPLEX,i,P_TAG, PI_COMM_WORLD); MPI_Send(q,variableNum,MPI_DOUBLE_COMPLEX,i, Q_TAG,MPI_COMM_WORLD);elseMPI_Recv(&variableNum,1,MPI_INT,0, V_TAG,MPI_COMM_WORLD, &status); MPI_Recv(p,variableNum,MPI_DOUBLE_COMPLEX,0, P_TAG, PI

26、_COMM_WORLD, &status); MPI_Recv(q,variableNum,MPI_DOUBLE_COMPLEX,0, Q_TAG,MPI_COMM_WORLD,&status);wLength = 2*variableNum;for(i = 0; i < wLength; i +)wi= complex<double>(cos(i*2*PI/wLength),sin(i*2*PI/wLength);everageLength = wLength / size;moreLength = wLength % size;startPos = mor

27、eLength + rank * everageLength;stopPos = startPos + everageLength - 1;if(rank = 0)startPos = 0;stopPos = moreLength+everageLength - 1;evaluate(p, 0, variableNum - 1, w, s,startPos, stopPos, wLength);evaluate(q, 0, variableNum - 1, w, r, startPos, stopPos, wLength);for(i = startPos; i <= stopPos ;

28、 i +) si = si*ri/(wLength*1.0); if (rank > 0)MPI_Send(s+startPos), everageLength,MPI_DOUBLE_COMPLEX, 0, S_TAG, MPI_COMM_WORLD); MPI_Recv(s,wLength, MPI_DOUBLE_COMPLEX,0, S_ TAG2,MPI_ COMM_WORLD,&status);else for(i = 1; i < size; i +) MPI_Recv(s+moreLength+i*everageLength),everageLength, MP

29、I_DOUBLE_COMPLEX, i,S_TAG, MPI_COMM_WORLD, &status);for(i = 1; i < size; i +) MPI_Send(s,wLength, MPI_DOUBLE_COMPLEX,i, S_TAG2,MPI_COM M_WORLD); for(int i = 1; i < wLength/2; i +)temp = wi;wi = wwLength - i; wwLength - i = temp;evaluate(s, 0, wLength - 1, w, r, startPos, stopPos, wLength);

30、if (rank > 0)MPI_Send(r+startPos), everageLength,MPI_DOUBLE_COMPLEX,0, R_TAG,MPI_COMM_WORLD);elsefor(i = 1; i < size; i +) MPI_Recv(r+moreLength+i*everageLength),everageLength, MPI_DOUBLE_COMPLEX,i, R_TAG,MPI_COMM_WORLD,&status);puts("nAfter FFT r(t)=p(t)q(t)"); printf("r(t)

31、 = "); print(r, wLength - 1); puts(""); printf("Use prossor size = %dn",size);MPI_Finalize();void evaluate(complex<double>* f, int beginPos, int endPos, const complex<double>* x,complex<double>* y, int leftPos, int rightPos, int totalLength)int i;complex<

32、;double> tempX2*MAX_N,tempY12*MAX_N,tempY22*MAX_N;int midPos = (beginPos + endPos)/2;if (beginPos > endPos)|(leftPos > rightPos)puts("Error in use Polynomial!"); exit(-1);else if(beginPos = endPos)for(i = leftPos; i <= rightPos; i +)yi = fbeginPos;else if(beginPos + 1 = endPos)

33、for(i = leftPos; i <= rightPos; i +) yi = fbeginPos + fendPos*xi;elseshuffle(f, beginPos, endPos); for(i = leftPos; i <= rightPos; i +) tempXi = xi*xi;evaluate(f, beginPos, midPos, tempX, tempY1, leftPos, rightPos,totalLength);evaluate(f, midPos+1, endPos, tempX, tempY2, leftPos, rightPos, tot

34、alLength);for(i = leftPos; i <= rightPos; i +)yi = tempY1i + xi*tempY2i;void shuffle(complex<double>* f, int beginPos, int endPos)complex<double> temp2*MAX_N;int i, j;for(i = beginPos; i <= endPos; i +)tempi = fi; j = beginPos;for(i = beginPos; i <= endPos; i +=2)fj = tempi; j +

35、;for(i = beginPos +1; i <= endPos; i += 2)fj = tempi; j +;void print(const complex<double>* f, int fLength)bool isPrint = false;int i;if (abs(f0.real() > EPS)printf(“%lf”, f0.real(); isPrint = true; for(i = 1; i < fLength; i +) if (fi.real() > EPS) if (isPrint) printf(" + &quo

36、t;); else isPrint = true; printf("%lft%d", fi.real(),i); else if (fi.real() < - EPS) if(isPrint) printf(" - "); else isPrint = true; printf("%lft%d", -fi.real(),i); if (isPrint = false) printf("0");printf("n");2. 運行實例編譯：mpicc o fft fft.cc 運行：使用如下命令運行程序mpirun np 1 fftmpirun np 2 fft