簡述MATLAB中傳統(tǒng)信號(hào)如何處理

振動(dòng)分析在工程領(lǐng)域很常見，汽車中的 NVH，風(fēng)機(jī)的模態(tài)分析，塔架、葉片等機(jī)械件的振動(dòng)頻率分析，發(fā)電機(jī)軸承振動(dòng)，飛機(jī)發(fā)動(dòng)機(jī)階次分析等等。經(jīng)常提到的知識(shí)點(diǎn)像信號(hào)處理，F(xiàn)FT，階次分析，頻譜估計(jì)，包絡(luò)譜分析，模態(tài)分析，頻率響應(yīng)函數(shù)估計(jì)等等。

傅里葉變換

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

容易證明

將時(shí)域信號(hào)投影到這組正交基上。傅里葉變換也就可以看作 f（t） 在 k 對應(yīng)的基 上的投影。對于離散傅里葉變換（DFT）

為了方便理解，我們先通過一個(gè)例子，利用 fft 函數(shù)（快速傅里葉變換，實(shí)現(xiàn)DFT的一種快速算法）將一段矩形信號(hào)（圖中深藍(lán)色信號(hào)）進(jìn)行變換。通常我們得到幅值頻率圖像，可以看到對應(yīng)頻率的幅值。

為了更好的理解，我們把離散頻率對應(yīng)的正交基向量（圖淺綠色）也繪制出來，對應(yīng)圖中淺綠色。方塊綠則為幅值頻率值，此處我們首先只選了前 20 個(gè)頻率對應(yīng)的基向量進(jìn)行信號(hào)重構(gòu)，得到圖中紅色信號(hào)。接下來，我們增加用于信號(hào)重建的基的數(shù)量，分別我們?nèi)∏?60，前 300，發(fā)現(xiàn)信號(hào)越來越接近原始信號(hào)。

圖表2

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

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

我們從上面推導(dǎo)可以知道頻率分辨率由采樣時(shí)長決定，如果采樣時(shí)長太短，DFT 沒有能力區(qū)分兩個(gè)接近的頻率值。

這里也給大家一個(gè)小練習(xí)，例如可以合成一段信號(hào)，由 90Hz 和 100Hz 兩個(gè)頻率組成，我們只選取 0.05s 長的數(shù)據(jù)：對應(yīng) DFT 的頻率分辨率為 20Hz 》 （100Hz – 90Hz），然后進(jìn)行 fft 變換，查看是否在頻域上能夠區(qū)分開這兩個(gè)信號(hào)成分。

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

圖表4

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

從時(shí)域角度解釋：我們?nèi)×舜翱跁r(shí)間 1s 長度的數(shù)據(jù)，對應(yīng) f1 的 100 個(gè)完整周期以及 f2 的 202.5 個(gè)周期，DFT 默認(rèn)信號(hào)是采樣時(shí)間窗內(nèi)周期出現(xiàn)的無限長信號(hào)，因此當(dāng)窗重復(fù)時(shí)對于f1來說正好是整數(shù)拼接，與原信號(hào)完全一致，而 f2 在上一個(gè) 1s 采樣時(shí)間窗結(jié)束時(shí)并非對應(yīng)信號(hào)自身周期結(jié)束，如下圖：

所以在上一個(gè)周期和下一個(gè)周期之間的信號(hào)過渡時(shí)有一個(gè)突變（周期邊界由黑色虛線表示）。這個(gè)突變導(dǎo)致信號(hào)中添加了除原始頻率之外的其他頻率。這就解釋了頻譜泄露現(xiàn)象。頻譜泄露會(huì)讓我們在計(jì)算頻域幅值時(shí)不夠準(zhǔn)確，因?yàn)槟芰糠当恢車念l率平分。

了解了頻譜泄露的原因，我們可以通過調(diào)整采樣時(shí)長，保證在每個(gè)時(shí)間窗口中都有原始信號(hào)的整數(shù)個(gè)周期，例如設(shè)置采樣窗口 T = 2s（保證所有信號(hào)都是整數(shù)周期截?cái)啵Ｈ欢ǔＮ覀兛赡芪幢啬軐φ鎸?shí)數(shù)據(jù)進(jìn)行延長采樣，可以考慮補(bǔ)零 zero-padding 的方式。

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

