連續(xù)子波變換的MATLAB實(shí)現(xiàn)分析

1、連續(xù)子波變換的MATLAB實(shí)現(xiàn)分析趙元英 趙元英，女，四川大學(xué)電子信息學(xué)院碩士生，研究方向：信號(hào)與信息處理。E-mail: zyyflyworld.，袁 曉，滕旭東，魏永豪（四川大學(xué)電子信息學(xué)院，成都 610064）摘 要：本文目的在于考察序列信號(hào)的連續(xù)子波變換（CWT）在MATLAB中的程序?qū)崿F(xiàn)。首先介紹子波變換的基本概念和MATLAB中的算法原理，對(duì)MATLAB中的cwt.m程序進(jìn)行深入解讀，闡述了各個(gè)參數(shù)的功能；而后發(fā)現(xiàn)差分運(yùn)算的存在使子波濾波器系數(shù)含有誤差，繼而對(duì)此誤差的產(chǎn)生原因進(jìn)行分析，解釋如何調(diào)整cwt.m算法中的參數(shù)，以克服差分運(yùn)算的影響，改善子波濾波器系數(shù)的取值，并且提出一種新

2、算法，改善了CWT變換的運(yùn)算效果；最后以頻移鍵控FSK信號(hào)為例，分別用兩種算法對(duì)之進(jìn)行連續(xù)子波變換，并對(duì)結(jié)果進(jìn)行分析和比較，證明了本文的改進(jìn)算法優(yōu)于MATLAB中的原算法。關(guān)鍵詞：子波濾波器；尺度；morlet子波；FSK信號(hào)1. 引言MATLAB是當(dāng)今世界上應(yīng)用最廣泛的數(shù)學(xué)軟件之一，其最引人入勝的功能是具有各領(lǐng)域?qū)Ｓ玫墓ぞ呦洌@些工具箱的出現(xiàn)大大增強(qiáng)了MATLAB的實(shí)用性。子波分析是在現(xiàn)代調(diào)和分析的基礎(chǔ)上發(fā)展起來(lái)的一門(mén)學(xué)科，具有深刻的理論和廣泛的應(yīng)用1。本文研究了子波工具箱中的連續(xù)子波變換(CWT)的實(shí)現(xiàn)，詳細(xì)解讀了MATLAB中的cwt.m程序，并對(duì)該算法中由于引入差分運(yùn)算而導(dǎo)致的誤差進(jìn)行

3、了解析，從而提出一種新算法，改善了CWT的運(yùn)算效果。2. MATLAB中的CWT算法2.1 算法原理及程序?qū)崿F(xiàn)子波變換是一種多分辨分析方法2，對(duì)于平方可積函數(shù)：，其連續(xù)子波變換定義為 (1)式中和分別稱(chēng)為尺度參數(shù)和平移參數(shù)。在MATLAB中，基于子波變換的定義式，為便于數(shù)值計(jì)算，將信號(hào)和子波函數(shù)都進(jìn)行了離散化。令，則有式（2）。由此可見(jiàn)，算法的本質(zhì)是信號(hào)與子波濾波器函數(shù)積分的差分進(jìn)行卷積，其中的差分序列就是子波變換時(shí)采用的子波濾波器系數(shù)，算法的優(yōu)劣取決于此系數(shù)的選取是否精確。 （2）在MATLAB中,基于上述算法思想，利用cwt命令對(duì)信號(hào)進(jìn)行連續(xù)子波變換。cwt命令的格式為function c

