參考：A Practical Guide to Wavelet Analysis
本文主要學習上述鏈接文章，并提出自己的一點思考和疑問。

目錄

• • 小波分析（Wavelet Analysis）基本原理
  • • 1.Introduce
    • 2.Wavelets
    • 3.Algorithms
    • 4.Monte Carlo
    • 5.References
  • 相關MATLAB代碼

小波分析（Wavelet Analysis）基本原理

1.Introduce

地球物理學中許多時間序列在統計上表現出非平穩性。雖然序列可能包含主導周期信號，但這些信號的振幅和頻率可以在很長一段時間內變化。
一個例子是赤道太平洋的海表溫度(圖1)。變率的主導“模式”是El Ni?o-Southern振蕩(ENSO)，在2-7年的時間尺度上顯示為高頻“尖峰”。疊加在這個信號上的是更長的年代際波動。年代際波動具有調節El Ni?o事件發生幅度和頻率的作用。

圖1 東太平洋(5°S-5°N, 90°W-150°W) NINO3區域的海表平均溫度較高。藍色曲線為低通濾波(>12個月)的海表溫度。黃色背景曲線是15年的方差，在15年周期的中間點。曲線被倒過來復制，以顯示方差的“包絡”

理想情況下，人們希望將較短周期的振蕩與較長周期的振蕩分開。
注:此處使用“振蕩”一詞來表示時間序列中任何重復的波動，無論這種波動是否有規律地重復。在我們的例子中，El Ni?o絕對是一個非周期(不規則)信號。

分析時間序列的非平穩性最簡單的方法是計算統計數據，如不同時間段的均值和方差，看看它們是否有顯著差異。

在圖1中，我們還繪制了運行的15年方差，作為信號相對于時間的固有總功率的測量。我們可以看到，在1880-1920年和1950年之后，ENSO有更多的變化，在1920-1950年期間有一個相對平靜的時期。
雖然運行方差告訴我們信號在特定時間的總體強度，但它存在兩個主要缺陷:
【時間定位】曲線的形狀高度依賴于所用窗口的長度。上面選擇的15年是平滑度過高(比如使用30年窗期)和過低(5年窗期)之間的折衷。理想的方法是根據感興趣的尺度允許不同的窗口大小。
【頻率定位】運行方差不包含周期信號的頻率信息，只包含它的振幅(并且只有當窗口被廣泛選擇時)。

一種可能是做一個窗口(或運行)傅里葉變換(WFT)，使用特定的窗口大小并在時間上滑動它，每次只使用窗口內的數據計算FFT。這將解決第二個問題(頻率本地化)，但仍然取決于所使用的窗口大小。
WFT的主要問題是對不同頻率的不一致處理:在低頻時，窗口內的振蕩太少，導致頻率局部化丟失;而在高頻時，振蕩太多，導致時間局部化丟失。最后，WFT依賴于信號可以分解為正弦分量的假設。

小波分析試圖通過同時將時間序列分解為時間/頻率空間來解決這些問題。我們既可以得到級數中任何“周期”信號的振幅信息，也可以得到振幅隨時間的變化情況。

2.Wavelets

在前一節中，我們看到了如何使用固定寬度窗口計算運行方差來衡量時間序列的平穩性。雖然我們指出了使用固定寬度窗口的缺點，但可以使用不同的窗口寬度重復分析。通過平滑地改變窗口寬度，我們就可以建立方差隨時間和窗口寬度變化的圖像。這種技術的明顯問題是窗口函數的簡單“車廂”形狀，這引入了邊緣效應，如振鈴。正如前面提到的，使用這樣一個黑匣子汽車，我們仍然不知道盒子里發生了什么，而只是恢復了平均能量。
在圖2a中，我們看到了一個波“包”的例子，它具有有限的持續時間和特定的頻率。你可以想象用這樣一個形狀作為我們方差分析的窗口函數。這種“小波”的優點是包含了一定周期的波，并且在范圍上是有限的。事實上，圖2a所示的小波(稱為Morlet小波)只不過是一個正弦波(圖2b中的綠色曲線)乘以一個高斯包絡線(紅色曲線)。
圖2 (a)沿x軸隨時間變化的任意寬度和振幅的Morlet小波。(b)構造Morlet小波(藍色虛線)作為正弦曲線(綠色)調制的高斯曲線(紅色)。

