范安東 李小偉 王 娜 肖思和

（成都理工大學(xué) 數(shù)學(xué)地質(zhì)四川省重點實驗室，成都610059）

通過將2個N點實序列x（n）和y（n）組合成一個N 點復(fù)序列z（n）＝x（n）＋jy（n）的方法就可以同時計算出它們各自對應(yīng)的離散傅里葉變換（discrete Fourier transform，DFT），并能夠節(jié)約計算 量 和 存 儲 空 間［1，2］。1999 年 和 2002 年，Moshe、Hertz［3］和 Gunther［4］又先后提出了同時計算實序列的DFT和逆離散傅里葉變換（inverse discrete Fourier transform，IDFT）的直接公式。Moshe和Hertz的算法還曾被推廣到二維實序列和有 限 域 Fq2 上 的 Mersenne變 換 上 去［5，6］；而Gunther提出的算法也被推廣到二維的情況［7］，并應(yīng)用于數(shù)字水印嵌入取得了較好的效果［8］。

由文獻(xiàn)［4］可知，對于給定的X（k）和y（n），通過計算復(fù)序列z（n）的DFT（記為Z（k）），x（n）（X（k）對應(yīng)的IDFT）和Y（k）（y（n）對應(yīng)的 DFT）可由表1中的任一行給出的公式直接求出。

這里X＊（n）表示X（n）的共軛，（－n）N表示－n模N的最小非負(fù)剩余。在文獻(xiàn)［8］中，范安東等指出了表1中第4組公式的錯誤和文獻(xiàn)［4］中同時計算2個N點實序列的DFT公式中的錯誤，并進(jìn)行了改正，給出了相關(guān)正確公式的嚴(yán)格證明（文獻(xiàn)［4］中只給出了表1中的結(jié)果，未提供任何證明）。

表1 同時計算N點實序列的DFT和N點實序列的DFT的IDFT的新公式Table 1 New formulas for simultaneously calculating DFT and IDFT of two real N－point sequences

1 DFT的一些性質(zhì)

本節(jié)介紹DFT的一些相關(guān)性質(zhì)，特別是實序列的 DFT 的一些性質(zhì)［4，8－10］。

性質(zhì)1設(shè)x（n）是一個N 點實序列，則X（k）是Hermitian對稱的，即

這里k＝0，1，…，N－1。

注1：如無特殊說明，本文中變量n和k的取值范圍均表示為：n，k＝0，1，…，N－1。

性質(zhì)2設(shè)x（n）是一個N 點實序列，x（n）的DFT為X（k），X（n）的 DFT為則

性質(zhì)3設(shè)a（n）是一個N點復(fù)序列，則a（n）可以分解為下面4個N點序列之和

性質(zhì)4

1）x（n）是實偶序列?X（k）是實偶序列。

2）x（n）是實奇序列?X（k）是純虛奇序列。

性質(zhì)5設(shè)x（n）是一個N 點實序列，把X（k）分解為形如（3）式，則

引理［8］（同時計算2個實序列的DFT的一個新公式）

設(shè)x（n）和y（n）是2個N 點實序列，且令z（n）＝x（n）＋jy（n），則

其中

這里X（k）、Y（k）和Z（k）分別是x（n）、y（n）和z（n）的DFT。

2 同時計算y（n）的DFT和X（k）的IDFT的新公式及其在卷積計算中的應(yīng)用

2.1 新公式

本節(jié)假定x（n）和y（n）是2個N 點實序列，X（k）是x（n）的DFT，y（n）和X（k）是給定的，要求通過計算一個N點復(fù)序列z（n）的DFT（設(shè)為Z（k））同時求出x（n）和Y（k）的值。由性質(zhì)3和5可知，y（n）可由（n）和（n）決定，X（k）可由（k）和（k）決定。

利用上面討論的DFT的一些性質(zhì)和引理不難得到下面的結(jié)果。

定理1（同時計算實序列的DFT和IDFT的新公式）設(shè)

如果

則

注2：本定理改正了文獻(xiàn)［4］表1中第2行所給出的公式，更正后公式的正確性可用上節(jié)給出的引理和DFT的一些性質(zhì)進(jìn)行證明，詳細(xì)證明請參閱文獻(xiàn)［12］。

2.2 傳統(tǒng)的卷積計算方法