時(shí)頻分析

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

第二個(gè) chirp 信號(hào)在 0.1s 到 0.75s 一直存在，頻率隨時(shí)間t的函數(shù)是

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

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

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

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

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

繪制 CWT 的時(shí)頻圖。時(shí)頻圖顏色對應(yīng)幅值，頻率軸用對數(shù)值，因?yàn)?CWT 中的頻率是對數(shù)的。從圖中可以清楚地看出信號(hào)中存在兩個(gè)雙曲 chirp 信號(hào)。使用連續(xù)小波變換，可以準(zhǔn)確地估計(jì)整個(gè)信號(hào)周期在某個(gè)瞬時(shí)出現(xiàn)的頻率成分，而不必手動(dòng)設(shè)置窗口長度。

有人可能會(huì)問白色虛線是什么？白色虛線是 COI（cone of influence （COI）），我們可以確定白色虛線以內(nèi)的數(shù)據(jù)是準(zhǔn)確的。在陰影區(qū)域的白線之外，時(shí)頻圖中的信息應(yīng)被視為不準(zhǔn)確的信息，因?yàn)榭赡艽嬖谶吘壭?yīng)?？梢詮臅r(shí)間分辨率的角度去看在低頻對應(yīng)較大尺度小波從而時(shí)間分辨率較差。

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

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

振動(dòng)分析

我們再回到振動(dòng)場景本身，在進(jìn)行時(shí)頻分析時(shí)，我們經(jīng)常使用階次分析來確定發(fā)生在旋轉(zhuǎn)機(jī)械中的頻譜成分，從時(shí)域波形中追蹤和提取階次，計(jì)算不同階次的平均譜值，通過估計(jì)頻響函數(shù)、固有頻率、阻尼比和模態(tài)振型進(jìn)行試驗(yàn)?zāi)B(tài)分析，用時(shí)間同步平均法相干去除噪聲，用包絡(luò)譜分析磨損，為疲勞分析生成高周期雨流計(jì)數(shù)等等。

我們通過一個(gè)示例，對直升機(jī)在加速和減速過程機(jī)艙內(nèi)的加速度傳感器數(shù)據(jù)進(jìn)行階次分析從而確定振動(dòng)源并進(jìn)行優(yōu)化來演示分析過程。當(dāng)設(shè)備的轉(zhuǎn)速隨時(shí)間一直變化時(shí)，階次分析可以用來計(jì)算噪聲或振動(dòng)。每個(gè)階數(shù)對應(yīng)某個(gè)參考轉(zhuǎn)速的倍頻。例如，一個(gè)信號(hào)的頻率如果等于發(fā)動(dòng)機(jī)轉(zhuǎn)頻的 2 倍，那階數(shù)就是 2。

我們通過仿真得到直升機(jī)在加速和減速過程機(jī)艙內(nèi)的加速度傳感器數(shù)據(jù)。直升機(jī)有很多旋轉(zhuǎn)部件，包括發(fā)動(dòng)機(jī)，齒輪箱，主旋翼和尾翼。每個(gè)部件都以一個(gè)和主發(fā)動(dòng)機(jī)固定的速比運(yùn)轉(zhuǎn)，每個(gè)部件都可能導(dǎo)致有害振動(dòng)。示例中的信號(hào)是一個(gè)時(shí)域電壓信號(hào) vib， 按 fs = 500Hz 的頻率進(jìn)行采樣。

數(shù)據(jù)還包含渦輪發(fā)動(dòng)機(jī)的轉(zhuǎn)速（rpm） vs 時(shí)間關(guān)系。通過數(shù)據(jù)可視化可以看到發(fā)動(dòng)機(jī)轉(zhuǎn)速在爬升和滑行過程變化，同時(shí)振動(dòng)加速度幅值也跟轉(zhuǎn)速有關(guān)。使用 RPM-Frequency 圖可視化數(shù)據(jù)為了能對振動(dòng)信號(hào)進(jìn)行時(shí)頻域的分析，我們可以使用 rpmfreqmap。這個(gè)函數(shù)計(jì)算信號(hào)的短時(shí)傅里葉變換生成 RPM-Frequency 圖譜。

