基于Fourier變換離散化的連續(xù)小波變換頻域算法

易 華

?

基于Fourier變換離散化的連續(xù)小波變換頻域算法

易 華

（井岡山大學(xué)數(shù)理學(xué)院，江西，吉安 343009）

對于固定的尺度, 小波變換是待分析信號與小波基函數(shù)的線性卷積。當(dāng)小波基函數(shù)的Fourier變換有顯式表達(dá)式時(shí), 利用其Fourier變換進(jìn)行線性卷積稱為小波變換的頻域計(jì)算方法。由于線性卷積的長度大于信號的長度, 因此, 選取線性卷積中的哪一部分作為小波變換的系數(shù)也是一個(gè)亟需回答的問題。本文利用Fourier變換的離散化和離散Fourier變換的關(guān)系由小波變換時(shí)域算法推導(dǎo)了小波變換頻域算法，證明了時(shí)域算法與頻域算法的等價(jià)性; 解釋了這兩種方法分別應(yīng)該選取線性卷積中的哪一部分作為小波變換的系數(shù); 分析了頻域算法產(chǎn)生邊界效應(yīng)的原因; 給出了頻域算法中參數(shù)的選取方法, 以便克服邊界效應(yīng)。時(shí)間復(fù)雜度分析以及數(shù)值實(shí)驗(yàn)均表明了頻域算法至少比時(shí)域算法減少了1/3的運(yùn)行時(shí)間。

連續(xù)小波變換；Fourier變換；離散時(shí)間Fourier變換；離散Fourier變換；線性卷積；周期卷積

0　引言

連續(xù)小波變換的算法大致分為兩類。第一類快速算法依賴于雙尺度方程，得到近似的低通、高通濾波器系數(shù)，從而小波變換系數(shù)可以迭代求解[1]。這類算法的時(shí)間復(fù)雜度為Ο(N)[1-2]。但是，此類算法限制了尺度的取值范圍，例如限制尺度取二進(jìn)尺度或者整數(shù)值。而且近似的濾波器系數(shù)使得其在某些應(yīng)用中，計(jì)算精度不夠高。

第二類快速算法直接對連續(xù)小波變換的積分表達(dá)式離散化，因此要求小波函數(shù)的解析表達(dá)式是顯式給出的。例如Morlet小波、Mexican hat小波等。此類方法的代表有直接數(shù)值積分法，基于Mellin變換的快速算法[3-5]，基于調(diào)頻Z變換的方法[6]等。

由于小波變換可以寫為信號與小波函數(shù)的線性卷積，線性卷積在一定條件下等價(jià)于循環(huán)卷積，循環(huán)卷積又可以轉(zhuǎn)化為頻域的乘法。而時(shí)頻域的轉(zhuǎn)換可以采用FFT、IFFT快速算法。因此，第二類連續(xù)小波變換算法的另外一個(gè)代表為基于FFT的快速算法，其時(shí)間復(fù)雜度為Ο(N log(N))?；贔FT的快速算法又可以分為時(shí)域算法與頻域算法。如果直接對小波函數(shù)在時(shí)域進(jìn)行采樣，再利用線性卷積函數(shù)與信號進(jìn)行卷積，此方法稱為時(shí)域算法。如果信號與小波函數(shù)進(jìn)行卷積時(shí)，直接利用小波函數(shù)的頻域表達(dá)式進(jìn)行采樣，再與信號的頻域表達(dá)式進(jìn)行點(diǎn)點(diǎn)相乘，再進(jìn)行IFFT，則稱此方法為頻域算法。

頻域算法的精度較高，且這類算法被廣泛用于各種開源軟件中。例如，Torrence[7]采用連續(xù)小波變換方法分析了ENSO數(shù)據(jù)，并在http://paos. colorado. edu/research/wavelets/給出了計(jì)算連續(xù)小波變換的源代碼。其他的包含連續(xù)小波變換算法的開源軟件有JLAB(http://www. jmlilly.net)[8]，以及Wavelab (http://www-stat.stanford.edu/~wavelab)[9-10]。

但是，頻域算法應(yīng)該注意處理邊界效應(yīng)。雖然文獻(xiàn)[7]給出了處理邊界效應(yīng)的一種方法，但是，根據(jù)筆者了解的文獻(xiàn)來看，邊界效應(yīng)產(chǎn)生的原因未有報(bào)道。

