振動分析在工程領域很常見，汽車中的 NVH，風機的模態分析，塔架、葉片等機械件的振動頻率分析，發電機軸承振動，飛機發動機階次分析等等。經常提到的知識點像信號處理，FFT，階次分析，頻譜估計，包絡譜分析，模態分析，頻率響應函數估計等等。

傅里葉變換

要談信號處理就離不開傅里葉變換。我們先從最常用的傅里葉變換開始。傅里葉變換的本質是從空間（希爾伯特空間）中找一組正交基向量 或者正交基函數，正交意思是內積為 0，其中希爾伯特空間內積定義為：對于向量 f， g， 有

容易證明

將時域信號投影到這組正交基上。傅里葉變換也就可以看作 f（t） 在 k 對應的基 上的投影。對于離散傅里葉變換（DFT）

為了方便理解，我們先通過一個例子，利用 fft 函數（快速傅里葉變換，實現DFT的一種快速算法）將一段矩形信號（圖中深藍色信號）進行變換。通常我們得到幅值頻率圖像，可以看到對應頻率的幅值。

為了更好的理解，我們把離散頻率對應的正交基向量（圖淺綠色）也繪制出來，對應圖中淺綠色。方塊綠則為幅值頻率值，此處我們首先只選了前 20 個頻率對應的基向量進行信號重構，得到圖中紅色信號。接下來，我們增加用于信號重建的基的數量，分別我們取前 60，前 300，發現信號越來越接近原始信號。

圖表2

值得一提的是很多類似的程序會通過for循環來遍歷各個基向量來實現上面的程序，但 MATLAB 中的矩陣思想更應該充分利用，既可以簡化代碼，又可以提升效率，例如程序中 baseSeries=cos（2*pi*Dfreq‘*t+Phi（1:length（Dfreq））’）計算得到的 baseSeries 變量是一個 n*m 矩陣，包含了所有的基向量n 是基向量的個數（k 的個數），m 是時域信號長度（t 的長度）。簡單理解了傅里葉變換，接下來我們就可以使用 fft 進行一個初步的轉動數據分析，我們先生成一段時長 2s，采樣頻率 600Hz 的數據，在這 2s 中，系統一直以 20Hz 的轉速運行。

通過 fft 我們可以得到信號在對應轉速的 20Hz 頻率上幅值最大。這部分有很多重要的概念，包括時域/頻域分辨率，頻譜泄露，功率譜/能量譜等等。頻率分辨率簡言之就是區分兩個不同頻率的能力。假設信號的采樣頻率為 Fs， 在時域上，tn= n*Fs 反映的是時間軸，T=N*Fs 反映的是信號總時長，在頻域上，我們可以得到 N 個頻率（基向量），可以理解為 Fs 被 N 個值平分，因此頻率的分辨率為

我們從上面推導可以知道頻率分辨率由采樣時長決定，如果采樣時長太短，DFT 沒有能力區分兩個接近的頻率值。

這里也給大家一個小練習，例如可以合成一段信號，由 90Hz 和 100Hz 兩個頻率組成，我們只選取 0.05s 長的數據：對應 DFT 的頻率分辨率為 20Hz 》 （100Hz – 90Hz），然后進行 fft 變換，查看是否在頻域上能夠區分開這兩個信號成分。

接下來我們通過一個由兩個頻率組成的正弦信號來介紹頻譜泄露。詳見 https://www.mathworks.com/help/releases/R2020b/signal/ug/amplitude-estimation-and-zero-padding.html這兩種正弦波的頻率分別為 f1 = 100Hz 和 f2 = 202.5 Hz。采樣率為 1000Hz，信號長度為 1000 個采樣點（T=1s）。我們期望的結果是頻譜幅值只在 f1 = 100Hz 和 f2 = 202.5Hz 處的結果不為零，其他位置的頻譜幅值為零。但真實結果是這樣的：

圖表4

也就是在 202.5Hz 的信號幅值被周圍的頻率值平分了，這種現象就是頻譜泄露。從頻譜分辨率的角度解釋通過上面的結論：我們得到DFT的頻率分辨率為 1/T = 1Hz，對于 f1 = 100Hz 正好落在分辨率的分倉邊界上，而 f2 = 202.5Hz 不能被分辨率 1Hz 取整。