循環(huán)卷積是數(shù)字信號處理中經(jīng)常用到的一個重要運算。由于循環(huán)卷積和離散傅里葉變換在數(shù)學(xué)定義形式上存在一定的內(nèi)在聯(lián)系，通過分析易知計算2個信號的卷積可通過計算2個DFT和1個IDFT來完成，而DFT和IDFT的計算又存在快速算法（FFT和IFFT），故在實際應(yīng)用中一般都是通過DFT來計算信號的卷積。

定理2［2，11］記長度為 N 的2個序列x（n）和h（n）（其中n＝0，1，…，N－1），它們的循環(huán)卷積為

這里（n－k）N表示n－k模N 的最小非負(fù)剩余。

記x（n）的 DFT 為 X（k），h（n）的 DFT 為H（k），Y（k）的IDFT為y（n），則有

上述定理稱為循環(huán)卷積定理。該定理說明只要先分別計算出x（n）和h（n）的DFT，然后將它們相乘即可得到x（n）和h（n）的卷積的DFT，最后再做一次IDFT即可得到x（n）和h（n）的卷積。

利用定理2，可以得到傳統(tǒng)的計算卷積的過程：

1）計算x（n）和h（n）的DFT

2）計算X（k）和 H（k）乘積

3）計算Y（k）的IDFT

從上面的計算過程可以看出求2個長為N的實序列的卷積需要計算2個一維N點DFT和1個一維N點IDFT。據(jù)此推算，當(dāng)計算m個長為N的實序列所構(gòu)成的信號簇通過一個濾波器的輸出時，使用傳統(tǒng)的方法計算此m個卷積需要計算m＋1個一維N點DFT和m個一維N點IDFT。

2.3 新公式在計算多個卷積中的應(yīng)用

從上一小節(jié)中可以看到，當(dāng)計算由m個長為N的序列所構(gòu)成的信號簇通過一個濾波器的輸出時，使用傳統(tǒng)的計算方法其計算量還是比較大的。本小節(jié)將引入同時計算實序列的DFT與IDFT的新公式，得到新的計算一維卷積的方法。

假設(shè)一簇信號｛xi（n），（i＝0，1，…，m－1）｝通過一個濾波器h（n），濾波器的輸出分別為｛yi（n），（i＝0，1，…，m－1）｝。根據(jù)引理和定理1可以得到新的計算m個信號組成的信號簇通過一個濾波器的一維循環(huán)卷積算法。新算法的流程如圖1所示。

圖1 m個一維卷積運算的新算法流程圖The flow chart of a new algorithm for computing the convolution of a m－serial

首先設(shè)信號簇｛xi（n）｝的 DFT為｛Xi（k）｝，濾波器h（n）的 DFT 為 H（k），濾波器的輸出｛yi（n）｝的 DFT為｛Yi（k）｝，則新的計算卷積的過程如下：

第1步設(shè)z（n）＝x0（n）＋jh（n），利用引理，可同時得到x0（n）和h（n）的DFT：

再計算

注3：H（k）在后面的計算過程中還要用到，由于濾波器的脈沖響應(yīng)總是相對固定的，因此可以在進(jìn)行正式的濾波處理之前將H（k）先計算出來備用。

第2步不妨記x（n）＝x1（n），Y（k）＝Y(jié)0（k），y（n）＝y(tǒng)0（n），X（k）＝X1（k），設(shè)

利用定理1可以得到

從而同時得到了y0（n）和X1（k），再計算

第3～（m－1）步可仿照第2步進(jìn)行。

第m步不妨記x（n）＝xm－1（n），Y（k）＝Y(jié)m－2（k），y（n）＝y(tǒng)m－2（n），X（k）＝Xm－1（k）。根據(jù)定理1和用第2步類似的方法構(gòu)造序列z（n），然后利用新公式計算可同時得到y(tǒng)m－2（n）和Xm－1（k），再計算

第m＋1步利用IDFT計算可得

新算法到此結(jié)束。

顯然，如果用上面給出的多個卷積計算的新方法來計算長度為N的m個信號通過單個濾波器的輸出時，僅僅需要m個長度為N 的DFT和1個長度為N的IDFT。與傳統(tǒng)的算法相比，減少了m－1個IDFT的計算，從而提高了計算效率。關(guān)于二維的情況與本文討論的一維的情形類似，詳細(xì)內(nèi)容可參閱文獻(xiàn)［12］。

2.4 傳統(tǒng)和改進(jìn)卷積計算方法對比分析

