做EEG頻譜分析，看這一篇文章就夠了！

所謂頻譜分析，又稱為功率譜分析或者功率譜密度（Power Spectral Density, PSD）分析，實際就是通過一定方法求解訊號的功率power隨著頻率變化曲線。筆者在這裡對目前常用的頻譜分析方法做一個總結，並重點介紹目前EEG分析中最常用的頻譜分析方法，並給出相應的Matlab程式。

1.頻譜分析的方法有哪些？
目前來說，功率譜分析的方法大致可以分為兩大類：第一類是經典的功率譜計算方法，第二類是現代功率譜計算方法，如圖1所示。
其中第一類經典功率譜分析方法，又可以分為直接法、間接法和改進的直接法。直接法又稱之為週期圖法，簡單地說，其直接利用訊號的傅立葉變換系數的幅度平方來計算訊號的功率譜。間接法又稱為自相關函式法，其先估算出訊號的自相關函式，然後對自相關函式求傅立葉變換從而得到訊號的功率譜。改進的直接法，是針對直接法存在的缺點改進而來的方法，包括Barlett法、Welch法和Nuttall法。
第二類現代功率譜計算方法，又可以分為基於引數建模的功率譜計算和基於非引數建模的功率譜計算。基於引數建模的功率譜計算方法又分為基於AR模型、MA模型、ARMA模型等方法；基於非引數建模的功率譜計算方法主要基於矩陣特徵分解的功率譜估計，主要包括基於MUSIC演算法的功率譜估計和基於特徵向量的功率譜估計。
本文，筆者主要對經典功率譜分析方法中的直接法（週期圖法）以及在EEG頻譜分析中最常用的改進直接法（Welch法）進行詳細的介紹，並給出相應的Matlab程式。

圖1

2.直接法計算PSD
直接法又稱之為週期圖法，是由Schuster於1899年提出，其方法很簡單：對於長度為N的序列x(n)，其傅立葉變換為X(k),那麼x(n)在每個頻率（或k）處的功率可以表示為

即傅立葉變換系數的幅值平方除以資料長度N。
根據直接法求解PSD的定義，可以直接通過呼叫Matlab中的fft函式（fft函式是計算訊號的傅立葉變換）進行計算；
此外，Matlab中有專門的函式periodogram實現直接法的PSD計算。
例1：按照直接法計算PSD的定義，利用Matlab中的fft函式直接計算訊號y=5sin(2pif1t)+3cos(2pif2*t)+δ（其中f1=20Hz，f2=35Hz，δ為一隨機噪聲）的PSD。
結果如圖2所示，Matlab程式可以在公眾號後臺輸入“PSDcode”進行下載，下載後可以直接在Matlab中執行出以下結果。
Matlab中有專門的函式periodogram，也可以實現直接法的PSD計算，關於其用法，這裡筆者就不再贅述，各位可以直接在Matlab中help一下這個函式即可，其使用方法很簡單。

圖2

3.直接法計算PSD
上述直接法計算PSD，雖然可以直接FFT，計算速度快，但是頻率解析度比較低。因此，研究者提出了改進的直接法來計算PSD。**其中Welch方法就是其中的一種。Welch方法的思路是：**先把長度為N的訊號分成L段，每段資料長為M，則N=LM；然後把窗函式w加到每段資料上，求出每段資料的功率譜；最後對每段資料的功率譜進行平均，得到整個訊號的功率譜。各個資料段之間可以有重疊，窗函式w可以選擇如Hanning、Hamming等任意一種視窗。
Welch方法是EEG的PSD計算中最常用的一種方法，在Matlab中直接採用pwelch函式可以實現Welch方法的功率譜估計，其一般呼叫格式如下：
[Pxx,F] = pwelch(X,WINDOW,NOVERLAP,NFFT,Fs)
[Pxx,F] = pwelch(X,WINDOW,NOVERLAP,NFFT,Fs,REQRANGE)
關於函式中的引數含義，各位可以在Matlab命令視窗中輸入“help pwelch”即可檢視其詳細用法，這裡筆者不再贅述。
例2：採用Welch方法計算訊號y=5sin(2pif1t)+3cos(2pif2*t)+δ（其中f1=20Hz，f2=35Hz，δ為一隨機噪聲）的PSD。
計算結果如圖3所示，與圖2利用直接法求得的PSD基本一樣；Matlab程式可以在公眾號後臺輸入“PSDcode”進行下載，下載後可以直接在Matlab中執行出以下結果。

圖3

4.總結
本文首先對目前進行PSD計算的不同方法進行了總結和簡單介紹，重點詳細介紹瞭如何利用直接法和改進的直接法（Welch法）來計算訊號的PSD，並給出了Matlab程式。其中Welch法是目前計算EEG的PSD最常用的方法之一，理解和學會使用Welch法進行頻率分析對於我們做EEG研究來說至關重要。

本文來自**河南悅影醫藥科技有限公司（簡稱“悅影科技”）**旗下以腦科學和神經科學為主題的公益性自媒體平臺——“腦之說”微信公眾號，歡迎各位朋友搜尋關注“腦之說”公眾號。本文為悅影科技原創技術文章，轉載請聯絡趙老師（微訊號：kervin_zhao；郵箱：yueyingkj@163.com）。