從時域角度解釋：我們取了窗口時間 1s 長度的數據，對應 f1 的 100 個完整周期以及 f2 的 202.5 個周期，DFT 默認信號是采樣時間窗內周期出現的無限長信號，因此當窗重復時對于f1來說正好是整數拼接，與原信號完全一致，而 f2 在上一個 1s 采樣時間窗結束時并非對應信號自身周期結束，如下圖：

所以在上一個周期和下一個周期之間的信號過渡時有一個突變（周期邊界由黑色虛線表示）。這個突變導致信號中添加了除原始頻率之外的其他頻率。這就解釋了頻譜泄露現象。頻譜泄露會讓我們在計算頻域幅值時不夠準確，因為能量幅值被周圍的頻率平分。

了解了頻譜泄露的原因，我們可以通過調整采樣時長，保證在每個時間窗口中都有原始信號的整數個周期，例如設置采樣窗口 T = 2s（保證所有信號都是整數周期截斷）。然而通常我們可能未必能對真實數據進行延長采樣，可以考慮補零 zero-padding 的方式。

你可以用補零來對 DFT 進行插值。補零可以獲得更準確的幅值估計。補零并不能提高 DFT 的譜（頻率）分辨率。根據這種方法，將 DFT 要處理的信號用 0 補長到 2000 個點。有了這個長度，DFT 分倉之間的間距是 F_s=0.5Hz，在這種情況下，202.5Hz 正弦波的能量直接落在 DFT 分倉中。獲得 DFT 并繪制振幅估計值得到更準確的幅值。

時頻分析

在實際的應用中，信號經常是非平穩的，他們的頻域是隨時間在變化的，因此如果只用 FFT 無法反映出在什么時間里頻域上的這種變化。為了回答在什么時刻出現什么頻率的信號，我們需要進行時頻分析。例如，給定一段信號，有兩個 chirp 信號組成，第一個 chirp 信號在 0.1s 到 0.68s 一直存在，頻率隨時間 t 的函數是

第二個 chirp 信號在 0.1s 到 0.75s 一直存在，頻率隨時間t的函數是

我們希望看到在不同時間軸上都各自有哪些頻率出現。通過傅里葉變換（圖表6），我們確實可以得到信號中包含的頻域成分，但它回答不了各頻率成分在什么時候出現。所以首先想到的是對原始數據加窗，每個窗口內做傅里葉變換。這就是短時傅里葉變換 STFT。STFT 提供了信號頻率以及對應出現的時間信息。但是選擇窗口大小是個問題。在短時傅里葉變換時頻分析中，根據上面傅里葉變換的結論，選擇較短的窗口大小可以得到較好的時間分辨率，但頻率分辨率差。相反，選擇較大的窗口會獲得較好的頻率分辨率，但時間分辨率差。而且一旦選擇了窗口大小，它將在整個分析中保持不變。

如果可以提前知道待估計信號中我們想要的的頻率分量，則可以根據這些需求來選擇一個合適的窗口大小。兩個 chirp 信號的在初始時間點的瞬時頻率約為5Hz和15Hz。我們先選定窗口大小為 200ms（對應的頻率分辨率為 1/0.2s = 5Hz）。瞬時頻率在信號的前段被分辨出來，但在后段就不那么好了（后段頻率變化快，時間分辨率不夠會導致頻帶太寬）。

我們把窗口縮短為 50ms（對應的頻率分辨率 1/0.05s = 20Hz），雖然在后面高頻率差中可以分辨（時間分辨率更好），但在初始時刻低頻率差中分辨不開（頻率分辨率差）。

對于這種非平穩雙曲線 chirp 信號，STFT 很難找到合適的時間窗口大小使得對整個頻域的各頻率都能分辨的很好。連續小波變換（CWT）可以解決 STFT 固有的分辨率問題。連續小波變換得時間窗口是可變的。如下圖：

圖表10如果要分析的信號主要是緩慢震蕩的低頻信號，而高頻信號只是在一些瞬態時長或突變過程出現，這種情況可以考慮連續小波變換。如果高頻信號持續時間很長并且是信號中的主要成分，這種情況使用連續小波變換就不合適了。