4、oefs = cwt(signal, scales, wname, plotmode, xlim)，signal是輸入信號(hào)，scales為子波變換的尺度范圍，wname是子波函數(shù)的名稱(chēng)，plotmode是系數(shù)coefs的繪圖模式，xlim用來(lái)設(shè)置橫坐標(biāo)的范圍，coefs即為對(duì)信號(hào)連續(xù)子波變換后的系數(shù)值。通過(guò)查找toolbox/wavelet/wavelet/cwt.m文件，得到其主干程序如下：precis = 10;/*設(shè)置子波樣本點(diǎn)的二進(jìn)階數(shù)，確定子波序列積分值PSI_INTEG的點(diǎn)數(shù)為個(gè)。*/signal = signal(:); /* 輸入信號(hào) */len = length(signal

5、); /* 求信號(hào)序列的長(zhǎng)度 */coefs = zeros(length(scales),len); /* 開(kāi)辟變換系數(shù)矩陣 */nbscales = length(scales); /* 求尺度序列的長(zhǎng)度 */psi_integ,xval=intwave (wname,precis); /* 將子波函數(shù)從到XVAL進(jìn)行積分，所得積分的返回值為PSI_INTEG。*/xval = xval-xval(1); /* 對(duì)子波的時(shí)間變量t定起始點(diǎn)為0: xval(1)=0 */dx = xval(2); /* 求取母波時(shí)間變量t的初始采樣間隔Tsdx =xval(2) */xmax = xval(

6、end); /* 求取母波時(shí)間變量t的終點(diǎn):max(t)=xval(end) */ind = 1; /* 子波變換系數(shù)矩陣的行指標(biāo) */for k = 1:nbscalesa = scales(k);j = 1+floor(0:a*xmax/(a*dx); /* 求取在尺度a下子波y(n/a)的整數(shù)采樣點(diǎn), 一維矢量j給定了y(n/a)在psi_integ中的位置,即可獲取在尺度a下的子波濾波器系數(shù)序列: psi_integ(j). */if length(j) = 1 , j = 1 1; endf = fliplr(psi_integ(j); /* 子波濾波器系數(shù)序列倒置,即完成y(-n/

7、a),然后賦值于f=y(-n/a) */coefs(ind,:)=-sqrt(a)*wkeep(diff(conv(signal,f), len); /* 利用conv函數(shù)將信號(hào)與子波濾波器系數(shù)f進(jìn)行卷積，再用diff函數(shù)進(jìn)行差分，其原理如式（2）所示，之后wkeep截取中間數(shù)值,忽略掉兩邊極小的值。乘以進(jìn)行能量歸一化。 */ind = ind+1;end2. 2 算法的優(yōu)點(diǎn)和缺陷MATLAB算法中，通過(guò)尺度a調(diào)節(jié)矢量j，對(duì)子波濾波器系數(shù)進(jìn)行采樣，對(duì)于不同的尺度a，濾波器系數(shù)的數(shù)據(jù)長(zhǎng)度不同，因此可以達(dá)到不同的分辨精度，從而實(shí)現(xiàn)了子波變換的多分辨分析思想。通過(guò)改變采樣精度來(lái)控制分辨率，這是MAT

8、LAB算法中一個(gè)重要的思想。但是，算法中存在的差分運(yùn)算，直接導(dǎo)致了誤差的產(chǎn)生。如式（2）所示，信號(hào)實(shí)質(zhì)上是與子波濾波器函數(shù)積分的差分進(jìn)行卷積。對(duì)子波濾波器函數(shù)的積分psi_integ采樣后得到的是一個(gè)反對(duì)稱(chēng)序列，通過(guò)對(duì)進(jìn)行差分，得到的序列即為連續(xù)子波變換中使用的子波濾波器系數(shù)，序列是一個(gè)對(duì)稱(chēng)序列，如圖1所示（以下均使用morlet子波）。差分運(yùn)算通過(guò)MATLAB中的diff函數(shù)完成，格式為Y = diff(X)，用來(lái)計(jì)算序列X中相鄰元素的差值，YX(2)-X(1) X(3)-X(2) . X(n)-X(n-1)。通過(guò)序列中相鄰項(xiàng)依次相減，得到子波濾波器系數(shù)序列，由于對(duì)子波濾波器函數(shù)的積分psi

9、_integ的均勻采樣而得到，使得序列中數(shù)值變化小的部分采樣密集，數(shù)值變化大的部分采樣稀疏，在差分操作后，由于在采樣密度變化較大處，相鄰數(shù)值間的差值相對(duì)較大，因此，子波濾波器系數(shù)序列在這些點(diǎn)處產(chǎn)生尖狀突起，引起波形畸變，如圖1所示。使用這樣的子波濾波器系數(shù)對(duì)信號(hào)進(jìn)行分析，勢(shì)必使分析結(jié)果存在誤差，因此，有必要調(diào)節(jié)算法中的參數(shù)以產(chǎn)生準(zhǔn)確的子波濾波器系數(shù)。3. 算法的改進(jìn)3.1 參數(shù)調(diào)節(jié)圖1中差分運(yùn)算引起的子波波形畸變，其本質(zhì)原因是子波濾波器函數(shù)的積分序列的取值不均勻，才使得差分之后序列的波形不平滑。因此改善的采樣數(shù)據(jù)，可以改善的波形，提高子波變換的分辨率。子波濾波器函數(shù)的積分序列g(shù)(k)對(duì)g(k)

10、差分后得到的子波濾波器系數(shù)圖1 MATLAB中通過(guò)差分運(yùn)算產(chǎn)生子波濾波器系數(shù)（尺度a30）在cwt.m程序中，precis變量是控制子波樣本點(diǎn)的采樣精度的，precis的取值不同，子波濾波器函數(shù)的積分psi_integ的數(shù)據(jù)長(zhǎng)度不同，psi_integ的個(gè)數(shù)等于。當(dāng)precis的值較大時(shí)，psi_integ的數(shù)據(jù)個(gè)數(shù)比較多，從而通過(guò)對(duì)psi_integ采樣得到，使得序列中相鄰兩個(gè)樣本點(diǎn)之間的距離分布比較均勻，繼而在對(duì)差分后，得到的子波濾波器系數(shù)序列相對(duì)平滑。圖2中顯示了在不同的precis取值下，所得到的子波濾波器系數(shù)。由圖2可見(jiàn)，增大采樣精度precis，可以改善的波形，因此，使用MATLA

11、B進(jìn)行連續(xù)子波變換時(shí)，不適合直接采用cwt命令對(duì)信號(hào)進(jìn)行分析，應(yīng)當(dāng)深入cwt.m程序，按不同的精度要求，對(duì)cwt.m程序進(jìn)行局部修改，才能得到精確的分析結(jié)果。對(duì)g(k)差分后得到的子波濾波器系數(shù)h(k),precis10對(duì)g(k)差分后得到的子波濾波器系數(shù)h(k),precis15圖2 不同precis取值下，差分運(yùn)算產(chǎn)生子波濾波器系數(shù)（尺度a30）3.2新算法及其程序?qū)崿F(xiàn)為避免差分引入的誤差，可以摒除差分過(guò)程，而依照子波變換的定義式(1)，將信號(hào)與子波濾波器函數(shù)直接卷積以求得其子波變換系數(shù)。通過(guò)wavefun函數(shù)得到子波濾波器系數(shù)，將信號(hào)與卷積，得到連續(xù)子波變換。wavefun是子波和尺度函

12、數(shù)，它可以求得任意子波的子波函數(shù)和相應(yīng)的尺度函數(shù)（尺度函數(shù)存在情況）。改進(jìn)后的主干程序如下：signal = signal(:);len = length(signal);nbscales = length(scales); coefs = zeros(nbscales,len);for k = 1:nbscalesa = scales(k);wtype = wavemngr(type,wname);/* 子波管理函數(shù)，可以判斷子波的類(lèi)型 */if wtype = 4psi,xval = wavefun(wname); /* 改進(jìn)之一：利用wavefun函數(shù)產(chǎn)生子波尺度函數(shù)phi和子波函數(shù)ps

13、i。*/elsephi,psi,xval=wavefun(wname);/* 如果wtype=4,說(shuō)明此子波不具有尺度函數(shù)phi。 */endxval = xval-xval(1);dx = xval(2);xmax = xval(end);j = 1+floor(0:a*xmax/(a*dx); /* 求取尺度a下子波y(n/a)的整數(shù)采樣點(diǎn)個(gè)數(shù),即可獲取子波濾波器系數(shù)序列：psi_integ(j)，決定子波變換的精度。 */if length(j) = 1 , j = 1 1; endf = fliplr(psi(j); /* 子波濾波器系數(shù)序列倒置,即完成y(-n/a),然后賦值于f=y

14、(-n/a) */coefs(k,:)=wkeep(conv(signal,f),len)/sqrt(a); /* 改進(jìn)之二：將信號(hào)與子波濾波器系數(shù)f直接卷積得到其子波變換系數(shù)，除以進(jìn)行能量歸一化。 */ endMATLAB中的子波濾波器系數(shù)新算法中的子波濾波器系數(shù)圖3 子波濾波器系數(shù)比較（尺度a30）圖3對(duì)MATLAB中和新算法中的子波濾波器系數(shù)進(jìn)行了比較，在光滑性和準(zhǔn)確度方面，新算法得到的子波濾波器系數(shù)明顯優(yōu)于MATLAB中的濾波器系數(shù)。采用此種算法，同樣可以通過(guò)尺度a控制分辨精度，實(shí)現(xiàn)子波變換的多分辨分析思想；同時(shí)摒除了差分運(yùn)算，通過(guò)wavefun函數(shù)得到了準(zhǔn)確的子波濾波器系數(shù)，在對(duì)信號(hào)

15、進(jìn)行子波變換時(shí)，提高了變換的精度。4. MATLAB中原算法與改進(jìn)算法的比較上節(jié)中將MATLAB中的算法利用兩種方法進(jìn)行了改進(jìn)，一是調(diào)整參數(shù)，二是去除差分運(yùn)算，通過(guò)wavefun函數(shù)得到子波濾波器系數(shù)。這兩種方法都能得到準(zhǔn)確的子波濾波器系數(shù)，所得的效果一致。因此只將去除差分運(yùn)算的新算法與MATLAB中的原算法進(jìn)行比較。4.1 FSK信號(hào)的連續(xù)子波變換輸入一個(gè)頻移鍵控FSK信號(hào)， （3）為信號(hào)能量，i為載波頻率，為初相，為單位幅度的矩形函數(shù)，其支撐區(qū)間為，是信號(hào)的碼元周期。將4FSK信號(hào)進(jìn)行連續(xù)子波變換，使用以上兩種算法分別計(jì)算，使用的子波為morlet子波，其變換結(jié)果如圖4所示。信噪比為9dB

16、，分析尺度為1到65。圖4 4FSK信號(hào)的連續(xù)子波變換(a) 使用MATLAB中的算法所得的變換 (b) 使用改進(jìn)的算法所得的變換 圖4中亮暗程度顯示了系數(shù)模值的大小，顏色亮的表示系數(shù)模值較大。對(duì)于FSK信號(hào)而言，其載波頻率是變化的，不同的碼元時(shí)刻，存在不同的頻率，因此其子波系數(shù)圖表現(xiàn)為階梯狀，能量集中在不同的尺度，且階梯的個(gè)數(shù)與載頻數(shù)相同，根據(jù)子波變換尺度與頻率的對(duì)應(yīng)關(guān)系，可以推斷出信號(hào)的載頻大小。4.2 結(jié)果分析圖4中(a)和（b）采用的是不同的算法，由圖可見(jiàn)，兩種算法都有效地提取了信號(hào)的頻率和幅度特征，但新算法的變換結(jié)果更為清晰。在子波變換圖的上部，較大尺度的區(qū)域內(nèi)，明顯可以看到，圖像的

17、模糊程度得到了改善。圖5著重展示了尺度為22至62區(qū)域內(nèi)兩種算法所得到的結(jié)果，如圖所示，在清晰度方面，新算法的效果更好。這在本質(zhì)上得益于其子波濾波器系數(shù)的準(zhǔn)確提取。因此，所提出圖5 4FSK信號(hào)在尺度2262內(nèi)的連續(xù)子波變換(a) 使用MATLAB中的算法所得的變換(b) 使用改進(jìn)的算法所得的變換的新算法優(yōu)于MATLAB中的算法。5. 結(jié)束語(yǔ)本文解讀了MATLAB中的cwt.m程序，并對(duì)該算法中由于引入差分運(yùn)算而導(dǎo)致的誤差進(jìn)行了解析，深入分析了差分運(yùn)算如何影響各個(gè)參數(shù)的取值，從而一步步地使連續(xù)子波變換中關(guān)鍵的子波濾波器系數(shù)產(chǎn)生了畸變。隨后通過(guò)調(diào)節(jié)精度參數(shù)precis和運(yùn)用新算法這兩種途徑，得到

18、更為準(zhǔn)確的子波濾波器系數(shù)。最后通過(guò)使用新算法和原算法對(duì)FSK信號(hào)的處理，對(duì)兩種算法進(jìn)行了比較和分析。結(jié)果表明，通過(guò)對(duì)原算法的改進(jìn)，新算法提高了連續(xù)子波變換的變換精確度。參考文獻(xiàn)：1 Oliver Rioul, Martin Vetterli. Wavelets and signal processing. IEEE SP Magazine,1991,10, pp.14-38.2 張曉平，田立聲，彭應(yīng)寧.從時(shí)頻分布到連續(xù)子波變換.電子與信息學(xué)報(bào),1994.6，pp.631-640.3 楊福生.小波變換的工程分析與應(yīng)用. 北京：科學(xué)出版社，1999.Realization of Continuous Wavelet Transform by MATLABZHAO Yuanying, YUAN Xiao, TENG Xudong, WEI Yonghao(School of Electronic and Information, Sichuan University, Chengdu 610064, China)Abstract：In this paper, the program realization of continuous wavelet transform (CWT) by MATLAB is an