rpmfreqmap 函數(shù)生成一幅時(shí)頻圖像同時(shí)還有對應(yīng)的轉(zhuǎn)速時(shí)間圖像，還有下面的幾個(gè)數(shù)值參數(shù)。圖中幅值默認(rèn)使用均方根（RMS）。當(dāng)然也可以通過選項(xiàng)設(shè)置來使用其他幅值指標(biāo)，例如峰值或功率值。瀑布圖按鈕可以生成三維的瀑布圖。

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

從現(xiàn)在的結(jié)果來看，在 RPM 峰值處對應(yīng)的低頻成分可以被分辨出來，但同時(shí)我們也發(fā)現(xiàn)在轉(zhuǎn)速快速變化時(shí)出現(xiàn)了嚴(yán)重的頻率拖尾現(xiàn)象。這是因?yàn)樵诿總€(gè)時(shí)間窗內(nèi)隨著發(fā)動(dòng)機(jī)轉(zhuǎn)速增加或變小，振動(dòng)頻率階數(shù)也在變化，覆蓋一個(gè)比較大的范圍，從而導(dǎo)致一個(gè)較大的頻譜帶寬。

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

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

階次圖包含每個(gè)階次的直線路徑，每一階對應(yīng)著發(fā)動(dòng)機(jī)轉(zhuǎn)速的倍數(shù)對應(yīng)的振動(dòng)。階次圖可以清楚的看到每個(gè)頻率成分與發(fā)動(dòng)機(jī)轉(zhuǎn)速的關(guān)系。與 RPM-頻率圖相比，頻率拖尾現(xiàn)象顯著被抑制。使用平均階次頻譜確定峰值階次接下來，確定階次圖的峰值位置。尋找主旋翼和尾翼階次的整數(shù)倍的階次，對應(yīng)著這些旋翼產(chǎn)生的振動(dòng)。函數(shù) rpmordermap 返回包含各階次的階次圖和與時(shí)間對應(yīng)的 RPM 值。通過分析數(shù)據(jù)以確定直升機(jī)艙內(nèi)高振幅振動(dòng)的對應(yīng)階次。計(jì)算并返回階次圖譜數(shù)據(jù)。

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

返回平均頻譜值，并調(diào)用 findpeaks 以返回兩個(gè)最高峰值的位置。

在圖中大約 0.05 階處可以看到兩個(gè)間隔很近的主峰。階次小于 1，因?yàn)檎駝?dòng)頻率低于發(fā)動(dòng)機(jī)轉(zhuǎn)速。分析峰值階次隨時(shí)間的變化接下來，使用 ordertrack 求峰值階次的幅值隨時(shí)間的變化。使用 map 作為輸入，通過不帶輸出參數(shù)調(diào)用 ordertrack 來繪制兩個(gè)峰值階次的振幅。

隨著發(fā)動(dòng)機(jī)轉(zhuǎn)速的增加，兩個(gè)階次的振幅都會(huì)增加。接下來，使用 orderwaveform 提取每個(gè)階次對應(yīng)的時(shí)域波形。提取出來的時(shí)域波形可以直接與原始振動(dòng)信號(hào)進(jìn)行比較。orderwaveform 使用 Vold-Kalman 濾波器提取特定階次的時(shí)域波形。將兩個(gè)峰值階次對應(yīng)的時(shí)域波形之和與原始信號(hào)進(jìn)行比較。

減少機(jī)艙振動(dòng)為了確定機(jī)艙振動(dòng)的來源，我們可以將振動(dòng)峰值對應(yīng)的階次和直升機(jī)旋轉(zhuǎn)部件的階次配合起來看。

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

主峰的振幅現(xiàn)在大幅降低。

總結(jié)

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

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

最大的振幅振動(dòng)分量出現(xiàn)在主旋翼旋轉(zhuǎn)頻率的四倍處，表明主旋翼葉片不平衡。第二大分量出現(xiàn)在尾旋翼的旋轉(zhuǎn)頻率處。調(diào)整旋翼后，振動(dòng)程度得以降低。以后，我們還將陸續(xù)為大家介紹此系列的其他兩個(gè)主題：模態(tài)分析與預(yù)測性維護(hù)：人工智能算法應(yīng)用。

編輯：jq