理論

先講一點兒理論相關的東西，在后面使用matlab實現，加深大家的理解。

功率譜密度

在物理學中，信號通常是波的形式表示，例如電磁波、隨機振動或者聲波。當波的功率頻譜密度乘以一個適當的系數后將得到每單位頻率波攜帶的功率，這被稱為信號的功率譜密度（power spectral density, PSD）；不要和 spectral power distribution（SPD） 混淆。功率譜密度的單位通常用每赫茲的瓦特數（W/Hz）表示，后者使用波長而不是頻率，即每納米的瓦特數（W/nm）來表示。

計算定義

盡管并非一定要為信號或者它的變量賦予一定的物理量綱，下面的討論中假設信號在時域內變化。
上面能量譜密度的定義要求信號的傅里葉變換必須存在，也就是說信號平方可積或者平方可加。一個經常更加有用的替換表示是功率譜密度（PSD），它定義了信號或者時間序列的功率如何隨頻率分布。這里功率可能是實際物理上的功率，或者更經常便于表示抽象的信號被定義為信號數值的平方，也就是當信號的負載為1歐姆(ohm)時的實際功率。此瞬時功率（平均功率的中間值）可表示為：

$$P=S^2(t)$$

由于平均值不為零的信號不是平方可積的，所以在這種情況下就沒有傅里葉變換。幸運的是維納-辛欽定理（Wiener-Khinchin theorem）提供了一個簡單的替換方法，如果信號可以看作是平穩隨機過程，那么功率譜密度就是信號自相關函數的傅里葉變換。

計算方法

信號的功率譜密度當且僅當信號是廣義的平穩過程的時候才存在。如果信號不是平穩過程，那么自相關函數一定是兩個變量的函數，這樣就不存在功率譜密度，但是可以使用類似的技術估計時變譜密度。
f(t) 的譜密度和 f(t) 的自相關組成一個傅里葉變換對（對于功率譜密度和能量譜密度來說，使用著不同的自相關函數定義）。
通常使用傅里葉變換技術估計譜密度，但是也可以使用如Welch法（Welch’s method）和最大熵這樣的技術。
傅里葉分析的結果之一就是Parseval（帕塞瓦爾）定理（Parseval’s theorem，其有時也被稱為瑞利能量定理，Rayleigh’s energy theorem），這個定理表明函數平方的和（或積分），也就是其能量，等于其傅里葉轉換式平方之和（或者積分）：
$${\int }_{?\infty }^{+\infty }|x(t){|}^{2}dt={\int }_{?\infty }^{+\infty }|x(f){|}^{2}df$$
其中 X(f) = F.T. { x(t) } 為x(t) 的連續傅立葉變換，f 是 x 的頻率分量。
上面的定理在離散情況下也是成立的 (DTFT 和 DFT)。另外的一個結論是功率譜密度下總的功率與對應的總的平均信號功率相等，它是逐漸趨近于零的自相關函數。

其他定義

功率譜密度譜是一種概率統計方法，是對隨機變量均方值的量度。一般用于隨機振動分析，連續瞬態響應只能通過概率分布函數進行描述，即出現某水平響應所對應的概率。
功率譜密度的定義是單位頻帶內的“功率”（均方值）
功率譜密度是結構在隨機動態載荷激勵下響應的統計結果，是一條功率譜密度值—頻率值的關系曲線，其中功率譜密度可以是位移功率譜密度、速度功率譜密度、加速度功率譜密度、力功率譜密度等形式。數學上，功率譜密度值—頻率值的關系曲線下的面積就是均方值
，當均值為零時均方值等于方差，即響應標準偏差的平方值。

互譜定義

