python - 如何在 Python 中计算傅里叶变换的 95% 置信度？

Question

在 Python/Scipy 中计算时间序列的快速傅里叶变换 (FFT) 之后，我试图绘制功率谱不同于红噪声或白噪声的 95% 置信水平，但还没有找到直接的方法这样做。我尝试关注这个线程：Python中的功率谱-显着性水平 并编写了以下代码来测试具有随机噪声的正弦函数：

import numpy as np
from scipy.stats import chi2
from scipy.fft import rfft, rfftfreq

x=np.linspace(0,10,500)

data = np.sin(20*np.pi*x)+np.random.rand(500) - 0.5

yf = rfft(data)
xf = rfftfreq(len(data), 1) 
n=len(data)

var=np.var(data)

### degrees of freedom
M=n/2
phi=(2*(n-1)-M/2.)/M       
###values of chi-squared
chi_val_99 = chi2.isf(q=0.01/2, df=phi) #/2 for two-sided test
chi_val_95 = chi2.isf(q=0.05/2, df=phi)



### normalization of power spectrum with 1/n
plt.figure(figsize=(5,5))
plt.plot(xf,np.abs(yf)/n, color='k')  
plt.axhline(y=(var/n)*(chi_val_95/phi),color='r',linestyle='--') 

但结果线位于所有功率谱下方，如图 1 所示。我做错了什么？有没有另一种方法来获得 FFT 功率谱的重要性？

标签： pythonfftscipy.stats

Answer 1

背景考虑

我没有阅读您链接到的答案中包含的全部参考资料（特别是Pankofsky 等人），但找不到公式的明确推导以及结果应用的确切条件。另一方面，我发现了一些其他参考资料，可以更容易地确认推导。

根据dsp.stackexchange.com 上对这个问题的回答，如果您只有具有单位方差的高斯白噪声，则每个傅立叶系数的平方幅度将具有自由度渐近为 2 的卡方分布（2 个高斯的总和） ，对于复傅立叶系数的实部和虚部各有一个，当n >> 1)。当噪声没有单位方差时，它遵循更一般的Gamma 分布（尽管在这种情况下，您可以简单地将其视为缩放生存函数）。对于范围内分布均匀且样本数量足够多的噪声，由于中心极限定理[-0.5,0.5]，分布也可以近似为 Gamma 分布.

为了说明和更好地理解这些分布，我们可以逐步了解更复杂的案例。

随机噪声的频域分布

为了与后面的均匀分布数据的情况进行比较，我们将使用具有匹配方差的高斯噪声。由于均匀分布数据的方差在 范围内[-0.5,0.5]，1/12这给了我们以下数据：

data = np.sqrt(1.0/12)*np.random.randn(500)

现在让我们检查一下功率谱的统计数据。如前所述，每个频率系数的平方幅度是具有近似 Gamma 分布的随机变量。形状参数是可以用于单位方差情况的卡方分布自由度的一半（在这种情况下为 1），并且比例参数对应于时域缩放的平方（从线性度来看，变量yf缩放为data，使得np.abs(yf)**2缩放为 的平方data）。我们可以通过绘制data概率密度函数的直方图来验证这一点：

yf = rfft(data)
spectrum = np.abs(yf)**2/len(data)
plt.figure(figsize=(5,5))
plt.hist(spectrum, bins=100, density=True, label='data')
z = np.linspace(0, np.max(spectrum), 100)
plt.plot(z, gamma.pdf(z, 1, scale=1.0/12), 'k', label='$\Gamma(1,{:.3f})$'.format(1.0/12))

如您所见，这些值非常一致：

回到频谱图：

# degrees of freedom
phi = 2
###values of chi-squared
chi_val_95 = chi2.isf(q=0.05/2, df=phi) #/2 for two-sided test

### normalization of power spectrum with 1/n
plt.figure(figsize=(5,5))
plt.plot(xf,np.abs(yf)**2/n, color='k')
# the following two lines should overlap
plt.axhline(y=var*(chi_val_95/phi),color='r',linestyle='--')
plt.axhline(y=gamma.isf(q=0.05/2, a=1, scale=var),color='b')

只需将 更改data为使用[-0.5,0.5]范围内的均匀分布（与data = np.random.rand(500) - 0.5）即可得到几乎相同的图，而置信水平保持不变。

带噪声信号的频域分布

要获得对应于 95% 置信区间的单个阈值，如果您可以将其与data包含正弦分量和噪声的噪声部分分开（或以其他方式表示为 95% 置信区间的零假设的置信区间，则噪声部分将下降data）白噪声），您将需要噪声的方差。在尝试估计此方差时，您可能会很快意识到正弦曲线在整体中贡献了不可忽略的部分data的方差。为了消除这种影响，我们可以利用正弦信号在频域中更容易分离的事实。因此，我们可以简单地丢弃频谱的 x% 最大值，假设这些最大值主要是由频域中正弦分量的尖峰贡献的。请注意，对于异常值，以下 95% 的选择有些随意：

# remove outliers
threshold = np.percentile(np.abs(yf)**2, 95) 
filtered = [x for x in np.abs(yf)**2 if x <= threshold]

然后我们可以使用Parseval 定理得到时域方差：

# estimate variance
# In time-domain variance ~ np.sum(data**2)/len(data))
# In frequency-domain, using Parseval's theorem we get np.sum(data**2)/len(data) = np.mean(np.abs(spectrum)**2)/len(data)
var = np.mean(filtered)/len(data)

请注意，由于频谱中值的动态范围，您可能更喜欢以对数刻度可视化结果：

plt.figure(figsize=(5,5))
plt.plot(xf,10*np.log10(np.abs(yf)**2/n), color='k')  
plt.axhline(y=10*np.log10(gamma.isf(q=0.05/2, a=1, scale=var)),color='r',linestyle='--')

另一方面，如果您试图获得与频率相关的 95% 置信区间，那么您需要考虑每个频率的正弦分量的贡献。为了简单起见，我们在这里假设正弦分量的幅度和噪声的方差是已知的（否则我们首先需要估计这些）。在这种情况下，分布会因正弦分量的贡献而发生变化：

signal = np.sin(20*np.pi*x)
data = signal + np.random.rand(500) - 0.5
Sf   = rfft(signal)  # Assuming perfect knowledge of the sinusoidal component
yf   = rfft(data)

noiseVar = 1.0/12    # Assuming perfect knowledge of the noise variance
threshold95 = np.abs(Sf)**2/n + gamma.isf(q=0.05/2, a=1, scale=noiseVar)
plt.figure(figsize=(5,5))
plt.plot(xf, 10*np.log10(np.abs(yf)**2/n), color='k')  
plt.plot(xf, 10*np.log10(threshold95), color='r',linestyle='--')

最后，虽然我以平方幅度单位保留最终图，但没有什么能阻止您取平方根并以幅度单位查看相应的阈值。

---

编辑：我使用了一个gamma(1,s)分布，它对于具有足够数量的样本的数据来说是一个渐近良好的分布n。对于非常小的数据大小，分布更接近地匹配 a gamma(0.5*(n/(n//2+1)),s)（由于 DC 和 Nyquist 系数是纯实数，因此与所有其他系数不同，具有 1 个自由度）。