功率譜密度的相關推導以及Python實現

本文主要介紹了離散信號功率譜密度的相關推導以及$Python$實現。特別是，很多教材默認采樣頻率為單位1，本文不做此默認相關推導更具一般性。文章內容安排如下：第一部分介紹基本概念和相關推導；第二部分分別利用現成的$matplotlib.pyplot.psd$庫和$numpy.fft.fft$庫計算離散信號的功率譜密度并驗證結果。湍流領域中的文獻常用預乘譜，其物理解釋可以參考這個網站。功率譜密度的現代估計方法可以參考這個網站。

一、理論推導

1.1 基本概念

對于一個連續隨機信號$x_c(t)$，由帕薩維爾等式可知：
$${\int }_{?\infty }^{+\infty }{x}_c(t{)}^{2}dt=\frac{1}{2\pi }{\int }_{?\infty }^{+\infty }|{\hat{x}}_c(\omega ){|}^{2}d\omega \text{\hspace{1em}(1)}$$ 其中，${\hat{x}}_c(\omega )$為隨機信號$x_c(t)$的傅里葉變換。類似于電磁學中用電流信號的平方表示功率($P=I^{2}R$)，上式左端被積函數$W=x_c(t{)}^2$表示隨機信號的功率。因此，(1)式左端表示隨機信號的總能量，右端被積函數$E(\omega )=|{\hat{x}}_c(\omega ){|}^2$表示單位頻率的能量，也即連續信號的能量譜密度。對于離散隨機信號$x(n)=x_c(n\Delta t),\hat{x}(k)={\hat{x}}_c({\omega }_k)$，有離散形式的帕薩維爾等式：
$$\sum _{n=0}^{N?1}{x}_c(n\Delta t{)}^2\Delta t=\frac{1}{2\pi }\sum _{k=0}^{N?1}|{\hat{x}}_c({\omega }_k){|}^2\Delta \omega \text{\hspace{1em}(2)}$$$${\omega }_k=2\pi f_k=\frac{2k\pi }{T}=\frac{2k\pi }{N\Delta t},\Delta \omega =2\pi \Delta f=\frac{2\pi }{N\Delta t}.$$其中，$\Delta t$為采樣周期，相應的采樣頻率為：$f_s=1/\Delta t.$ 將上式整理得到：
$$\sum _{n=0}^{N?1}x(n{)}^2\Delta t=\sum _{k=0}^{N?1}\frac{|f_s?\hat{x}(k){|}^2}{f_s^2}?\Delta f\text{\hspace{1em}(3)}$$式中${\hat{x}}_s(k)=f_s?\hat{x}(k)$是離散傅里葉變換的結果(推導見下節內容)。于是上式可以寫成：
$$\sum _{n=0}^{N?1}x(n{)}^2\Delta t=\sum _{k=0}^{N?1}\frac{|{\hat{x}}_s(k){|}^2}{f_s^2}?\Delta f\text{\hspace{1em}(4)}$$對于離散信號，我們可以定義離散信號的能量譜密度為離散信號離散傅里葉變換的模的平方除以采樣頻率的平方，即：$E(k)=|{\hat{x}}_s(k){|}^2/f_s^2.$ 此時，離散信號的平均功率為：
$$\overline{W}=\frac{1}{N}\sum _{n=0}^{N?1}x(n{)}^2=\sum _{k=0}^{N?1}\frac{|{\hat{x}}_s(k){|}^2}{N{f}_s}?\Delta f.\text{\hspace{1em}(5)}$$同理，可以定義離散信號的功率譜密度為離散信號的能量譜密度除以總時間T，即：$P(k)=\frac{E(k)}{N\Delta t}=|{\hat{x}}_s(k){|}^2/(f_{s}N)$。一般來說，我們需要對原始信號加窗$w(n)$。加窗后，為使其平均功率保持不變需要乘上一個能量恢復系數$K$[1]，也即有：$$\frac{1}{N}\sum _{n=0}^{N?1}x(n{)}^2=K?\frac{1}{N}\sum _{n=0}^{N?1}w(n{)}^{2}x(n{)}^2.$$為計算方便，我們令$x(n)=1$得到加函數$w(n)$窗的恢復系數為：$$K=\frac{N}{{\sum }_{n=0}^{N?1}w(n{)}^2}=\frac{N}{||w|{|}^2}.$$加窗以后的功率譜密度為：$$P(k)=K?\frac{E(k)}{N\Delta t}=\frac{|{\hat{x}}_s(k){|}^2}{f_s||w|{|}^2}.$$

1.2 離散傅里葉變換

