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

如何解决如何在 Python 中计算傅立叶变换的 95% 置信度？

在 Python/Scipy 中计算时间序列的快速傅立叶变换 (FFT) 后，我试图绘制功率谱不同于红噪声或白噪声的 95% 置信水平，但尚未找到一种直接的方法。我尝试关注此线程：Power spectrum in python - significance levels 并编写了以下代码来测试具有随机噪声的正弦函数：

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 功率谱的显着性？

解决方法

背景注意事项

我没有阅读整个 references included in the answer you linked to（尤其是 Pankofsky et. al.），但找不到公式的明确推导以及结果适用的确切条件。另一方面，我发现了一些其他参考资料，可以更容易地确认推导。

基于the answer to this question on dsp.stackexchange.com，如果只有具有单位方差的高斯白噪声，则每个傅立叶系数的平方幅度将具有自由度渐近为2的卡方分布（ 2 个高斯的总和，当 n >> 1 时，复傅立叶系数的实部和虚部各一个。当噪声没有单位方差时，它遵循更一般的 Gamma distribution（尽管在这种情况下您可以简单地将其视为缩放生存函数）。对于在 [-0.5,0.5] 范围内均匀分布且样本数量足够多的噪声，由于 Central Limit Theorem，该分布也可以近似为 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's theorem得到时域方差：

# 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,10*np.log10(np.abs(yf)**2/n),color='k')  
plt.axhline(y=10*np.log10(gamma.isf(q=0.05/2,scale=var)),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,scale=noiseVar)
plt.figure(figsize=(5,color='k')  
plt.plot(xf,10*np.log10(threshold95),linestyle='--')

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

---

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