Matlab小波工具箱也有實(shí)現(xiàn)連續(xù)小波變換算法的函數(shù)“cwt”，該函數(shù)是時(shí)域算法的典型代表。但是，其計(jì)算精度低。文獻(xiàn)[11]摒除了“cwt”函數(shù)中冗余且?guī)碚`差的差分運(yùn)算，計(jì)算精度有所提高。但是，文獻(xiàn)[11]中的連續(xù)小波變換方法以及“cwt”函數(shù)均涉及到小波基函數(shù)的某種近似。文獻(xiàn)[12]進(jìn)一步改進(jìn)了文獻(xiàn)[11]的方法，計(jì)算精度較文獻(xiàn)[11]有顯著提高。

文獻(xiàn)[12]中的數(shù)值實(shí)驗(yàn)結(jié)果表明，時(shí)域方法和頻域方法的計(jì)算精度幾乎完全一樣，揭示了這兩種方法的某種等價(jià)性，但是相應(yīng)的等價(jià)性理論還不夠完備。

本文在第1節(jié)介紹了基本的定義、定理以及之前的工作。在第2節(jié)，定理2.1揭示了時(shí)域算法、頻域算法的等價(jià)關(guān)系。同時(shí)也說明了時(shí)域算法、頻域算法分別應(yīng)該選取線性卷積的哪些系數(shù)作為小波變換的結(jié)果。定理2.2指出了頻域算法克服邊界效應(yīng)的方法。定理2.3對比了時(shí)域算法、頻域算法的時(shí)間復(fù)雜度。第3節(jié)給出的數(shù)值實(shí)驗(yàn)驗(yàn)證了上述理論結(jié)果。第4節(jié)為本文的小結(jié)。

1 基本的定義、定理以及之前的工作

連續(xù)小波變換可以寫為卷積的形式[9]

由(4)式離散化，得到了連續(xù)小波變換的時(shí)域算法。

線性卷積的計(jì)算可以根據(jù)線性卷積的定義進(jìn)行計(jì)算，也可以采用基于FFT的快速算法。

則

利用Fourier Parseval公式，由(1)定義的連續(xù)小波變換可以寫為頻域積分的形式：

這里

(8)

附注1.4 該算法中有一個(gè)參數(shù)需要預(yù)先給定，算法的推導(dǎo)未給出選取的方法。文獻(xiàn)[10]將選為信號的長度，該方法的好處在于信號的頻譜可以直接由信號的離散Fourier變換得到，同時(shí)離散Fourier反變換后得到的結(jié)果長度也為。只是，當(dāng)數(shù)據(jù)本身不具有周期性質(zhì)時(shí)，該方法得到的連續(xù)小波變換具有邊界效應(yīng)。

在本文的下一小節(jié)，將由連續(xù)小波變換時(shí)域算法推導(dǎo)頻域算法。

2 本文的工作

另外一方面，

所以，(14)得到了證明。

引理2.3表明小波變換時(shí)域算法可以由（15）式進(jìn)行計(jì)算。

則

從而

進(jìn)一步，

由定理1.2，可知

則

從而

對比時(shí)域算法(15)式與頻域算法(29)式，可以知道時(shí)域算法與頻域算法的濾波器的離散Fourier變換分別為

下文的數(shù)值實(shí)驗(yàn)也表明該選取方法可以避免邊界效應(yīng)的產(chǎn)生。

3 數(shù)值實(shí)驗(yàn)

圖1 Matlab產(chǎn)生的長度為1024的信號以及時(shí)頻圖(三個(gè)時(shí)頻圖分別采用了時(shí)域算法、Wavelab中的算法[10]以及頻域算法)

圖1中的(b，c，d)是對(a)中信號進(jìn)行連續(xù)小波變換后的時(shí)頻圖。所采用的基本小波為Morlet小波

圖2為頻域算法與時(shí)域算法得到的小波系數(shù)模值的誤差圖。這里最大的誤差為3.2872e-15。

圖2 圖1中(d)與(b)的誤差，即頻域算法與時(shí)域算法小波系數(shù)的誤差。這里最大的誤差為3.2872e-15。

為了比較時(shí)域算法與頻域算法的運(yùn)行時(shí)間，需要將這兩種算法在相同的運(yùn)行環(huán)境下獨(dú)立運(yùn)行多次。再取每次運(yùn)行的平均時(shí)間進(jìn)行比較。

表1 采用頻域算法與時(shí)域算法運(yùn)行圖1中的程序所耗費(fèi)的時(shí)間對比

4 小結(jié)

[1] Leigh G M. Fast FIR Algorithms for the Continuous Wavelet Transform From Constrained Least Squares[J]. IEEE Transactions on Signal Processing, 2013, 61(1):28-37.

[2] Vrhel M J, Lee C, Unser M. Fast continuous wavelet transform: A least-squares formulation[J]. Signal Processing, 1997,57(2):103-119.

[3] Alotta G, Paola M D, Failla G. A Mellin transform approach to wavelet analysis[J]. Communications in Nonlinear Science & Numerical Simulation, 2015,28(1-3):175-193.

[4] 張彤,楊福生. 基于Mellin變換的連續(xù)小波變換快速算法[J]. 信號處理, 1996(4):342-349.

[5] 王東偉,應(yīng)懷樵. 一種基于Mellin變換的連續(xù)小波變換快速算法[J]. 信號處理, 1999(3):278-280.

[6] Jones D L, Baraniuk R G. Efficient approximation of continuous wavelet transforms[J]. Electronics Letters, 1991,27(9):748-750.

[7] Torrence C, Compo G P. A Practical Guide to Wavelet Analysis[J]. Bulletin of the American Meteorological Society, 1998,79(79):61-78.

[8] Lilly J M, Olhede S C. Higher-Order Properties of Analytic Wavelets[J]. IEEE Transactions on Signal Processing, 2009,57(1):146-160.

[9] Mallat S. A wavelet tour of signal processing the sparse way[M]. Academic:Academic press, 2008.

[10] Buckheit J B, Donoho D L. Wavelab and reproducible research[M]. Springer:Springer press, 1995.

[11] 趙元英,袁曉,魏永豪.序列信號連續(xù)子波變換的MATLAB方法實(shí)現(xiàn)[J].四川大學(xué)學(xué)報(bào):自然科學(xué)版, 2006,43(2):325-329.

[12] Yi H, Chen Z, Cao Y. High Precision Computation of Morlet Wavelet Transform for Multi-period Analysis of Climate Data[J]. Journal of Information and Computational Science, 2014,11(17):6369-6385.

[13] Yi H, Shu H. The improvement of the Morlet wavelet for multi-period analysis of climate data[J]. Comptes Rendus Geoscience, 2012,344(10):483-497.

[14] Mallat S G. A Wavelet Tour of Signal Processing[J]. The Sparse Way, 2009,31(3):83-85.

[15] Proakis J G. Digital signal processing: principles, algorithms, and application-3/E[M]. Prentice: Prentice- Hall of India,1996.

[16] Frazier M W. An Introduction to Wavelets Through Linear Algebra[M]. Springer:springer press, 1999.

The Frequency Domain Algorithm of Continuous Wavelet Transform Based on the Discretization of Fourier Transform

YI Hua

（School of Mathematics and Physics, Jinggangshan University, Ji’an, Jiangxi 343009, China）

For a fixed scale, continuous wavelet transform(CWT) is a linear convolution of signal and wavelet function. Because the length of a linear convolution is larger than that of a signal, which part of linear convolution should be chosen as the wavelet coefficients is a question that need an answer. If the Fourier transform of wavelet function has explicit expression, the expression can be exploited to calculate the linear convolution, which method is known as the frequency-domain algorithm of CWT. In this paper, the frequency-domain algorithm of CWT is deduced from the time-domain algorithm of CWT by exploiting the relationship of the discretization of Fourier transform and the discrete Fourier transform. Four results are given. Firstly, the equivalent relation of these two algorithms are proved; Secondly, for time-domain method, the middle coefficients of linear convolution, while for frequency-domain method, the first coefficients of linear convolution are equal, and are the right wavelet coefficients. Thirdly, how to choose the parameter of the frequency domain algorithm is given by the theoretical derivation, which can conquer the boundary effect of frequency-domain method. In the end, the time complexity of the frequency domain algorithm of CWT is lower than that of the time domain algorithm, which is demonstrated by theoretical analysis and numerical experiments.

continuous wavelet transform; Fourier transform; discrete Fourier transform; linear convolution; circular convolution

1674-8085(2019)03-0001-08

O244

A

10.3969/j.issn.1674-8085.2019.03.001

2019-01-18；

2019-03-02

江西省自然科學(xué)基金項(xiàng)目(20161BAB201017)

易 華(1973-)，男, 湖北松滋人，講師，博士，主要從事小波分析及其應(yīng)用研究(E-mail:876145777@qq.com).