對連續信號$x_c(t)$進行采樣得到離散信號$x(t)$的過程可以描述為：
$$x(t)=x_c(t)s(t),s(t)=\sum _{n=0}^{N?1}\delta (t?n\Delta t),\delta (t?t_0)=\{\begin{array}{rlr}1,& & t=t_0\\ 0,& & t\ne t_0\end{array}\text{\hspace{1em}(6)}$$其中，$s(t)$稱之為采樣函數。根據連續函數的傅里葉變換可知：
$${\hat{x}}_c(\omega )={\int }_{?\infty }^{+\infty }{x}_c(t)e^{?i\omega t}dt.\text{\hspace{1em}(7)}$$將隨機信號$x(t)$代入上式即可得到離散時間傅里葉變換如下：
$$\hat{x}(\omega )={\int }_{?\infty }^{+\infty }{x}_c(t)s(t)e^{?i\omega t}dt=\sum _{n=0}^{N?1}{\int }_{?\infty }^{+\infty }{x}_c(t)\delta (t?n\Delta t)e^{?i\omega t}dt.\text{\hspace{1em}(8)}$$將上述積分寫成離散形式可得：
$$\hat{x}(\omega )=\sum _{n=1}^{N?1}\sum _{k=0}^{N?1}{x}_c(t)\delta (k\Delta t?n\Delta t)e^{?i\omega k\Delta t}\Delta t=\sum _{n=1}^{N?1}{x}_c(n\Delta t)e^{?i\omega n\Delta t}\Delta t.\text{\hspace{1em}(9)}$$將頻率離散后(見上一節的處理)代入上述方程即可得到離散傅里葉變換如下：
$$\hat{x}({\omega }_k)=\sum _{n=0}^{N?1}{x}_c(n\Delta t)e^{?i{\omega }_{k}n\Delta t}\Delta t.\text{\hspace{1em}(10)}$$整理得到：
$${\hat{x}}_s(k)=f_s?\hat{x}(k)=\sum _{n=0}^{N?1}x(n)e^{?i\frac{2\pi kn}{N}}.\text{\hspace{1em}(11)}$$等式左端${\hat{x}}_s(k)$即為離散傅里葉變換的結果。

二、python代碼

導入基本庫

import numpy as np from numpy.fft import fft from numpy.linalg import norm import matplotlib.pyplot as plt

離散信號的生成

pi = np.pi # 圓周率 dt = 0.1 # 采樣周期 fs = 1/dt # 采樣頻率f1 = 0.05 # 信號的特征頻率1 f2 = 0.40 # 信號的特征頻率2 f3 = 1.0 # 信號的特征頻率3N = 2**12 # 離散信號的長度 tn = np.arange(0,dt*N,dt) # 時間序列 x = 2*np.cos(2*pi*f1*tn)+2*np.cos(2*pi*f2*tn)+2*np.cos(2*pi*f3*tn) # 生成離散信號 x = x*(1+0.1*np.random.randn(N)) # 加入隨機噪聲 nfft = 256 # psd的窗長# 將離散信號輸出 plt.figure(figsize=(22,2)) plt.plot(tn,x,'c-',linewidth=0.5) plt.xlim(0,dt*N) plt.show()

輸出結果如下圖所示：

利用plt.psd計算功率譜密度

[Pxx1,f1] = plt.psd(x, # 隨機信號NFFT=nfft, # 每個窗的長度Fs=fs, # 采樣頻率detrend='mean', # 去掉均值window=np.hanning(nfft), # 加漢尼窗noverlap=int(nfft*3/4), # 每個窗重疊75%的數據sides='twosided') # 求雙邊譜 plt.xscale('log')

輸出結果如下圖所示：

利用numpy庫中的fft計算功率譜密度

Pxx = [] for i in range(N//nfft): # 窗與窗之間數據不重疊w = np.hanning(nfft) # 加漢尼窗p = np.abs(fft(w*x[i*nfft:(i+1)*nfft]))**2/fs/norm(w)**2 # 計算功率譜Pxx.append(p) # 每個窗的結果 Pxx2 = np.mean(Pxx,axis=0) # 將所有窗平均得到最終結果f2 = np.array([k*fs/nfft for k in range(nfft//2)]) # 相應的頻率：fk

將結果可視化并與plt.psd的結果進行比較。

# 雙邊譜只有半邊有效(采樣定理) plt.plot(f2,Pxx2[:nfft//2],color='orange',linestyle='-',label='fft') # 利用fft計算的結果 plt.plot(f1[nfft//2:],Pxx1[nfft//2:],'b--',label='plt.psd') # 利用plt.psd計算的結果plt.xscale('log') plt.yscale('log') plt.xlabel('f/Hz',fontsize=12) plt.ylabel('PSD', fontsize=12) plt.legend(frameon=False,fontsize=12) plt.tick_params(top=True,right=True,direction='in',which='both') plt.show()

結果驗證
按照前面的定義可知，功率譜密度與頻率圍成的面積表示離散隨機信號的平均功率。

print('平均功率：',sum(x**2)/N) # 輸出離散隨機信號的的平均功率 print(' plt.psd：',sum(Pxx1)*fs/nfft) # 輸出plt.psd計算得到的平均功率 print(' fft：',sum(Pxx2)*fs/nfft) # 輸出fft計算得到的平均功率

輸出結果為：

若有錯誤歡迎在下方評論指出，謝謝！
[1]焦新濤,丁康.加窗頻譜分析的恢復系數及其求法[J].汕頭大學學報(自然科學版),2003(03):26-30+38.

總結

以上是生活随笔為你收集整理的功率谱密度的相关推导以及Python实现的全部內容，希望文章能夠幫你解決所遇到的問題。

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