如何计算Python中傅立叶变换的95%置信度？
How to calculate 95% confidence level of Fourier transform in Python?
计算 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 (and in particular Pankofsky et. al.)，但找不到公式的明确推导以及结果应用的确切条件。另一方面，我找到了一些其他参考资料，可以更容易地确认推导。
基于the answer to this question on dsp.stackexchange.com, if you only had white gaussian noise with unit variance, the squared-amplitude of each Fourier coefficients would have Chi-squared distribution with degree of freedom asymptotically 2 (sum of 2 Gaussians, one for each of the real and imaginary parts of the complex Fourier coefficient, when n >> 1). When the noise does not have unit variance, it follows a more general Gamma distribution (although in this case you can simply think of it as scaling the survival function). For noise with a uniform distribution in the [-0.5,0.5] range, and a sufficiently large number of samples, the distribution can also be approximated by a Gamma distribution thanks to the Central Limit Theorem.
为了说明和更好地理解这些分布，我们可以逐步研究更复杂的案例。
随机噪声的频域分布
为了与后一种均匀分布数据的情况进行比较，我们将使用具有匹配方差的高斯噪声。由于均匀分布数据的方差在 [-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,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='--')
最后，虽然我将最终绘图保留为幅度平方单位，但没有什么能阻止您取平方根并以幅度单位查看相应的阈值。
Edit ：我使用了 gamma(1,s) 分布，它是具有足够数量样本 n 的数据的渐近良好分布。对于非常小的数据大小，分布更接近 gamma(0.5*(n/(n//2+1)),s)（由于 DC 和 Nyquist 系数是纯实数，因此与所有其他系数不同，具有 1 个自由度）。
计算 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 (and in particular Pankofsky et. al.)，但找不到公式的明确推导以及结果应用的确切条件。另一方面，我找到了一些其他参考资料，可以更容易地确认推导。
基于the answer to this question on dsp.stackexchange.com, if you only had white gaussian noise with unit variance, the squared-amplitude of each Fourier coefficients would have Chi-squared distribution with degree of freedom asymptotically 2 (sum of 2 Gaussians, one for each of the real and imaginary parts of the complex Fourier coefficient, when n >> 1). When the noise does not have unit variance, it follows a more general Gamma distribution (although in this case you can simply think of it as scaling the survival function). For noise with a uniform distribution in the [-0.5,0.5] range, and a sufficiently large number of samples, the distribution can also be approximated by a Gamma distribution thanks to the Central Limit Theorem.
为了说明和更好地理解这些分布，我们可以逐步研究更复杂的案例。
随机噪声的频域分布
为了与后一种均匀分布数据的情况进行比较，我们将使用具有匹配方差的高斯噪声。由于均匀分布数据的方差在 [-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,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='--')
最后，虽然我将最终绘图保留为幅度平方单位，但没有什么能阻止您取平方根并以幅度单位查看相应的阈值。
Edit ：我使用了 gamma(1,s) 分布，它是具有足够数量样本 n 的数据的渐近良好分布。对于非常小的数据大小，分布更接近 gamma(0.5*(n/(n//2+1)),s)（由于 DC 和 Nyquist 系数是纯实数，因此与所有其他系数不同，具有 1 个自由度）。