參考：A Practical Guide to Wavelet Analysis
本文主要學(xué)習(xí)上述鏈接文章，并提出自己的一點(diǎn)思考和疑問。

目錄

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

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

1.Introduce

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

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

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

分析時(shí)間序列的非平穩(wěn)性最簡單的方法是計(jì)算統(tǒng)計(jì)數(shù)據(jù)，如不同時(shí)間段的均值和方差，看看它們是否有顯著差異。

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

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

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

2.Wavelets

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

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

3.Algorithms

4.Monte Carlo

5.References

相關(guān)MATLAB代碼

數(shù)據(jù)：Nino3 區(qū)海表溫度（SST）序列，1870-1996共126年的季尺度數(shù)據(jù)

求滑動(dòng)窗口為十年時(shí)的方差包絡(luò)線,如下圖所示:
可知,SST為非平穩(wěn)時(shí)間序列

相關(guān)代碼如下:

[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函數(shù)用于獲取滑動(dòng)窗口為n時(shí)的均值和方差 % X為輸入數(shù)據(jù)序列 % n為滑動(dòng)窗口長度 N = length(X); if N<nprintf("錯(cuò)誤：滑動(dòng)窗口長度大于數(shù)據(jù)序列長度。") 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

對時(shí)間序列進(jìn)行連續(xù)小波變換,得到如下結(jié)果:

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

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) ;

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

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) ;

根據(jù)該網(wǎng)站相應(yīng)MATLAB代碼得到以下結(jié)果：

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。