圖3 (a) El Ni?o海表溫度時間序列。(b)小波功率譜，使用Morlet小波。
x軸是小波在時間上的位置。y軸為以年為單位的小波周期。黑色輪廓是10%顯著性區域，使用紅色噪聲背景光譜。紅色區域表明，1880-1920年和1965-現在的高厄爾尼諾Ni?o活動發生，而1920-1960年則相對平靜。

3.Algorithms

4.Monte Carlo

5.References

相關MATLAB代碼

數據：Nino3 區海表溫度（SST）序列，1870-1996共126年的季尺度數據

求滑動窗口為十年時的方差包絡線,如下圖所示:
可知,SST為非平穩時間序列

相關代碼如下:

[E_SSTA , D_SSTA] = GetED( sst ,10*4 ); figure(1) hold on H1=bar(5*4+1:5*4+length(D_SSTA),D_SSTA); H2=bar(5*4+1:5*4+length(D_SSTA),-D_SSTA); set(H1,'BarWidth',1 ,'facecolor',[1 1 0]); set(H2,'BarWidth',1 ,'facecolor',[1 1 0]); plot(5*4+1:5*4+length(D_SSTA),D_SSTA,'k-','linewidth',1); plot(5*4+1:5*4+length(D_SSTA),-D_SSTA,'k-','linewidth',1); a = bar(sst,1); xlabel('Time (year)') ylabel('NINO3 SST (degC)') axis([1 length(sst) -3 4]) set(a,'BarWidth',1 ,'facecolor',[0 0 0.8]); x = 1.5:5*4:126*4; y = cell(1,length(x)); y(1:4:end)={'1870','1890','1910','1930','1950','1970','1990'}; set(gca, 'XTick',x,'XTickLabel',y) ;function [ E , D ]=GetED(X,n) % GetED函數用于獲取滑動窗口為n時的均值和方差 % X為輸入數據序列 % n為滑動窗口長度 N = length(X); if N<nprintf("錯誤：滑動窗口長度大于數據序列長度。") endE = zeros(N-n+1,1); D = zeros(N-n+1,1);E(1) = mean( X(1:n) ); D(1) = var( X(1:n) ); for i=1:N-nE( i+1) = E(i)+ (X(i+n)-X(i) )/n ;D(i+1) = D(i) +(E( i+1) -E(i) )/n*( (X(i+n)+X(i))*(n-1)-2*n*E(i) + 2*X(i)); endend

對時間序列進行連續小波變換,得到如下結果:

相應代碼如下:(使用了工具箱)

wt(sst); xlabel('Time (year)') ylabel('Period(year)') x = 1.5:5*4:126*4; y = cell(1,length(x)); y(1:4:end)={'1870','1890','1910','1930','1950','1970','1990'}; set(gca, 'XTick',x,'XTickLabel',y) ;

相應代碼如下:(同上,使用了工具箱)

figure(1) subplot('position',[0.1 0.35 0.85 0.6]) wt(sst); ylabel('Period(year)') x = 1.5:5*4:126*4; y = cell(1,length(x)); set(gca, 'XTick',x,'XTickLabel',y) ; subplot('position',[0.1 0.1 0.745 0.25]) plot(1:length(sst),sst,'k-'); xlabel('Time (year)') ylabel('NINO3 SST (degC)'); axis([0 length(sst) -2.8 3.5] ); x = 1.5:5*4:126*4; y = cell(1,length(x)); y(1:4:end)={'1870','1890','1910','1930','1950','1970','1990'}; set(gca, 'XTick',x,'XTickLabel',y) ;

根據該網站相應MATLAB代碼得到以下結果：

總結

以上是生活随笔為你收集整理的小波分析及相应MATLAB实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。