功率譜密度的相關(guān)推導(dǎo)以及Python實(shí)現(xiàn)

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

一、理論推導(dǎo)

1.1 基本概念

對(duì)于一個(gè)連續(xù)隨機(jī)信號(hào)$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 )$為隨機(jī)信號(hào)$x_c(t)$的傅里葉變換。類似于電磁學(xué)中用電流信號(hào)的平方表示功率($P=I^{2}R$)，上式左端被積函數(shù)$W=x_c(t{)}^2$表示隨機(jī)信號(hào)的功率。因此，(1)式左端表示隨機(jī)信號(hào)的總能量，右端被積函數(shù)$E(\omega )=|{\hat{x}}_c(\omega ){|}^2$表示單位頻率的能量，也即連續(xù)信號(hào)的能量譜密度。對(duì)于離散隨機(jī)信號(hào)$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$為采樣周期，相應(yīng)的采樣頻率為：$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)$是離散傅里葉變換的結(jié)果(推導(dǎo)見(jiàn)下節(jié)內(nèi)容)。于是上式可以寫(xiě)成：
$$\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)}$$對(duì)于離散信號(hào)，我們可以定義離散信號(hào)的能量譜密度為離散信號(hào)離散傅里葉變換的模的平方除以采樣頻率的平方，即：$E(k)=|{\hat{x}}_s(k){|}^2/f_s^2.$ 此時(shí)，離散信號(hào)的平均功率為：
$$\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)}$$同理，可以定義離散信號(hào)的功率譜密度為離散信號(hào)的能量譜密度除以總時(shí)間T，即：$P(k)=\frac{E(k)}{N\Delta t}=|{\hat{x}}_s(k){|}^2/(f_{s}N)$。一般來(lái)說(shuō)，我們需要對(duì)原始信號(hào)加窗$w(n)$。加窗后，為使其平均功率保持不變需要乘上一個(gè)能量恢復(fù)系數(shù)$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.$$為計(jì)算方便，我們令$x(n)=1$得到加函數(shù)$w(n)$窗的恢復(fù)系數(shù)為：$$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 離散傅里葉變換

對(duì)連續(xù)信號(hào)$x_c(t)$進(jìn)行采樣得到離散信號(hào)$x(t)$的過(guò)程可以描述為：
$$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)$稱之為采樣函數(shù)。根據(jù)連續(xù)函數(shù)的傅里葉變換可知：
$${\hat{x}}_c(\omega )={\int }_{?\infty }^{+\infty }{x}_c(t)e^{?i\omega t}dt.\text{\hspace{1em}(7)}$$將隨機(jī)信號(hào)$x(t)$代入上式即可得到離散時(shí)間傅里葉變換如下：
$$\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)}$$將上述積分寫(xiě)成離散形式可得：
$$\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)}$$將頻率離散后(見(jiàn)上一節(jié)的處理)代入上述方程即可得到離散傅里葉變換如下：
$$\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)$即為離散傅里葉變換的結(jié)果。

二、python代碼

導(dǎo)入基本庫(kù)

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

離散信號(hào)的生成

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

輸出結(jié)果如下圖所示：

利用plt.psd計(jì)算功率譜密度

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

輸出結(jié)果如下圖所示：

利用numpy庫(kù)中的fft計(jì)算功率譜密度

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

將結(jié)果可視化并與plt.psd的結(jié)果進(jìn)行比較。

# 雙邊譜只有半邊有效(采樣定理) plt.plot(f2,Pxx2[:nfft//2],color='orange',linestyle='-',label='fft') # 利用fft計(jì)算的結(jié)果 plt.plot(f1[nfft//2:],Pxx1[nfft//2:],'b--',label='plt.psd') # 利用plt.psd計(jì)算的結(jié)果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()

結(jié)果驗(yàn)證
按照前面的定義可知，功率譜密度與頻率圍成的面積表示離散隨機(jī)信號(hào)的平均功率。

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

輸出結(jié)果為：

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

總結(jié)

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

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。