非平穩(wěn)信號的廣義小波分析及其工程應用報告

1、非平穩(wěn)信號的廣義小波分析及其工程應用報告徐文豪1. 窗口傅里葉變換1.1 原理分析眾所周知，傅里葉變換可以獲取信號的全局頻譜，但很多時候，我們需要的是信號的瞬時頻譜，如機械故障檢測、地震信號瞬時屬性提取等。為了達到這個目的，一種很直觀的思路就是截取信號的一小段做傅里葉變換，并將得到的頻譜作為該小段中心位置的頻譜。這種變換稱之為窗口傅里葉變換（WFT），其示意圖如下：圖1 窗口福利葉變換示意圖從數(shù)學上描述WFT應為：對平方可積信號和平方可積窗函數(shù)，定義 窗口傅里葉變換如下：值得強調(diào)的是，式（1）中窗函數(shù)的支撐大小對時頻譜的有效性有著很大的影響，當窗寬過大或過小時都會使時頻譜出現(xiàn)較嚴重的假象。作者

2、在實際編程時發(fā)現(xiàn)，當信號的長度為，取一個較小的正值作為值零閾值，則取窗寬點數(shù)半徑為一般能達到較理想的效果。將視為整體，假設并通過逆傅里葉變換可得逆窗口傅里葉變換如下：式（2）對推導WFT值域的性質(zhì)有著很重要的作用，如重建核方程等。但由于使用雙重積分，其在實現(xiàn)時稍顯繁瑣，殷勤業(yè)的時頻分析及其在工程中的應用講義中給出了一個更簡單的重構(gòu)公式如下：其推導如下：1.2 程序?qū)崿F(xiàn)附錄1給出了與式（1）對應的WFT程序，附錄2給出了與式（3）對應的IWFT程序,對如下的四段調(diào)頻信號：使用WFT程序和IWFT程序，得到的時頻譜圖和誤差圖如下：圖2 四段跳頻信號WFT時頻譜圖及重構(gòu)誤差圖2. 連續(xù)小波變換2.1

3、 原理分析窗口傅里葉變換的缺點在于窗寬是固定的，因而時頻譜容易出現(xiàn)假頻現(xiàn)象，如下圖所示：圖3 WFT缺點示意圖圖3中，對快變信號采用大窗或?qū)β冃盘柌捎眯〈岸急厝粫斐杉兕l。解決這一問題的思路是使窗寬可變，這就是連續(xù)小波變換(CWT)的思路。在CWT中我們一般把窗函數(shù)稱為小波函數(shù)，其支撐如下所示：圖4 小波函數(shù)的支撐從數(shù)學上描述CWT應為：對平方可積信號和平方可積小波，定義連續(xù)小波變換如下：其中系數(shù)的引入是為了小波在平移伸縮之后范數(shù)不變，即對式（5）中的連續(xù)小波變換，常見的逆變換公式如下：其中稱為時間尺度原子。與WFT類似，式（7）對推導CWT值域的性質(zhì)有著很重要的作用，如重建核方程等，但其一

4、般只能通過二重數(shù)值積分進行計算，速度慢且誤差大。Holschneider在Wavelets: An Analysis Tool中給出了當范數(shù)變量時的一個簡單重構(gòu)公式：如果能求出的值，則可利用數(shù)值積分計算式（8）的值。遺憾的是式（9）是不一定收斂的，例如對常用的Morlet小波()，式（9）在0處發(fā)散。高老師的講義中給出了當范數(shù)時的一個簡單逆變換公式，但其有效性我還沒有搞懂，其計算過程如下：式（11）中表示單位脈沖信號在0時刻尺度為處CWT的值。2.2 程序?qū)崿F(xiàn)附錄3給出了與式（8）對應的CWT程序，附錄4給出了與式（10）對應的ICWT程序，對地球物理中常用的Ricker信號(取主頻為50Hz

5、)，取Gauss小波(，參考劉乃豪師兄提供的資料)使用CWT程序和ICWT程序，并利用頻率尺度轉(zhuǎn)換關系，得到其時頻譜和重構(gòu)誤差如下：圖5 50Hz主頻Ricker信號CWT時頻譜圖及重構(gòu)誤差圖從圖5可以看出，對Ricker信號，在選取與其相配的小波函數(shù)之后，不用考慮窗寬，連續(xù)小波變換就能得到很好的時頻譜和重構(gòu)精度。CWT具有自適應性的優(yōu)點可以從其Heisenberg盒得到進一步驗證：圖6 連續(xù)小波變換的Heisenberg盒從圖6可以看出，CWT既具有顯微鏡的功能，又具有望遠鏡的功能。即當尺度較小時，有可能捕捉小范圍的信號變化，而當尺度較大時，有可能捕捉大范圍的信號變化。然而，福兮禍之所倚，禍