互譜(cross-power spectrum) 互功率密度譜的簡稱，在頻域內描述兩個不同信號之間統計相關程度的一種方法。
設有兩個平穩隨機信號x(t)與y(t)，根據隨機過程理論，它們之間的統計相關特性，應該用其互相關函數表達。對x(t)與y(t)的互相關函數進行傅里葉變換，獲得其頻域中的功率密度譜，即稱為互功率密度譜，也稱互頻譜。可見，互譜與互相關函數是分別從頻域和時域描述兩個信號統計相關的兩種不同表示，它們互為傅里葉變換。互譜也適用于確定性信號分析。互譜在通信及信號處理領域中有重要用途，可用來測定一個未知參數的線性系統的頻率響應。這時主要要測出系統輸入和輸出信號之間的互譜。互譜也可以用于系統時延，如聲納接收信號等時延估計。
在實用中，通常利用快速傅里葉變換來計算和測量互譜，這是因為實際要求提高測量運算速度而提出來的，已經生產了許多測量功率譜密度函數的儀器，它們也可以用于互譜的測量。
現代譜分析是20世紀70年代發展起來的新興學科。新的方法、新的軟件以及新的VLSI譜分析專用硬件不斷出現，并在越來越多的領域中獲得成功的應用。在現代譜分析理論和方法的研究中，ARMA參數模型方法和特征分析方法是人們十分關心的研究課題。在參數模型方法中對AR方法的研究比較成熟。雖然從理論上說ARMA方法應當具有比AR方法更優越的性能，但現在還沒有有效的ARMA譜估計方法。已經提出的方法不是性能不夠理想，就是計算太復雜，距實際應用相差甚遠。現代譜估計中的特征分析方法都有優越的譜分析性能。特征方法和子空間理論的研究在陣列信號處理中對提高方位估計的分辨力和估計精度均有重要意義。
由于實際信號的時變、非平穩等復雜情況，研究多種信號的韌性(robust)譜估計，以及研究時變譜估計是譜分析研究的一個重要方面。對于實際應用中經常遇到的如噪聲中單個或多個正弦信號的頻率估計，窄帶信號和寬帶信號譜分析，以及淹沒在噪聲中的信號的譜分析等問題也是今后研究的重要課題。
為了加速現代譜分析方法的實際應用，人們重視快速算法和有效的算法結構的研究，尤其是能夠用超大規模集成電路硬件實現的實時譜分析方法。多維譜估計和高階譜估計的研究受到重視，其中人們十分感興趣的研究課題是對雙譜和三譜估計的研究，可用于估計信號中的相位和描述時間序列的非線性特性等。

互譜密度函數

互譜密度函數(cross-spectral density function)是互相關函數的傅立葉變換。互譜密度函數一般與互相關函數具有同樣的應用，但它提供的結果是頻率的函數而不是時間的函數。這—事實大大開拓了使用范圍，因此在可以應用相關分析的工程問題中大大增加了互譜方法的應用。互譜密度函數是有重要用途的，頻譜分析中能用互譜的測量結果來識別動力系統的特性以及計算頻響函數的振幅比和相位角

互譜密度函數定義

互譜密度函數的定義，數學上可描述為
$$S_{xy}(f)=\underset{T?>\infty }{\lim }\frac{1}{T}[X_T^{?}(f)Y_T(f)],?\infty <f<+\infty (1)$$

由于互譜密度函數的推導方法與自譜密度函數相同，它們的差別只是$$S_{xx}$$是信號x(t)的自乘，而$$S_{xy}$$是信號$$x(t),y(t)$$的互乘。應當注意的是，因為$$X_T^{?}$$與$$Y_T$$一般不是互為共軛，所以$$S_{xy}$$為復數性質。
互譜的單邊譜為
$$G_{xy}(f)=2S_{xy}(f),f\ge 0,0,f<0$$
由式(2)互譜密度函數也可以描述為
$$S_{xy}(f)=\underset{T?>\infty }{\lim }\frac{2}{T}[X_T^{?}(f)Y_T(f)].$$

總結

自功率譜密度函數Sxx（f）：反映相關函數在時域內表達隨機信號自身與其他信號在不同時刻的內在聯系。
1、當隨機信號均值為零時，自相關函數和自功率譜密度函數互為傅立葉變換對。
2、自功率譜密度有明確的物理含義：當tao＝0時，Sxx（f）曲線與頻率軸f所包圍的面積就是信號的平均功率。另外，Sxx（f）還表明了信號的功率密度沿頻率軸的分布狀況，因此稱Sxx（f）為自功率譜密度函數。

Matlab 實現

cpsd 使用welch方法計算互功率譜密度

