如果我的数据是 1/f 噪声，如何从功率谱密度确定？答案

【问题标题】：How to determine from the Power Spectral Density if my data is 1/f noise?如果我的数据是 1/f 噪声，如何从功率谱密度确定？
【发布时间】：2017-02-28 14:15:48
【问题描述】：

假设你有以下一系列值：

import numpy as np
data1=np.array([1,  0.28571429,  0.20408163,  0.16326531,  0.14285714,
        0.10204082,  0.08163265,  0.06122449,  0.04081633,  0.02040816])

您现在想要使用numpy.fft 绘制data1 的Spectral Density：

ps1 = np.abs(np.fft.fft(data1))**2
time_step = 1 / 30
freqs1 = np.fft.fftfreq(data1.size, time_step)
idx1 = np.argsort(freqs1)
fig,ax=plt.subplots()
ax.plot(freqs1[idx1], ps1[idx1],color='red',label='data1')
plt.grid(True, which="both")
ax.set_yscale('log')
ax.set_xscale('log')
plt.xlabel("Frequency")
plt.ylabel("Power")
plt.xlim(0,15)

我的问题：如果该图代表我的系列的信号，如果这是1/f（或任何其他）噪音？

【问题讨论】：

你怎么知道数据中有噪音？相反，您怎么知道数据中除了噪声之外还有其他内容？而且，重要的是，您对噪音的定义是什么？
我不知道是否存在噪音。我想弄清楚。这可能是带有1/f 噪音或其他东西的系列。我怎样才能弄清楚？
这取决于，但它肯定是我们可以使用的。我最初认为数据中存在有用的信号和噪声，而您想提取噪声。如果数据实际上是噪音，并且您想知道它是否真的是您认为的那种噪音，请查看我在下面发布的更长答案。
我冒昧地简单地编辑了这个问题。请确保措辞仍然符合您的兴趣并且是正确的（我不是母语人士）。哦，检查更新的答案。它似乎真的有 1/f 噪声的频谱，除了 DC 部分，无论如何这很棘手。
表示数据平均不为0。DC部分是数据的均值，相当于0Hz的谱峰。

标签： python numpy signal-processing fft frequency-analysis

【解决方案1】：

如果我看这张图片，我很想得出结论，这个系列有类似幂律行为，但要说这个，你需要证明它有 1/f 噪声。

让我们看看我们是否可以验证这个假设：

为 1/f 噪声构建噪声模型：scale / (1 + f ** alpha)。这是一个 1/f niose 的数学模型，参数为scale（幅度）和alpha（描述高频与低频的比率）
将噪声模型拟合到数据（在本例中将其拟合到功率谱密度）
检查结果。模型看起来是否很好地描述了数据？如果没有 - 尝试找到不同的 niose 模型。但要小心overfitting！

截图：

from scipy.optimize import curve_fit
import numpy as np
import matplotlib.pyplot as plt


data1 = np.array([1,  0.28571429,  0.20408163,  0.16326531,  0.14285714,
                  0.10204082,  0.08163265,  0.06122449,  0.04081633,  0.02040816])


def one_over_f(f, alpha, scale):
    return scale / f ** alpha    


time_step = 1 / 30
ps1 = np.abs(np.fft.fft(data1)) ** 2
freqs1 = np.fft.fftfreq(data1.size, time_step)

# remove the DC bin because 1/f cannot fit at 0 Hz
ps1 = ps1[freqs1 != 0]
freqs1 = freqs1[freqs1 != 0]

params, _ = curve_fit(one_over_f, np.abs(freqs1), ps1)
ps2 = one_over_f(freqs1, *params)

idx1 = np.argsort(freqs1)

fig, ax = plt.subplots()
ax.plot(freqs1[idx1], ps1[idx1], color='red', label='data1')
ax.plot(freqs1[idx1], ps2[idx1], color='blue', label='fit')
plt.grid(True, which="both")
ax.set_yscale('log')
ax.set_xscale('log')
plt.xlabel("Frequency")
plt.ylabel("Power")
plt.xlim(0,15)
plt.legend()

print('Alpha:', params[0], '\nScale:', params[1])

从视觉上看，该模型似乎非常适合 1/f 频谱。这不是证据。很难证明数据有一定的分布。但您可以自行判断所选噪声模型是否在视觉上符合您的需求。

【讨论】：