互譜密度

理論

先講一點兒理論相關的東西，在後面使用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），這個定理表明函數平方的和（或積分），也就是其能量，等于其傅裡葉轉換式平方之和（或者積分）：

∫+∞−∞∣x(t)∣2dt=∫+∞−∞∣x(f)∣2df

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

互譜密度函數定義

互譜密度函數的定義，數學上可描述為

Sxy​(f)=T−>∞lim​T1​[XT∗​(f)YT​(f)],−∞<f<+∞(1)

由于互譜密度函數的推導方法與自譜密度函數相同，它們的差别隻是

是信号x(t)的自乘，而

是信号

的互乘。應當注意的是，因為

​    

與

​    

一般不是互為共轭，是以

為複數性質。

互譜的單邊譜為

由式(2)互譜密度函數也可以描述為

總結

自功率譜密度函數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）Pxy=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 noise

       y = cos(2*pi*t*100)+randn(size(t));  % A cosine of 100Hz plus noise

       cpsd(x,y,[],[],[],Fs,'twosided');    % Uses default window, overlap & NFFT. 

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

--------------------- 

作者：Big_quant 

來源：CSDN 

原文：https://blog.csdn.net/lvsehaiyang1993/article/details/84348788 

版權聲明：本文為部落客原創文章，轉載請附上博文連結！