本數(shù)值模擬分析是在Intel P4 2.8G＋512M＋ WinXP SP2＋Visual C＋＋6.0環(huán)境下進(jìn)行的，模擬數(shù)據(jù)分別為128個長度為256的實序列x［128］［256］與1個長度為256的實序列h［256］。分別采用2.2中的傳統(tǒng)卷積計算方法和2.3中的新方法來計算128個長為256的實序列分別與1個長為256的實序列的卷積，即相當(dāng)于計算128個信號通過1個濾波器的輸出。用傳統(tǒng)和改進(jìn)方法進(jìn)行計算的運行時間對比如表2所示。

表2 傳統(tǒng)與改進(jìn)方法運行時間對比Table 2 Contrast of the running time of the traditional and improved methods

由于軟件的運行時間會受到很多因素的影響，表2中的結(jié)果是將程序反復(fù)運行多次后出現(xiàn)的1個較穩(wěn)定的數(shù)據(jù)。從表2所示的運行時間對比來看，改進(jìn)方法的運行時間大約比傳統(tǒng)方法節(jié)約1／3，這也充分證明了將同時計算實序列的DFT和IDFT新公式應(yīng)用于卷積計算確實有較大的優(yōu)勢。

3 結(jié)論

對于給定的X（k）和y（n）（這里y（n）是一個實序列，X（k）是另一個實序列x（n）的 DFT），Gunther（2002）給出的4組同時計算x（n）和Y（k）的直接公式中第2、4組都有錯誤，作者在文獻(xiàn)［8］更正了第4組公式的基礎(chǔ)上，改正了第2組公式中的錯誤，并以此公式為基礎(chǔ)提出了一種快速計算多個卷積的新算法。

用一個復(fù)序列的DFT同時計算一個實序列的DFT和另一個實序列的DFT的IDFT的工程意義為有可能將傳統(tǒng)數(shù)字信號處理中的編碼和解碼過程同時進(jìn)行，從而極大地提高數(shù)字信號處理的速度。本文第2部分討論了計算m個長為N的序列所構(gòu)成的信號簇通過一個濾波器的輸出問題，這是數(shù)字信號處理中的典型應(yīng)用。分析結(jié)果表明，新方法能較大地減少計算量，提高計算效率。

［1］Smith W W，Smith J M.Handbook of Real－Time Fast Fourier Transform［M］.Piscataway NJ：IEEE Press，1995.

［2］蔣長錦，蔣勇.快速傅里葉交換及其C程序［M］.合肥：中國科學(xué)技術(shù)大學(xué)出版社，2004.

［3］Moshe S，Hertz D.On computing DFT of real N－point vector and IDFT of DFT－transformed real N－point vector via single DFT ［J］.IEEE Signal Processing Letters，1999，6（6）：141.

［4］Gunther Jacob H.Simultaneous DFT and IDFT of real N－point sequences ［J］.IEEE Signal Processing Letters，2002，9（8）：245－246.

［5］Sun Q，Tang Y Y.Some new explicit formulas for computing 2－D DFT of a real matrix and IDFT of 2－D DFT of another real matrix by a single 2－D DFT and their applications［J］.Journal of Sichuan University：Natural Science Edition，2000，37（6）：819－828.

［6］Liu L，Sun Q.A problem on mersenne transform over Fq2［J］.Advances in Mathematics，2004，33（4）：500－501.

［7］范安東，孫琦.同時計算二維實序列的DFT和二維實序列的DFT的IDFT的新公式［J］.成都理工大學(xué)學(xué)報：自然科學(xué)版，2007，34（6）：691－696.

［8］范安東，孫琦.同時計算實序列的DFT和IDFT的新公式及其在數(shù)字音頻水印中的應(yīng)用［J］.四川大學(xué)學(xué)報：工程科學(xué)版，2008，40（2）：96－100.

［9］Gasquet C，Witomski P.Fourier Analysis and Applications［M］.北京：世界圖書出版公司，2005.

［10］胡廣書.數(shù)字信號處理［M］.北京：清華大學(xué)出版社，2003.

［11］Yuan H，Chen H F，Yao D Z.A new general linear convolution model for fMRI data process［J］.JESTC，2005，3（1）：68－71.

［12］范安東.離散傅里葉變換及其在密碼學(xué)中的應(yīng)用［D］.成都：四川大學(xué)檔案館，2008.