本專欄中講了很多時頻域分析的知識，不過似乎還沒有講過時頻域分析是怎樣引出的。

所以本篇將回歸本源，講一講從傅里葉變換→短時傅里葉變換→小波分析的過程。

為了讓大家更直觀得理解算法原理和推導(dǎo)過程，這篇文章將主要使用圖片案例。

一、頻譜分析？——還不夠

頻譜分析可以告訴我們一個信號包含哪些頻率的信息以及這些頻率的強度。

通過頻譜分析我們可以將信號從其原始的時間域（即隨時間變化的形式）轉(zhuǎn)換到頻域（即按頻率分布的形式）。

如果對此你不懂的話，那么我推薦你讀一下這篇文章：Heinrich：如果看了這篇文章你還不懂傅里葉變換，那就過來掐死我吧。

頻譜分析是如此常用的工具，如果你在做信號處理的研究領(lǐng)域，幾乎不可能沒聽過它。

然而，傳統(tǒng)的傅里葉變換有一個限制：它假設(shè)信號是平穩(wěn)的，即其頻率成分不會隨時間改變。

對于那些隨時間變化其頻率特性的信號（比如音樂或語音），這種假設(shè)并不成立，因此傳統(tǒng)的頻譜分析無法準確地描繪這些信號的頻率隨時間的變化情況。

舉個可能不太恰當?shù)睦樱斈愠皩m廷玉液酒”和“玉液宮廷酒”時，聲音信號可能在頻譜上差別不大。

我們可以用一張圖演示這個現(xiàn)象：

“宮廷玉液酒”和“玉液宮廷酒”——兩個不同的信號，他們的頻譜是類似的

可見當信號的頻率成分隨時間顯著變化時，即所謂的非平穩(wěn)信號，傳統(tǒng)頻譜分析就不太合適了。舉幾個例子：

• 研究信號的局部特性：如爆炸聲、機器的故障噪聲等，這些信號的特性在短時間內(nèi)變化很大。
• 語音處理：在自動語音識別和語音合成中，語音的特征如音高和音量隨著時間而變化。
• 雷達和無線電信號分析：雷達信號的特性依賴于時間和觀測的角度，需分析這些信號隨時間的變化。
• 生物醫(yī)學(xué)信號分析：比如心電圖（ECG）和腦電圖（EEG）這樣的生理信號，它們的頻率特性隨時間改變。

二、直觀的解決思路——給頻譜分析加窗

傳統(tǒng)的頻譜分析是全面的，它給出了整個信號的平均頻率內(nèi)容，卻無法告訴我們這些頻率成分是怎樣隨時間變化的。

這就好比是我們閉上眼睛聽整首曲子，然后試圖回憶每個樂器何時加入。這個任務(wù)對于我們的大腦來說是非常困難的，對于數(shù)學(xué)工具來說也同樣困難。

這里就需要一種方法來解決這個問題，而“加窗”提供了一個直觀的解決思路。

加窗，簡單來說，就是在分析信號時不是一次性看整個信號，而是只看信號的一個小段，這個小段就是一個“窗口”。

通過在信號上滑動這個窗口，我們可以逐漸查看整個信號，就像是通過一系列快照來觀察它的變化。

這個過程可以用一個日常的比喻來說明：想象你正在看一場煙花表演，你決定用相機拍攝整個表演。如果使用長時間曝光，你最終得到的照片將是所有煙花的疊加，這相當于傳統(tǒng)的頻譜分析；然而，如果你使用一連串短時間曝光的快照，每張照片將捕捉到表演的一個瞬間，這就類似于“加窗”的頻譜分析。

在實際應(yīng)用中，我們選擇一段時間內(nèi)的信號，對它進行傅里葉變換，這樣得到的頻譜只代表了那段時間內(nèi)的信號。然后，我們移動窗口到下一個時間段，再做一次變換，如此反復(fù)。最后，我們將所有的頻譜組合起來，就能夠得到一個動態(tài)的頻譜圖，這個圖不僅展示了信號的頻率成分，還展示了它們是如何隨時間變化的——這個過程就是短時傅里葉變換（Short-Time Fourier Transform， STFT）。

其計算過程可以用一張動圖生動地展示出來：

短時傅里葉變換STFT計算過程

STFT圖就是頻譜圖在時間軸上的連續(xù)展開，這個圖有三個軸：

• 時間軸：表示信號分析的持續(xù)時間。
• 頻率軸：它表示了信號中存在的各種頻率成分。
• 幅度軸：幅度軸垂直于時間-頻率平面，通常通過顏色的深淺來表示。強度表示在特定時間和頻率下信號的能量大小。

STFT圖也有二維形式，這種形式中是使用顏色來代表幅度大小的，其生成過程如下：

二維STFT圖繪制過程