繪制 CWT 的時頻圖。時頻圖顏色對應幅值，頻率軸用對數值，因為 CWT 中的頻率是對數的。從圖中可以清楚地看出信號中存在兩個雙曲 chirp 信號。使用連續小波變換，可以準確地估計整個信號周期在某個瞬時出現的頻率成分，而不必手動設置窗口長度。

有人可能會問白色虛線是什么？白色虛線是 COI（cone of influence （COI）），我們可以確定白色虛線以內的數據是準確的。在陰影區域的白線之外，時頻圖中的信息應被視為不準確的信息，因為可能存在邊緣效應。可以從時間分辨率的角度去看在低頻對應較大尺度小波從而時間分辨率較差。

詳細解釋請查看：Boundary Effects and the Cone of Influencehttps://www.mathworks.com/help/releases/R2020b/wavelet/ug/boundary-effects-and-the-cone-of-influence.html通過 3D 可視化查看小波幅值的增長速率，同時也可以將 CWT 的結果和真實瞬時頻率結果進行對比可以發現：CWT 結果和真實瞬時頻率一致性較好。

MATLAB 中提供了除了上述時頻變換，還有例如 Constant-Q Gabor Transform，Empirical Mode Decomposition and Hilbert-Huang Transform 等等。歡迎大家查看文檔搜索 https://ww2.mathworks.cn/help/。

振動分析

我們再回到振動場景本身，在進行時頻分析時，我們經常使用階次分析來確定發生在旋轉機械中的頻譜成分，從時域波形中追蹤和提取階次，計算不同階次的平均譜值，通過估計頻響函數、固有頻率、阻尼比和模態振型進行試驗模態分析，用時間同步平均法相干去除噪聲，用包絡譜分析磨損，為疲勞分析生成高周期雨流計數等等。

我們通過一個示例，對直升機在加速和減速過程機艙內的加速度傳感器數據進行階次分析從而確定振動源并進行優化來演示分析過程。當設備的轉速隨時間一直變化時，階次分析可以用來計算噪聲或振動。每個階數對應某個參考轉速的倍頻。例如，一個信號的頻率如果等于發動機轉頻的 2 倍，那階數就是 2。

我們通過仿真得到直升機在加速和減速過程機艙內的加速度傳感器數據。直升機有很多旋轉部件，包括發動機，齒輪箱，主旋翼和尾翼。每個部件都以一個和主發動機固定的速比運轉，每個部件都可能導致有害振動。示例中的信號是一個時域電壓信號 vib， 按 fs = 500Hz 的頻率進行采樣。

數據還包含渦輪發動機的轉速（rpm） vs 時間關系。通過數據可視化可以看到發動機轉速在爬升和滑行過程變化，同時振動加速度幅值也跟轉速有關。使用 RPM-Frequency 圖可視化數據為了能對振動信號進行時頻域的分析，我們可以使用 rpmfreqmap。這個函數計算信號的短時傅里葉變換生成 RPM-Frequency 圖譜。

rpmfreqmap 函數生成一幅時頻圖像同時還有對應的轉速時間圖像，還有下面的幾個數值參數。圖中幅值默認使用均方根（RMS）。當然也可以通過選項設置來使用其他幅值指標，例如峰值或功率值。瀑布圖按鈕可以生成三維的瀑布圖。

可以看到瀑布途中很多頻率對應的幅值會隨發動機轉速變化而變化。這說明這些頻率是發動機轉頻的階次頻率。在 RPM 峰值也對應著振動高幅值的成分，這些成分主要集中在 20-30Hz 之間。可以看到在當前默認的頻率分辨率（fs/128=3.906Hz，rpmfreqmap 函數默認將采樣率做 128 等分作為分辨率）不足以分辨一些轉速峰值處的低頻成分，因此，我們嘗試將頻率分辨率調小一些，設置為 1Hz。從結果可以看出1Hz分辨率的 RPM-frequency 圖可以清晰分辨出這些成分。

從現在的結果來看，在 RPM 峰值處對應的低頻成分可以被分辨出來，但同時我們也發現在轉速快速變化時出現了嚴重的頻率拖尾現象。這是因為在每個時間窗內隨著發動機轉速增加或變小，振動頻率階數也在變化，覆蓋一個比較大的范圍，從而導致一個較大的頻譜帶寬。

