基于 MATLAB 進(jìn)行小波分析

先前，我曾發(fā)布過一篇名為《基于 Python 進(jìn)行小波分析》的推文，但其中的方法存在一定問題，目前該推文已刪除。在使用小波變換進(jìn)行分析時(shí)，需要將序列延長，以消除或減小邊界效應(yīng)．小波函數(shù)變換后再從中截取出實(shí)際長度列的小波系數(shù)進(jìn)行分析，而先前的方法并沒有延長序列相關(guān)的操作。故本文參考小文 Vince 的視頻，通過 MATLAB 進(jìn)行小波分析，從而解決這一問題。

先前，我曾發(fā)布過一篇名為《基于 Python 進(jìn)行小波分析》的推文，但其中的方法存在一定問題，目前該推文已刪除。在使用小波變換進(jìn)行分析時(shí)，需要將序列延長，以消除或減小邊界效應(yīng)．小波函數(shù)變換后再從中截取出實(shí)際長度列的小波系數(shù)進(jìn)行分析，而先前的方法并沒有延長序列相關(guān)的操作。故本文參考小文 Vince 的視頻（視頻鏈接見文末），通過 MATLAB R2016a 進(jìn)行小波分析，從而解決這一問題。

1 數(shù)據(jù)來源及下載方式

西北農(nóng)林科技大學(xué)的彭守璋研究員在國家青藏高原科學(xué)數(shù)據(jù)中心公開發(fā)布了氣溫、降水、干燥度等氣象數(shù)據(jù)，本文所使用的數(shù)據(jù)為基于其中的中國1km分辨率逐月降水量數(shù)據(jù)集 轉(zhuǎn)換得到的某地區(qū)歷年降水量數(shù)據(jù)，可通過 FTP 進(jìn)行下載，數(shù)據(jù)文獻(xiàn)引用格式詳見官網(wǎng)說明。

以 Filezilla 桌面版軟件為例，輸入官網(wǎng)中給定的主機(jī)號(hào)、用戶名、密碼、端口等信息后即可將數(shù)據(jù)從遠(yuǎn)程站點(diǎn)下載至本地。

2 計(jì)算小波系數(shù)

2.1 數(shù)據(jù)格式的轉(zhuǎn)化

為便于數(shù)據(jù)導(dǎo)入，可將原始數(shù)據(jù)在 Excel 表中保存為一列。通過單擊 MATLAB 操作界面上方【導(dǎo)入數(shù)據(jù)】按鈕，找到對(duì)應(yīng) Excel 表，以列矢量的形式將數(shù)據(jù)導(dǎo)入 MATLAB。

在工作區(qū)中右鍵單擊導(dǎo)入的數(shù)據(jù)，在彈出窗口中選擇【打開所選內(nèi)容】即可將數(shù)據(jù)以表格形式顯示，如下圖所示：

在工作區(qū)中右鍵單擊導(dǎo)入的數(shù)據(jù)，選擇【另存為...】將導(dǎo)入 MATLAB 中的數(shù)據(jù)另存為 mat 格式，以備后續(xù)操作使用。

2.2 數(shù)據(jù)序列延長

在命令行窗口中輸入【wavemenu】命令后回車，以打開小波分析工具箱，如下圖所示：

單擊工具箱右下角【Signal Extension】按鈕，進(jìn)入數(shù)據(jù)延伸窗口，通過窗口左上角的【File】→【Load Signal】導(dǎo)入上一步導(dǎo)出的 mat 格式表格。

導(dǎo)入數(shù)據(jù)后建議將右側(cè)的【Desired Lengt】修改為【Next Power of 2】的兩倍，【Extension Mode】建議設(shè)置為 【Symmetric (Whole-Point)】，修改完成后單擊【Extend】按鈕延伸數(shù)據(jù)，延伸結(jié)果如下圖所示：