返回到上邊“宮廷玉液酒”和“玉液宮廷酒”的例子，我們畫一下其STFT圖：

在STFT中可以看到兩個信號的區(qū)別了

對比第一個信號和第二個信的STFT圖，可以看到在不同時間段內(nèi)頻率成分的變化。

通過比較這些圖，可以直觀地看出STFT與傳統(tǒng)傅里葉變換的區(qū)別。STFT揭示了信號的局部時頻特性。這就是為什么STFT特別適用于顯示信號中頻率成分隨時間的變化。

三、成也窗口，敗也窗口——從短時傅里葉變換到小波分析

在短時傅里葉變換（STFT）中，窗口函數(shù)及其大小選擇是分析的關(guān)鍵。窗口函數(shù)決定了在任何給定時間點，信號的哪一部分被用于分析。窗口大小的選擇直接影響了分析結(jié)果的時間分辨率和頻率分辨率，這是進行有效STFT分析的最重要的權(quán)衡。

窗口大小的時間分辨率影響：時間分辨率與窗口的寬度密切相關(guān)。一個窄窗口提供較高的時間分辨率，因為它捕捉了信號在很短時間內(nèi)的變化。這對于分析包含快速變化的瞬態(tài)事件，如敲擊聲或爆炸聲，是非常有用的。然而，較小的窗口將限制頻率分辨率，因為頻率分析需要足夠的周期來準確估計。

窗口大小的頻率分辨率影響：頻率分辨率與窗口的寬度呈反比。一個寬窗口覆蓋了信號的較長時間段，提供了較高的頻率分辨率。這是因為更多的周期可以在窗口內(nèi)被分析，從而更準確地確定低頻成分。但是，這會犧牲時間分辨率，因為窗口中的信號被假定在這段時間內(nèi)是平穩(wěn)的。

只靠文字可能不夠直觀，還是來看圖吧：

• 小窗口長度（例如64樣本）： 小窗口提供較高的時間分辨率，這意味著可以更精確地看到頻率變化的時間位置（即0.5s-0.6s）。例如，短暫的頻率跳變將在圖中更清晰地顯示為突然的顏色變化。然而，頻率分辨率較低，這意味著看不清楚具體是哪個頻率發(fā)生了變化，頻率成分可能會顯得模糊。
• 中等窗口長度（例如128樣本）： 中等大小的窗口提供了時間分辨率和頻率分辨率之間的折衷。會看到頻率跳變不如小窗口那樣清晰，但頻率的分辨率會更好一些，使得不同的頻率成分更加明顯。
• 大窗口長度（例如256樣本）： 大窗口提供較高的頻率分辨率，您將能夠看到信號中不同頻率成分的細節(jié)，就像200Hz的頻率可以被很好得識別出來。然而，時間分辨率較低，這意味著短暫事件的精確時間定位可能不太清晰，圖中跳變的時間段都延展到0.4-0.6s了，而實際上它是0.5-0.6s。

總結(jié)一句話，上述矛盾體現(xiàn)了所謂的測不準原理，即時空域和頻域無法同時準確。

在實際的信號處理工作中，選擇窗口大小常常是基于對信號特性的了解以及對分析精度需求的權(quán)衡。如果事件的時間特性是分析的關(guān)鍵，那么小窗口可能更合適；如果頻率分辨率是關(guān)鍵，那么大窗口可能更優(yōu)。

那有沒有一種可能，窗口大小是可調(diào)的呢？

是的——小波分析就可以實現(xiàn)！

四、小波分析到底是怎么計算的

4.1 又雙叒叕說卷積

“卷積”的概念在本專欄里是明星選手了，還不熟悉卷積的同學(xué)可以看這里：Mr.看海：這篇文章能讓你明白卷積。

我們回憶一下卷積的重要特性：頻域相乘等于時域的卷積。

那么頻域分析我們可以從一種獨特的視角來看：想象一下，我們從頻率為0的正弦波開始，逐漸增加這個頻率，同時記錄每個頻率的正弦波與原始信號卷積的結(jié)果。這樣，我們可以繪制出一張圖，其橫軸表示頻率，縱軸表示卷積結(jié)果的強度。請問得到的結(jié)果是什么？

是的，就（近似）是頻譜！

為什么？因為1Hz的頻譜為一個在1Hz處的脈沖，這個脈沖乘以原始信號的頻譜得到的就是原始信號1Hz處的頻譜幅值，轉(zhuǎn)換到時域上就是原始信號與這個1Hz的正弦信號的卷積。

動圖又來了，看過后相信你就能明白：

頻譜計算的一種理解視角——頻譜可以被視為一系列不同頻率的正弦波與原始信號進行卷積的結(jié)果。

然而，這種方法只能提供頻率上的信息，畢竟中間圖中的正弦無論如何平移，都不會影響到頻譜的結(jié)果。

聰明的你應(yīng)該想到解決方案了。

