使用scipy.signal.welch找不到合适的能量

Question

对于x(t)具有间距dt(等于1/fs,fs即采样率)的给定离散时间信号,能量为:

E[x(t)] = sum(abs(x)**2.0)/fs

然后我做了一个DFT x(t):

x_tf = np.fft.fftshift( np.fft.fft( x ) ) / ( fs * ( 2.0 * np.pi ) ** 0.5 )

并再次计算能量:

E[x_tf] = sum( abs( x_tf ) ** 2.0 ) * fs * 2 * np.pi / N

(这里因子fs*2*np.pi/N=脉动间距dk,文档fftfreq提供了有关频域间距的更多细节),我有相同的能量:

E[x(t)] = E[x_tf]

但是...当我计算的功率谱密度x(t)使用scipy.signal.welch,我无法找到合适的能量.scipy.signal.welch返回频率f和能量的矢量Pxx(或每个频率的能量,取决于scaling我们在参数中输入的能量scipy.signal.welch).

我怎样才能找到E[x(t)]和E[x_tf]使用相同的能量Pxx？我试着计算:

E_psd = sum(Pxx_den) / nperseg

其中nperseg是Welch算法的每个段的长度,因子fs和np.sqrt(2*np.pi)被抵消的因素,并且重新缩放E [x(t)] nperseg,但没有任何成功(数量级小于E[x(t)])

我使用以下代码生成我的信号:

#Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by 0.001 V**2/Hz of white noise sampled at 10 kHz.

fs = 10e3   #sampling rate, dt = 1/fs
N = 1e5
amp = 2*np.sqrt(2)
freq = 1234.0
noise_power = 0.001 * fs / 2
time = np.arange(N) / fs
x = amp*np.sin(2*np.pi*freq*time)
x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)

我做了以下工作来获得功率谱密度:

f, Pxx_den = signal.welch(x, fs )

Answer 1

解决这种明显的差异在于仔细理解和应用

• 连续与离散傅里叶变换,和
• 给定信号的能量,功率和功率谱密度

我也一直在努力解决这个问题,所以我将尽力在下面的讨论中尽可能明确.

离散傅立叶变换(DFT)

满足某些可积性条件的连续信号x(t)具有傅里叶变换X(f).然而,当使用离散信号x [n]时,通常使用离散时间傅里叶变换(DTFT).我将DTFT表示为X_ {dt}(f),其中dt等于相邻样本之间的时间间隔.回答问题的关键要求您认识到DTFT 不等于相应的傅里叶变换!事实上,这两者是相关的

X_ {dt}(f)=(1/dt)*X(f)

此外,离散傅立叶变换(DFT)仅仅是DTFT 的离散样本.当然,DFT是Python在使用时返回的np.fft.fft(...).因此,您的计算DFT 不等于傅里叶变换!

功率谱密度

scipy.signal.welch(..., scaling='density', ...)返回离散信号x [n] 的功率谱密度(PSD)的估计值.对PSD的全面讨论有点超出了本文的范围,但对于简单的周期性信号(例如在您的示例中),PSD S_ {xx}(f)给出为

S_ {xx} = | X(f)| ^ 2/T.

其中| X(f)| 是信号的傅立叶变换,T是信号的总持续时间(时间)(如果你的信号x(t)是一个随机过程,我们必须在系统的许多实现中采用整体平均值. ..).信号中的总功率只是S_ {xx}在系统频率带宽上的积分.使用上面的代码,我们可以写

import scipy.signal

# Estimate PSD `S_xx_welch` at discrete frequencies `f_welch`
f_welch, S_xx_welch = scipy.signal.welch(x, fs=fs)

# Integrate PSD over spectral bandwidth
# to obtain signal power `P_welch`
df_welch = f_welch[1] - f_welch[0]
P_welch = np.sum(S_xx_welch) * df_welch

要与您的np.fft.fft(...)计算(返回DFT)联系,我们必须使用上一节中的信息,即

X [k] = X_ {dt}(f_k)=(1/dt)*X(f_k)