6、兮福之所伏，圖6的Heisenberg盒也反映出了CWT的一大缺點，即低頻可能看不見(頻率分辨率過高，而離散尺度有限)，高頻可能看不清（頻率分辨率過低），這種現(xiàn)象可以從式（4）的四段調(diào)頻信號反應出來，同樣取高斯小波，四段調(diào)頻信號的CWT時頻譜和重構(gòu)誤差如下：圖7四段跳頻信號CWT時頻譜圖及重構(gòu)誤差圖從圖7可以很明顯的看出，在CWT時頻譜中，信號頻率分辨率隨著頻率的增大而降低。這里需要指出的是，緩解這一問題的一個有效方法是二進制采樣，其能夠在低頻部分采相對較多的點并在高頻部分采相對較少的點。令一方面，圖7也反應出圖5的高重構(gòu)精度并不適合所有信號，其原因可能有兩個，一是Ricker信號的變化較平緩

7、而四段調(diào)頻信號的變化較劇烈，二是高斯小波可能與四段調(diào)頻信號不太匹配。3. 最優(yōu)對偶標架變換3.1 原理分析標架是指Hilbert空間中的一個完備序列，其滿足對任意，存在，使得其中為的由內(nèi)積誘導的范數(shù)。和分別稱為標架下界和標架上界，特別地，當時，稱這個標架為緊標架，并稱為緊標架的冗余度。已經(jīng)證明對任意標架，存在對偶標架，使得對于任意，有而當為緊標架時，Daubechies已證明其對偶標架為。對采樣間隔為，長度為的周期離散信號(將有限離散信號周期化可以得到更好的邊界重構(gòu)精度)，定義其對偶標架變換及逆變換如下：其中為時間采樣點數(shù)，為波數(shù)采樣點數(shù)，為由分析函數(shù)的平移調(diào)制生成的標架，為由綜合函數(shù)的平移調(diào)

8、制生成的的對偶標架，的定義如下：其中分別為時間和波數(shù)初始采樣點（引入這兩個常量是為了在標架相空間中得到更多的空間波數(shù)信息），分別為時間和波數(shù)采樣間隔點數(shù)且滿足，為波數(shù)離散采樣間隔。通過給定綜合函數(shù)的離散采樣序列，由式（14）和式（15）可推出分析函數(shù)的離散采樣序列所滿足的方程組：其中，。式（17）所對應方程組的系數(shù)矩陣是行列的，故若且系數(shù)矩陣滿秩時可以解出唯一的，然而Daubechies已經(jīng)證明這種情況下重構(gòu)公式是不穩(wěn)定的，一般稱這種情形下的采樣為臨界采樣。錢世鍔在著作Joint Time Frequency Analysis中提出取， 并在方程組的解空間中尋找與最接近的解，令方程組的值向量為

9、，并假設則可構(gòu)造對應的優(yōu)化問題如下：對式（18）進行簡單推導得式（19）意味著要在的解空間中尋找模最小的解。由廣義逆矩陣理論，可給出所需解的表達式如下：其中表示矩陣的Moore-Penrose逆。對于計算，錢世鍔書中是通過判斷是否行滿秩并通過SVD分解將行不滿秩的情形轉(zhuǎn)換為行滿秩情形，這樣不僅計算繁瑣還會影響計算精度。實際上有直接計算矩陣MP逆的快速算法，如Greville遞推法，而在matlab中可直接調(diào)用庫函數(shù)pinv實現(xiàn)。然而，式（20）只能處理規(guī)模相對較小的矩陣，且由于使用復數(shù)運算會影響計算精度。錢世鍔書中在以上的基礎上通過將大系數(shù)矩陣分解為多個小系數(shù)矩陣，將復數(shù)運算化為實數(shù)運算，給出