分辨率越高，要求時間窗采樣時長也越長，其中包含的頻率覆蓋范圍越廣，模糊現象越明顯。在直升機起飛或減速滑行階段，通過提升分辨率會導致頻率拖尾現象更嚴重。在這時，我么可以考慮使用階次圖來避免這種分辨率和包含頻率帶寬大小的矛盾。

使用 RPM 階次圖可視化數據函數 rpmordermap 生成階次頻譜圖與 RPM 圖，用于階次分析。由于每個階次是參考旋轉速度的固定倍數，因此階次圖包含一條條直的階次路徑，都是 RPM 的函數。函數 rpmordermap 與 rpmfreqmap 使用相同的參數，對應的頻域軸現在是階次，而不是頻率。使用 rpmfreqmap 可視化直升機數據的階次圖。指定階次分辨率為 0.005 倍基頻。

階次圖包含每個階次的直線路徑，每一階對應著發動機轉速的倍數對應的振動。階次圖可以清楚的看到每個頻率成分與發動機轉速的關系。與 RPM-頻率圖相比，頻率拖尾現象顯著被抑制。使用平均階次頻譜確定峰值階次接下來，確定階次圖的峰值位置。尋找主旋翼和尾翼階次的整數倍的階次，對應著這些旋翼產生的振動。函數 rpmordermap 返回包含各階次的階次圖和與時間對應的 RPM 值。通過分析數據以確定直升機艙內高振幅振動的對應階次。計算并返回階次圖譜數據。

接下來，使用 orderspectrum 計算并繪制 map 的平均譜。該函數接受 rpmordermap 生成的階次圖表作為輸入，并對時域取平均。

返回平均頻譜值，并調用 findpeaks 以返回兩個最高峰值的位置。

在圖中大約 0.05 階處可以看到兩個間隔很近的主峰。階次小于 1，因為振動頻率低于發動機轉速。分析峰值階次隨時間的變化接下來，使用 ordertrack 求峰值階次的幅值隨時間的變化。使用 map 作為輸入，通過不帶輸出參數調用 ordertrack 來繪制兩個峰值階次的振幅。

隨著發動機轉速的增加，兩個階次的振幅都會增加。接下來，使用 orderwaveform 提取每個階次對應的時域波形。提取出來的時域波形可以直接與原始振動信號進行比較。orderwaveform 使用 Vold-Kalman 濾波器提取特定階次的時域波形。將兩個峰值階次對應的時域波形之和與原始信號進行比較。

減少機艙振動為了確定機艙振動的來源，我們可以將振動峰值對應的階次和直升機旋轉部件的階次配合起來看。

可以看到最大振幅分量的頻率是主旋翼頻率的四倍。而主旋翼有四個葉片，很可能是這種振動的來源，通常對于每個旋翼有 N 葉片的直升機來說，主轉速的 N 階振動是很常見的。同樣，第二大分量位于尾旋翼速度的一階次處，表明振動可能源于尾旋翼。在對主旋翼和尾旋翼進行軌跡和平衡調整后，采集新數據集。加載新數據集并比較調整前后的階次頻譜。

主峰的振幅現在大幅降低。

總結

此示例使用階次分析來確定直升機的主旋翼和尾翼是否為機艙內高振幅振動的潛在來源。首先，使用了 rpmfreqmap 和 rpmordermap 對階次進行可視化。RPM-階次圖在整個 RPM 范圍內實現了階次分離，且消除了在 RPM-頻率圖中出現的頻率拖尾。

rpmordermap 最適合可視化在發動機加速和減速期間低 RPM 時的振動分量。接下來，該示例使用 orderspectrum 確定峰值階次，使用 ordertrack 可視化峰值階次的振幅隨時間的變化情況，使用 orderwaveform 提取峰值階次的時域波形。

最大的振幅振動分量出現在主旋翼旋轉頻率的四倍處，表明主旋翼葉片不平衡。第二大分量出現在尾旋翼的旋轉頻率處。調整旋翼后，振動程度得以降低。以后，我們還將陸續為大家介紹此系列的其他兩個主題：模態分析與預測性維護：人工智能算法應用。

編輯：jq