4.2 當小波遇上卷積

是的，我們可以將無限延展周期性的正弦波換成時變的信號，比如這樣的：

Meyer小波的圖像

這樣一個信號，就同時可以做拉伸（頻率上變換）和平移（延遲時間上變換）了。

依然采用上述卷積計算的方式，來求小波變換的結(jié)果，此時畫出來的圖就不是頻譜這樣的二維圖，而是一個三維圖了，就像下邊這樣：

小波變換過程（示意圖，計算可能不太嚴謹）

我們重點看中間的小波的圖像，它在時間軸上進行平移，每平移一次進行一次卷積運算，計算結(jié)果在第三張圖上按照對應(yīng)的時間（由位置決定）和頻率（由尺度決定）畫出計算結(jié)果。仔細觀察小波的變換，可以發(fā)現(xiàn)其特點：高頻部分具有較高的時間分辨率和較低的頻率分辨率，而低頻部分具有較高的頻率分辨率和較低的時間分辨率，這就恰好解決了STFT的痛點。（PS：小波是復(fù)信號，上圖中只畫了實部投影）

再回到反映STFT中時空域和頻域測不準原理這張圖，如果我們用小波分析處理同樣的信號，得到的結(jié)果是這樣的：

增加了最后一張圖——CWT變換結(jié)果

在最下邊這張圖中可以清晰地顯示信號中的所有頻率分量及其隨時間的變化。三個突變信號的邊緣清晰銳利、三個主要頻率的分辨也較為清晰可辨，在時間分辨率和頻率分辨率的取舍中找到了完美的平衡點。

最后以一個優(yōu)美（誤）的比喻來結(jié)束理論講解部分吧：

頻譜分析是長曝光的藝術(shù)，揭示了振動的連續(xù)舞蹈；短時傅里葉變換像是定焦的高速連拍，捕捉時間的每一跳動。而小波分析賦予了相機變焦的魔力，讓每幀都能自適應(yīng)地聚焦于故事的精髓，無論是宏偉的篇章還是微妙的細節(jié)。

五、短時傅里葉變換（STFT）和小波分析（CWT）的MATLAB代碼實現(xiàn)

筆者對STFT和CWT算法畫圖程序分別進行了封裝，同學(xué)們在調(diào)用的時候可以方便輕松很多。

5.1 STFT畫圖（二維圖+三維圖）

短時傅里葉變換畫圖常用設(shè)置參數(shù)就四個：窗口類型、窗口大小、窗口重疊大小、每個窗口上執(zhí)行FFT的點數(shù)。選擇合適的參數(shù)對于獲取準確和有用的頻譜分析至關(guān)重要。以下是這四個參數(shù)的詳細說明及其對STFT結(jié)果的影響：

1. 窗口類型：窗口類型決定了信號在時頻分析中的平滑程度。常用的窗口類型包括漢明窗、漢寧窗和矩形窗。不同的窗口類型會影響頻譜的泄露和分辨率。例如，矩形窗提供最高的頻率分辨率，但其頻譜泄漏也最大；而漢明窗和漢寧窗則提供更好的頻譜平滑性，適用于波形變化不是非?？斓男盘?。
2. 窗口大小：窗口大小決定了頻率分辨率和時間分辨率之間的權(quán)衡。較大的窗口提供更高的頻率分辨率但較低的時間分辨率，而較小的窗口則相反。
3. 窗口重疊大小：窗口重疊用于提高時間分辨率并減少數(shù)據(jù)丟失。較大的重疊程度可以提供更平滑的頻譜表示，但計算成本也更高。
4. 每個窗口上執(zhí)行FFT的點數(shù)：這決定了頻譜的頻率分辨率。通常，F(xiàn)FT點數(shù)至少應(yīng)等于窗口大小，它通常取2的整數(shù)次冪，比如128、256這種。

好了，知道要設(shè)置的參數(shù)，代碼實現(xiàn)就很簡單了，就像這樣：

% 設(shè)置STFT參數(shù)
windowType = @hamming;  % 設(shè)置窗口類型為 Hamming 窗口
windowSize = 100;       % 設(shè)置窗口大小為 100 樣本
overlap = 98;           % 設(shè)置窗口重疊為 98 樣本
FFTLength = 128;        % 設(shè)置在每個窗口上執(zhí)行FFT的點數(shù)為 128

% 調(diào)用 pSTFT2D 函數(shù)
% 使用上述設(shè)置對數(shù)據(jù)進行短時傅里葉變換，并繪制二維圖形
pSTFT2D(data, Fs, windowType, windowSize, overlap, FFTLength);

此時就可以畫出這樣的結(jié)果：

STFT的二維結(jié)果圖

如果想畫三維結(jié)果也可以，只需要寫：