1Pxy=cpsd（X，Y）返回使用welch方法計算的兩個離散時間序列信號的互譜功率密度Pxy，
默認參數中，X和Y會被分成八段，每段重復50%，每段通過漢明窗加窗，同時計算八個加權的譜密度并做平均。可以使用“Help pwelch”和“Help cpsd” 來查看完整的細節
X和Y可以是向量或者二維的矩陣。如果它們都是矩陣，他們必須有相同的尺寸。CPSD會按照列的維度計算，即Pxy（：，n）=cpsd（X(:,n),Y(:,n)）.如果其中一個是矩陣另一個是向量的話，這個向量會被轉變成一個列向量，并且進行擴展。因此，輸入都必須有相同數目的列。
Pxy是單位頻率能量的譜分布。對于實信號，cpsd反饋的是單邊的互譜密度；對于復數信號，它返回的是雙邊的互譜密度。注意這里單邊的互譜密度包含了輸入信號的全部能量。
Pxy=CPSD（X，Y，window），當window是一個向量時，將x和y的每一列分成由Pxy = cpsd(X,Y,WINDOW), 每一塊兒的向量長度和窗口數相同的不重合的區域。如果窗口是一個整數，那么每一列會被分割成和window數一樣長的區段，每一個區段都會使用漢明窗進行加窗。如果窗口參數是空的，則會將x和y分割成八個區域，并使用漢明窗進行加窗。
Pxy=cpsd（x,y,wiodow,noverlap）使用noverlap參數 控制每個區段之間重復的樣本數。
Noverlap 必須是一個小于窗口長度的整數，如果窗口數是一個向量的話。如果是一個標量的話，則必須比窗口小。默認的Noverlap等于50%。
【Pxy，W】=cpsd（x,y,window,noverlap,nfft）確定了用了計算CPSD的FFT的點數。對于每一個實信號，Pxy的長度為（NFFT/2+1），如果NFFT是偶數。如果NFFT是奇數則長度為（NFFT+1）/2.對于復信號，Pxy長度始終為NFFT。如果NFFT被設置為空，則NFFT點數或者是256或者是比x和y大的最小的2的冪。 如果NFFT比每段的數據都長，則會對數據段進行補零操作。如果數據段比NFFT點數大，則會使用wrapped使得數據的長度等于NFFT。這個產生了正確的FFT點數，在NFFT比區段長度短時。W是記載了歸一化后psd計算頻率點的向量。W是單位角度每樣本數。對于實信號，W擴張至【0,pi】這個區間當NFFT點數是偶數的時候，【0，pi）當NFFT是奇數的時候。
對于復數信號，w擴張至【0，2pi）這個區間。
[Pxy,W]=cpsd(X,Y,Window,Noverlap,W)計算在向量w中的標準角度頻率存在的雙邊互譜密度，w必須至少有兩個元素。
[Pxy,F]=cpsd(x,Y,window,noverlap,nfft,Fs)返回作為實際頻率函數的互譜密度。Fs是用HZ表示的采樣頻率。如果Hz數沒寫的話，則默認是1HZ.
F是存儲Pxy計算點的頻率向量（單位為Hz）。對于實信號，F會擴張到間隔[0,Fs/2]當NFFT
是偶數的時候，擴張到[0,Fs/2),當NFFT是奇數的時候。對于復數信號，F
則擴張至[0,Fs).
[Pxy,F]=cpsd(X,Y,window,noverlap,F，Fs)計算存儲在頻率向量F中存放的點的互譜密度。F必須使用至少兩個元素來表示，單位是HZ
[…]=cpsd(…,Freqrange)返回在特定頻率范圍計算的互譜密度，使用的參數FREQrange限定的范圍。
‘onesided’-返回輸入的實信號x和y的單邊的互譜密度。如果NFFT是偶數，Pxy的長度為NFFT/2+，計算的間隔為[0,pi].如果NFFT是奇數，Pxy的長度則為(NFFT+1)/2和同時在頻率范圍[0,pi]計算,當Fs被設置時，間隔為[0,fs/2]或者[0,Fs/2).
‘twosided’-返回雙邊的互譜密度信號，對于實或者復數的輸入信號x和y。Pxy有長度NFFT同時在頻率范圍[0,2pi)內計算。當Fs被設置時，間隔為[0,Fs）.
‘centered’-返回中心化的雙邊互譜密度對于實或者復的輸入信號x和y。Pxy長度為NFFT同時在間隔(-pi,pi]進行計算。
如果不設置輸出參數，則畫出互譜密度（單位是分貝/每頻率）
例子

Fs = 1000; t = 0:1/Fs:.296;x = cos(2*pi*t*200)+randn(size(t)); % A cosine of 200Hz plus noisey = cos(2*pi*t*100)+randn(size(t)); % A cosine of 100Hz plus noisecpsd(x,y,[],[],[],Fs,'twosided'); % Uses default window, overlap & NFFT.

結果如圖所示，100Hz有個高峰。

總結

以上是生活随笔為你收集整理的谱分析：互谱密度的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。