10、了求解方程組的快速算法：其中，。這樣大復系數(shù)矩陣就變成了個小實系數(shù)矩陣，再利用類似推導式（20）的過程即可得到所需解。對采樣間隔為的信號，一般取如下歸一化的高斯函數(shù)作為綜合函數(shù)：其中的取值是錢世鍔書中給出的和間取得最小誤差的條件。必須注意到的是，的支撐不應也不必等于，否則當很大時，求解對偶標架將相對比較耗時。由殷勤業(yè)講義3-4節(jié)的內(nèi)容和作者編程時的經(jīng)驗，一般取時可以達到很好的效果，且若讓在的值為在的值并在其它點處為，則可以利用式（21）只計算長度為的對偶函數(shù)向量。對偶標架變換的重構(gòu)誤差精度與時頻采樣點數(shù)與信號長度之比有關，這可以被稱為對偶標架采樣冗余度，定義如下：3.2 程序?qū)崿F(xiàn)附錄5給出了與

11、式（21）對應的快速計算最優(yōu)對偶標架的程序，附錄6給出了與式（14）對應的最優(yōu)對偶標架變換的程序，附錄7給出了與式（15）對應的逆最優(yōu)對偶標架變換程序。同樣對式（4）的四段跳頻信號，取采樣冗余度，使用最優(yōu)對偶標架變換程序和逆最優(yōu)對偶標架變換程序得到時頻譜和重構(gòu)誤差如下：圖8四段跳頻信號最優(yōu)對偶標架時頻譜圖及重構(gòu)誤差圖從圖8可以看出，最優(yōu)對偶標架變換只用了WFT存儲量的，就得到了有效的譜圖和非常高的重構(gòu)精度，這暗示了最優(yōu)對偶標架變換在處理較大數(shù)據(jù)的潛力。實際上，對長度為131072的音樂片段(許嵩憂傷還是快樂截取，音樂采樣頻率為44100Hz)，使用冗余度的最優(yōu)對偶標架變換得到時頻圖和重構(gòu)誤差如

12、下：圖9 音樂片段最優(yōu)對偶標架變換時頻譜圖及重構(gòu)誤差圖考慮到該段音樂為鋼琴片段，圖9中的時頻譜是比較合理的，且程序運行時間只有26.65s(其中有限支撐最優(yōu)對偶標架的計算時間僅為0.004s)，再加上非常高的重構(gòu)精度，完全可以斷言最優(yōu)對偶標架變換具有處理實際數(shù)據(jù)的能力。最優(yōu)對偶標架變換另一個比較好的特點是，其自動選擇了一個比較好的窗寬(見式（22）)。首先對變化較慢的頻率為7.8125Hz的正弦信號，使用冗余度的最優(yōu)對偶標架變換得到時頻譜圖和重構(gòu)誤差圖如下：圖10 慢變信號最優(yōu)對偶標架變換時頻譜圖及重構(gòu)誤差圖其次，對變化劇烈的單位脈沖信號，仍使用冗余度的最優(yōu)對偶標架變換，得到其時頻譜圖和重構(gòu)誤

13、差圖如下：圖11 快變信號最優(yōu)對偶標架變換時頻譜圖及重構(gòu)誤差圖從圖10和圖11可以看出，最優(yōu)對偶標架自動選擇的窗寬可以同時較好地處理慢變信號和快變信號，這給使用帶來了較大的便利。最后需要說明的是冗余度和重構(gòu)精度的關系，事實上求解對偶標架的方程組式（17）的系數(shù)矩陣規(guī)模為，而。因此越大，對偶標架滿足的方程組就越欠定，所得的最優(yōu)對偶標架就接近緊標架，因而有較高的重構(gòu)精度(標架正反變換都線性有界(連續(xù))的充要條件是其為緊標架)。而當達到一定程度，使得解空間中已經(jīng)有很好的對偶標架時，重構(gòu)精度就應當保持穩(wěn)定了。上述討論可以從最優(yōu)對偶標架變換重構(gòu)精度隨冗余度的變化驗證，對長度為1024的100Hz正弦信號