注：示例數(shù)據(jù)為 122 年的降水量數(shù)據(jù)，故 Length 為 122，122 之后下一個(gè) 2 的指數(shù)冪為 128，長度修改為 256 后，其前后各延伸出 (256-122)/2=67 個(gè)數(shù)據(jù)；如果 Length 為單數(shù)，以 123 年為例，則左側(cè)延伸出 67 個(gè)數(shù)據(jù)，右側(cè)延伸出 68 個(gè)數(shù)據(jù)。

通過窗口左上角的【File】→【Save Transformed Signal】將延伸后的數(shù)據(jù)以 mat 格式導(dǎo)出。

2.3 計(jì)算小波系數(shù)

單擊小波分析工具箱左上角的【Complex Continuous Wavelet 1-D】按鈕，進(jìn)入小波分析窗口計(jì)算小波系數(shù)，通過窗口左上角的【File】→【Load Signal】導(dǎo)入上一步導(dǎo)出的 mat 格式延長后數(shù)據(jù)表格。

【W(wǎng)avelet】函數(shù)建議使用【cmor】，修改完成后單擊 Analyze 按鈕運(yùn)行，運(yùn)行結(jié)果如下圖所示：

通過窗口左上角的【File】→【Save Coefficients】將小波變換結(jié)果數(shù)據(jù)以 mat 格式導(dǎo)出后，關(guān)閉小波分析工具箱。

2.4 計(jì)算小波系數(shù)的實(shí)部、模、模方及方差

將小波系數(shù) mat 文件導(dǎo)入 MATLAB，導(dǎo)入后工作區(qū)中的 coefs 數(shù)據(jù)即為小波系數(shù)。此時(shí)的小波系數(shù)包含延長后的數(shù)據(jù)，故需將其數(shù)據(jù)左右兩側(cè)各 67 列數(shù)據(jù)刪除。右鍵單擊 coefs 數(shù)據(jù)，在彈出窗口中選擇【打開所選內(nèi)容】以表格形式顯示數(shù)據(jù)。選中需要?jiǎng)h除的數(shù)據(jù)列后單擊操作界面上方【變量】選項(xiàng)卡下【刪除】下拉菜單中的【刪除列】即可刪除選中列，如下圖所示：

在命令行窗口中通過 real()、abs () 等函數(shù)計(jì)算小波系數(shù)的實(shí)部、模、模方及方差：

coefs_real = real(coefs);
coefs_abs = abs(coefs);
coefs_abs2 = coefs_abs.^2;
coefs_variance = sum(coefs_abs2, 2);

以表格形式顯示工作區(qū)中生成的各數(shù)據(jù)，全選表格中的數(shù)據(jù)，將其復(fù)制到 Excel 中。

3 繪制小波實(shí)部等值線圖及方差圖

為 Excel 中的小波系數(shù)實(shí)部添加橫縱坐標(biāo)，本例橫坐標(biāo)為年索引，縱坐標(biāo)為時(shí)間尺度，如下圖所示：

將實(shí)部表格所有數(shù)據(jù)復(fù)制粘貼到 Origin 中進(jìn)行圖表的制作。在 Origin 中全選數(shù)據(jù)，繪制等高線圖，等高線圖參數(shù)設(shè)置如下：

在Excel中將方差列前插入時(shí)間尺度列，將兩列表格復(fù)制粘貼到 Origin，選中方差列制作折線圖。

對(duì)圖表風(fēng)格進(jìn)行修飾后，得到小波系數(shù)實(shí)部等值線圖及方差曲線圖示例如下圖所示：

參考資料：

• Bilibili up 主：小文Vince，視頻鏈接
• 左斌斌, 徐宗學(xué), 任梅芳, 等. 北京市通州區(qū)1966—2016年降水特性研究[J]. 北京師范大學(xué)學(xué)報(bào)(自然科學(xué)版), 2019, 55(5): 556-563.
• 王勁峰, 廖一蘭, 劉鑫. 空間數(shù)據(jù)分析教程(第二版)[M]. 北京: 科學(xué)出版社, 2019.