因此,为了从FFT计算中计算功率谱密度(或总功率),我们需要认识到这一点

S_ {xx} = | X [k] | ^ 2*(dt ^ 2)/ T.

# Compute DFT
Xk = np.fft.fft(x)

# Compute corresponding frequencies
dt = time[1] - time[0]
f_fft = np.fft.fftfreq(len(x), d=dt)

# Estimate PSD `S_xx_fft` at discrete frequencies `f_fft`
T = time[-1] - time[0]
S_xx_fft = ((np.abs(Xk) * dt) ** 2) / T

# Integrate PSD over spectral bandwidth to obtain signal power `P_fft`
df_fft = f_fft[1] - f_fft[0]
P_fft = np.sum(S_xx_fft) * df_fft

您P_welch和和的值P_fft应该非常接近,并且接近信号中的预期功率,可以计算为

# Power in sinusoidal signal is simply squared RMS, and
# the RMS of a sinusoid is the amplitude divided by sqrt(2).
# Thus, the sinusoidal contribution to expected power is
P_exp = (amp / np.sqrt(2)) ** 2 

# For white noise, as is considered in this example,
# the noise is simply the noise PSD (a constant)
# times the system bandwidth. This was already
# computed in the problem statement and is given
# as `noise_power`. Simply add to `P_exp` to get
# total expected signal power.
P_exp += noise_power

注意: P_welch并且P_fft不会完全相等,甚至可能在数值精度内不相等.这可归因于存在与功率谱密度的估计相关联的随机误差的事实.为了减少这些错误,Welch的方法将您的信号分成几个段(其大小由nperseg关键字控制),计算每个段的PSD,并平均PSD以获得更好的信号PSD估计(平均的段数越多,产生的随机误差就越小.实际上,FFT方法仅相当于对一个大段进行计算和平均.因此,我们期望P_welch和之间存在一些差异P_fft,但我们应该期望P_welch更准确.

信号能量

正如您所说,信号能量可以从Parseval定理的离散版本中获得

# Energy obtained via "integrating" over time
E = np.sum(x ** 2)

# Energy obtained via "integrating" DFT components over frequency.
# The fact that `E` = `E_fft` is the statement of 
# the discrete version of Parseval's theorem.
N = len(x)
E_fft = np.sum(np.abs(Xk) ** 2) / N

我们现在想要了解S_xx_welch上面计算的如何与信号中scipy.signal.welch(...)的总能量相关E.从上面,S_xx_fft = ((np.abs(Xk) * dt) ** 2) / T.重新排列此表达式中的术语,我们看到了这一点np.abs(Xk) ** 2 = (T / (dt ** 2)) * S_xx_fft.进一步,

从上面,我们知道,np.sum(S_xx_fft) = P_fft / df_fft那P_fft和P_welch大致相等.而且,P_welch = np.sum(S_xx_welch) / df_welch这样我们就可以获得

np.sum(S_xx_fft) = (df_welch / df_fft) * np.sum(S_xx_welch)

而且,S_xx_fft = ((np.abs(Xk) * dt) ** 2) / T.代S_xx_fft入上面的等式并重新排列术语,我们到达

np.sum(np.abs(Xk) ** 2) = (T / (dt ** 2)) * (df_welch / df_fft) * np.sum(S_xx_welch)

现在,上述等式中的左侧(LHS)看起来非常接近于从DFT分量计算的信号中总能量的表达式.现在,请注意,信号中的采样点数量T / dt = N在哪里N.除以N,我们现在有一个LHS,根据定义,它等于E_fft上面计算的.因此,我们可以从韦尔奇的PSD中获得信号的总能量

# Signal energy from Welch's PSD
E_welch = (1. / dt) * (df_welch / df_fft) * np.sum(S_xx_welch)

E,E_fft并且E_welch应该都非常接近价值:)正如前一节末尾所讨论的那样,我们确实期望E_welch与E和之间存在一些细微差别E_fft,但这可归因于从Welch方法得到的值减少了随机误差的事实(即更准确).