14、，結(jié)果如下表所示：表1 最優(yōu)對偶標架變換重構(gòu)精度隨冗余度變化表13e-321e-1541e-1588e-16168e-16從表1中可以看出，當冗余度時，解空間中應該就已經(jīng)包含了很接近緊標架的對偶標架了。4. Grossman對偶標架變換4.1 原理分析不同于構(gòu)造方程組，Grossman給出了一個極其簡單的構(gòu)造性求解對偶標架的方法。為便于該方法的敘述，稱和為常數(shù)1的函數(shù)集為單位1的分割，定義平移算子為，定義調(diào)制算子為，則有：定理1：設是對單位1的分割，任意選擇，如果令，那么定義了一個Gabor分析標架，相應的綜合標架可為，這里定義為。Gross對偶標架變換實現(xiàn)的關鍵是單位1的分割，但實際上這可以

15、直接通過將一個窗函數(shù)離散向量求和歸一化實現(xiàn)。4.2 程序?qū)崿F(xiàn)附錄8給出了Grossman對偶標架變換的實現(xiàn)(其需要調(diào)用附錄6和附錄7中的程序)，對主頻為50Hz的Ricker信號，取采樣冗余度，使用Grossman對偶標架變換程序得到其時頻譜圖和重構(gòu)誤差如下：圖12 50Hz主頻Ricker信號Grossman對偶標架變換時頻圖和重構(gòu)誤差圖相對于最優(yōu)對偶標架變換，Grossman對偶標架變換節(jié)省了求解最優(yōu)對偶標架的時間，然而，當采樣冗余度較小時，Grossman對偶標架變換的重構(gòu)誤差精度較低，對長度為1024的100Hz正弦信號，結(jié)果如下表所示：表2 Grossman對偶標架變換重構(gòu)精度隨冗余

