功率谱密度的相关推导以及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\cdot \hat{x}(k){|}^2}{f_s^2}\cdot \Delta f\text{\hspace{1em}(3)}$$式中${\hat{x}}_s(k)=f_s\cdot \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}\cdot \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}\cdot \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\ast \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\ast \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\cdot \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代码

1. 导入基本库

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

2. 离散信号的生成

pi = np.pi  # 圆周率
dt = 0.1    # 采样周期
fs = 1/dt   # 采样频率

f1 = 0.05   # 信号的特征频率1
f2 = 0.40   # 信号的特征频率2
f3 = 1.0    # 信号的特征频率3

N  = 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()

输出结果如下图所示：

3. 利用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')

输出结果如下图所示：

4. 利用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()

5. 结果验证
   按照前面的定义可知，功率谱密度与频率围成的面积表示离散随机信号的平均功率。

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.