% 調(diào)用 pSTFT3D 函數(shù)
% 使用相同的設(shè)置對數(shù)據(jù)進行短時傅里葉變換，并繪制三維圖形
pSTFT3D(data, Fs, windowType, windowSize, overlap, FFTLength);

此時可以畫出：

STFT的三維結(jié)果圖

是不是既簡單又酷炫。上邊代碼中的pSTFT2D和pSTFT3D是筆者封裝的函數(shù)，文末可以獲取哈。

5.2 CWT畫圖（二維圖+三維圖）

連續(xù)小波變換的圖需要設(shè)置的參數(shù)就比較少了，這就很適合選擇困難癥的同學(xué)們。

通常情況下，只需要選擇小波類型就可以，甚至類型也不需要費心選，使用默認的Morse就可以。

不過同學(xué)們也許會在論文中看到CWT的圖是這樣的：

連續(xù)小波變換圖

在這個圖上，有一個錐形或者半圓形的邊界。這種邊界通常表示“錐形區(qū)域”（Cone of Influence, COI），用于指出小波變換結(jié)果中的邊緣效應(yīng)。

在小波分析中，信號的開始和結(jié)束部分可能會因為缺乏足夠的數(shù)據(jù)而產(chǎn)生不可靠的變換結(jié)果。由于小波變換涉及到信號的縮放和平移，信號的邊緣部分在變換過程中可能會“跑出”數(shù)據(jù)的實際范圍，導(dǎo)致這些區(qū)域的變換結(jié)果被信號之外的零填充或其他形式的邊界處理影響。錐形區(qū)域邊界內(nèi)的變換結(jié)果被認為是可靠的，而邊界外的變換結(jié)果可能會受到邊緣效應(yīng)的影響，因此在解釋這部分結(jié)果時需要格外小心。

畫這樣一張圖需要四行代碼：

% 情況1：繪制錐形區(qū)域
plotCone = true;     % 設(shè)置繪制錐形區(qū)域為 true
freqScale = 'log'; % 當繪制錐形區(qū)域時，freqScale 參數(shù)不影響輸出，都會為對數(shù)坐標
waveletName = 'morse'; % 設(shè)置小波名稱為 Morse
% 調(diào)用 pCWT 函數(shù)，傳入上述參數(shù)并獲取結(jié)果
[timeAxis, freqAxis, cwtData] = pCWT(data, waveletName, Fs, plotCone, freqScale);

有些同學(xué)表示，不想要錐形區(qū)域，就要單純的CWT結(jié)果，這個我也給出了選擇余地——把plotCone設(shè)置為false即可：

% 情況2：不繪制錐形區(qū)域，使用對數(shù)頻率軸
plotCone = false;     % 設(shè)置繪制錐形區(qū)域為 false
freqScale = 'log';   % 設(shè)置頻率軸類型為對數(shù)
waveletName = 'morse'; % 設(shè)置小波名稱為 Morse
% 調(diào)用 pCWT 函數(shù)，傳入上述參數(shù)并獲取結(jié)果
[timeAxis, freqAxis, cwtData] = pCWT(data, waveletName, Fs, plotCone, freqScale);

去掉了錐形區(qū)域

有細心的同學(xué)發(fā)現(xiàn)了，上邊兩張圖的頻率軸是取對數(shù)的結(jié)果，但是我就是想要線性的——也可以！把freqScale設(shè)置為'linear'：

% 情況3：不繪制錐形區(qū)域，使用線性頻率軸
plotCone = false;    % 設(shè)置不繪制錐形區(qū)域
freqScale = 'linear';% 設(shè)置頻率軸類型為線性
waveletName = 'morse'; % 設(shè)置小波名稱為 Morse
% 調(diào)用 pCWT 函數(shù)，傳入上述參數(shù)并獲取結(jié)果
[timeAxis, freqAxis, cwtData] = pCWT(data, waveletName, Fs, plotCone, freqScale);

線性頻率軸

有同學(xué)又表示，我要像STFT一樣要三維圖！好吧，我也準備了（而且可以選擇頻率軸線性坐標或者對數(shù)坐標）：

% 情況1：使用線性頻率軸繪制三維圖
freqScale = 'linear';
waveletName = 'morse'; % Morse 小波
[timeAxis, freqAxis, cwtData] = pCWT3D(data, waveletName, Fs, freqScale);


% 情況2：使用對數(shù)頻率軸繪制三維圖
freqScale = 'log';
waveletName = 'morse'; % Morse 小波
[timeAxis, freqAxis, cwtData] = pCWT3D(data, waveletName, Fs, freqScale);

頻率軸為線性坐標的三維CWT圖

頻率軸為對數(shù)坐標的三維CWT圖

至此，大家的需求應(yīng)該都可以滿足了吧！

獲取代碼

上邊的測試代碼和封裝函數(shù)，包括工具箱都可以在下邊的鏈接獲?。?/strong>