16、度變化表16e-123e-146e-283e-3165e-6322e-11645e-161285e-162565e-16對比表1和表2并結(jié)合計算量可得出結(jié)論，當采樣冗余度時，應選擇最優(yōu)對偶標架變換，當采樣冗余度時，應選擇Grossman對偶標架變換。5. 附錄附錄1：WFT的matlab程序%窗口傅里葉變換程序%s為長度為N的離散信號，w為窗函數(shù)指針%S為s的加窗變換,其分量S(i,j)(jlength(hTildeVec) error(采樣間隔過大，無法恢復信號！);end%計算所需常量L=length(hTildeVec);N=L/dN; %計算對偶標架頻率采樣點數(shù)M=L/dM; %計算對

17、偶標架時間采樣點數(shù)%初始化變量gammaTildeVec=zeros(1,L);hKMat=zeros(dN,M);muKBarVec=zeros(dN,1);muKBarVec(1)=1/N;%循環(huán)k計算最優(yōu)對偶標架函數(shù)for k=0:dM-1 %構(gòu)造dN行M列的系數(shù)矩陣 for q=0:dN-1 for l=0:M-1 hKMat(q+1,l+1)=hTildeVec(mod(k+l*dM+q*N,L)+1); end end %利用MP逆求解能量最小的解 gammaKBarVec=(pinv(hKMat)*muKBarVec); %注意這里包含共軛了 gammaTildeVec(k+(0

18、:M-1)*dM+1)=gammaKBarVec;endgammaTildeVec=real(gammaTildeVec);end附錄6：對偶標架變換的matlab程序%對偶標架變換程序%sTilde為長度為L的離散信號變量命名基于殷勤業(yè)講義%gammaTildeVec為最優(yōu)對偶標架函數(shù)離散向量%dM為對偶標架時間采樣間隔點數(shù)%dN為對偶標架頻率采樣間隔點數(shù)%m0為空間采樣起點，要求其值為0,dM-1間的整數(shù)%n0為波數(shù)采樣起點，要求其值為-L/2,-L/2+dN-1間的整數(shù)%cTildeMat為對偶標架系數(shù)矩陣function cTildeMat=dualFrameTransform(sTi

19、lde,gammaTildeVec,dM,dN,m0,n0)%判斷輸入的有效性if length(sTilde)=length(gammaTildeVec) error(要求對偶標架離散向量和信號維數(shù)相同！);endif mod(length(sTilde),dM)=0 error(dM必須為整數(shù)且能整除信號長度);endif mod(length(sTilde),dN)=0 error(dN必須為整數(shù)且能整除信號長度);endif rem(m0,1)=0 | m0dM-1 error(m0必須為0,dM-1-1間的整數(shù));endif rem(n0,1)=0 | n0-length(sTild

20、e)/2+dN-1 error(n0必須為-L/2,-L/2+dN-1間的整數(shù));end%計算所需常量L=length(sTilde);M=L/dM; %計算對偶標架空間采樣點數(shù)N=L/dN; %計算對偶標架波數(shù)采樣點數(shù)%計算對偶標架系數(shù)矩陣cTildeMatcTildeMat=zeros(M,N);for m=0:M-1 cache=fft(sTilde.*conj(gammaTildeVec(mod(0:(L-1)-m0-m*dM,L)+1); cTildeMat(m+1,:)=cache(mod(n0+(0:N-1)*dN,L)+1);endend附錄7：對偶標架逆變換的matlab程序

21、%逆對偶標架變換程序%cTildeMat為對偶標架系數(shù)矩陣，變量命名基于殷勤業(yè)講義%hTildeVec為綜合函數(shù)離散向量%dM為對偶標架時間采樣間隔點數(shù)%m0為空間采樣起點，要求其值為0,dM-1間的整數(shù)%n0為波數(shù)采樣起點，要求其值為-L/2,-L/2+dN-1間的整數(shù)%sTilde0為對偶標架重構(gòu)信號function sTilde0=inverseDualFrameTransform(cTildeMat,hTildeVec,dM,dN,m0,n0)%判斷輸入的有效性if rem(m0,1)=0 | m0dM-1 error(m0必須為0,dM-1-1間的整數(shù));endif rem(n0,1

22、)=0 | n0-length(hTildeVec)/2+dN-1 error(n0必須為-L/2,-L/2+dN-1間的整數(shù));end%計算所需常量M,N=size(cTildeMat);L=length(hTildeVec);unitImaginary=sqrt(-1);%重構(gòu)信號sTilde0=zeros(1,L);ifftCTildeMat=ifft(cTildeMat,2);for m=0:M-1 sTilde0=sTilde0+hTildeVec(mod(0:L-1)-m0-m*dM,L)+1).*ifftCTildeMat(m+1,mod(0:L-1,N)+1);endsTild

23、e0=N*exp(unitImaginary*n0*(0:L-1)*2*pi/L).*sTilde0;end附錄8：Grossman對偶標架變換的matlab程序%主程序clear;clc;%初始化變量dt=0.001;L=210;sTilde=sin(2*pi*100*(0:L-1)*dt).1); %構(gòu)造正弦信號%sTilde=sin(400*(0:L/4-1)*pi*dt),sin(800*(L/4:2*L/4-1)*pi*dt),sin(200*(2*L/4:3*L/4-1)*pi*dt),sin(600*(3*L/4:L-1)*pi*dt); %構(gòu)造四段跳頻信號%sTilde=(1-

24、2*(pi*(50)*(-L/2:L/2-1)*dt).2).*exp(-(pi*(50)*(-L/2:L/2-1)*dt).2); %ricker信號;%sTilde,fs,=wavread(截取憂傷還是快樂.wav);sTilde=transpose(sTilde(:,1);sTilde=sTilde(1:L);dt=1/fs; %實際音頻信號%計算Grossman對偶標架if fix(log2(length(sTilde)=log2(length(sTilde) error(信號的長度必須為2的整數(shù)次方);endL=length(sTilde);C1=256; %設置對偶標架變換的冗余度dM=2(ceil(log2(sqrt(L/C1);dN=L/(C1*dM); %模擬緊標架變換deltaX與