如何使用Python自动确定时间序列的PSD/FFT中是否没有季节性？

Question

我有大约1000个不同的时间序列，对于每一个序列，我想自动确定是否存在季节性。

假设存在季节性，可以通过FFT或PSD确定周期性。

但是如何基于FFT或PSD自动判断信号中不存在季节性或周期性呢？

def psd_time_series(y):
   yAC = np.correlate(Y-np.mean(Y), Y-np.mean(Y), mode='full')
   yAC = yAC/np.max(yAC) # not necessary, but scales large values
   fft_yAC= np.fft.fft(yAC)
   freqs = np.arange(0,len(fft_yAC))/len(fft_yAC)
   psd = 10*np.log10(np.abs(fft_yAC)/max(np.abs(fft_yAC))
   return psd,freqs

def determine_if_seasonal(psd):
    ### part I need help with

def detect_seasonality(y):

   psd,freqs = psd_time_series(y)
  
   seasonality = ... #### do some check of PSD to determine if seasonal

   if seasonality:
       periodicity = round(1/freqs[psd.argsort()[::-1]][0])
   else:
       periodicity = None
   return periodicity

自动确定单个尖峰或高斯噪声是否基于时间序列的FFT或PSD具有季节性的方法是什么？有没有关于PSD幅度阈值的经验法则？峰值的突出程度？峰值高度？

例如，单个尖峰的PSD图可能如下所示：

一次单脉冲的快速傅里叶变换（FFT）

高斯噪声的PSD可能看起来像这样。

高斯噪声的FFT

实际信号的PSD可能如下所示，具有周期性。

同一信号的FFT

感谢任何意见或洞见。

- kspr

---

我会投票关闭这个问题，但是由于它有赏金，我不能这样做。对我来说，它似乎不像一个编程问题，而更像一个关于信号处理的问题。想象一下，如果你要解决潜在的概念信号处理问题——“如何使用FFT发现这种类型信号的(a)周期性？”——还会有什么问题吗？关于Python的部分似乎是不相关的。也许信号处理堆栈交换是一个很好的选择？ - gspr

@gspr也许你是对的，我想我不能取消悬赏并关闭它？我以前从未使用过悬赏。 - kspr

不好意思，我不知道。 - gspr

我认为赏金是不可退款的。如果我有时间，在本周内可能会尝试给你一个答案。现在我可以告诉你的是，只有当你有几个周期时，PSD峰值才会显著。你是否可以对数据进行多个周期的采样？ - Bob

Answer 1

这个答案附带一个免责声明，我已经超过15年没有进行这种类型的时间序列分析了，因此非常生疏。我非常欢迎任何纠正，因为我是根据记忆工作的，可能会忘记一些事情。但是，希望它能给您一些想法，如果您的数据允许，您可以完善这种方法。（如果您这样做，请回复或编辑答案！）

正如我在上面的评论中提到的那样，如果您可以将数据集分成多个段，您就可以开始在数据集中不同频率上获得一些能量置信区间。当然，这意味着您需要长时间的数据集，这可能取决于您正在做什么。但是，像任何统计技术一样，如果您没有很多数据，很难对您发现的内容感到自信。我将描述如何确定是否存在信号，您可以倒转它以找到哪里没有信号。

我会做的第一件事是将数据集分成几个部分。你可以分成多少部分取决于你的数据集相对于你试图找到（或没有找到）的信号周期有多长。你可以做的一件事是重叠它们以获得更多的片段，但为了节省能量，你需要使用适当的窗口进行滤波。在下面的代码中，我假设50%的重叠并应用 Hanning 窗口。我忘记了不同窗口和重叠的细节，所以我很抱歉无法提供更多关于为什么选择这个重叠和这个窗口的信息。

希望从你的数据集中，你可以得出一个分布或期望，了解你的背景状态是什么。也许它只是一个白噪声地板或者更复杂的东西，是红移还是蓝移。功率谱背景水平应该给你关于噪声水平的信息，因为功率谱能量与噪声相关，公式为 2 * noise**2 / (fs * L)，其中 fs 是采样频率，L 是窗口长度。如果我没记错的话，这就是全部内容，就像我说的，我非常生疏。

首先，让我们来看看噪声和噪声与信号叠加后的功率谱。

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

def detrend(y):
    # This function removes the linear trend from the time series
    # by fitting y = mx + b to x and then subtracting it
    x = np.arange(len(y))
    A = np.vstack([x, np.ones(len(x))]).T
    m, b = np.linalg.lstsq(A, y, rcond=None)[0]
    return(y - m * x - b)

def subset(f, L, overlap, fs):
    # Break the function f into sections of length L.
    # Normally you do some overlapping of sections,
    # and you want to use an appropriate window
    # to conserve the energy in the time series signal.
    # Here I use a Hanning window.
    N = len(f)
    x = []
    for i in range(0, int(N-L+1), int(L * overlap)):
        window = np.hanning(L)
        y = detrend(f[i:i+L]) * window
        fft = np.fft.fft(y)
        fft = fft[:len(fft)//2+1]
        x.append(np.real(fft * np.conj(fft)) / (fs * L))
    return(x)

def noise(N):
    # Create a white noise field
    return(np.random.random(N) - 0.5)

N = 10**4 # length of time series
dt = 1    # sampling period
fs = 2 * np.pi / dt # sampling frequency
nyquist = fs / 2 # nyquist frequency
t = np.linspace(dt, N * dt, N) # times of data collection


L = N // 100 # window size used for subsampling
overlap = 0.5 # overlap in windows

# frequencies that will come out of the fft
# you can also use fp.fft.fftfreq and adjust
freqs = np.linspace(0, nyquist, L//2+1) 

# Create background white noise
noise_amplitude = 3
energy_scale = 2 * noise_amplitude**2 / (fs * L)
WhiteNoise = np.array(subset(noise_amplitude * noise(N), L, overlap, fs) ) / energy_scale

# Create another timeseries with a signal and noise
omega = 2 * np.pi / 10 # frequency in dataset
signal = (noise_amplitude / 4) * np.cos(t * omega)  # signal in dataset
data = signal + noise_amplitude * noise(N) # dataset
DataSet    = np.array(subset(data,  L, overlap, fs)) / energy_scale

# Percentiles of White Noise Power Spectra, for plotting
Wprct = np.percentile(WhiteNoise, [2.5, 25, 50, 75, 97.5], axis = 0)

# Percentiles of Data Power Spectra, for plotting
Dprct = np.percentile(DataSet, [2.5, 25, 50, 75, 97.5], axis = 0)

fig = plt.figure()
# Plot the spectra from all the subsets of the noise
for w in WhiteNoise:
    plt.semilogy(freqs[1:], w[1:], 'k', alpha = 0.05)

# Plot the median and confidence intervales
lines, = plt.semilogy(freqs[1:], Wprct[2,1:], label = 'Median of noise floor')
plt.fill_between(freqs[1:], Wprct[1,1:], Wprct[3,1:], alpha = 0.8, color = lines.get_color(), label = '50% CI of noise floor')
plt.fill_between(freqs[1:], Wprct[0,1:], Wprct[1,1:], alpha = 0.4, color = lines.get_color(), label = '95% CI of noise floor')
plt.fill_between(freqs[1:], Wprct[3,1:], Wprct[4,1:], alpha = 0.4, color = lines.get_color())

lines, = plt.plot(freqs[1:], Dprct[2,1:], label = 'Median signal-to-noise estimate')

plt.plot([omega, omega], [10**-3, 100], 'k:', label = 'Signal in dataset')
plt.plot([freqs[1], freqs[-1]], [1, 1], 'r:', label = 'Noise floor')
plt.legend()
 
plt.xlabel('Angular frequency')
plt.ylabel('Energy above noise floor [db]')

很明显，信号穿透了噪声，但我认为我们可以想出比“目测”更好的方法。在下面的代码中，我只是查看我的“数据集”每个频率上的点数，这些点数高于噪声分布并相应地进行缩放。从噪声中提取值会有50%的机会使您的值大于噪声，因此您希望您的信号远高于50%。（另一方面，如果该值小于50％，则表明信号低于噪声水平，这将是有趣的。）还有更高级的方法，我在这里没有涉及，但我认为会得到相同的结果。

plt.figure()
for i in range(1, len(freqs)):
    x = 0
    for j in range(len(DataSet[:,i])):
        x += np.sum(DataSet[j,i] > WhiteNoise[:,i])
    x /= len(DataSet[:,i])**2
    plt.plot(freqs[i], x, 'k.')
plt.plot([omega, omega], [0.4, 1], 'k:', label = 'Signal in dataset')
plt.xlabel('Angular frequency')
plt.ylabel('Probability signal is greater than noise')
plt.legend()

我认为有人可以使用这些方法来说明不同频率下的振荡出现和不存在的情况，这应该能回答你对数据集季节性的问题。

- ramzeek

---

同意@gspr的看法，这更像是个信号问题而不是编码问题，但是根据上面的描述，您的赏金可能是不可退款的。因此，我想给您一些建议。请在信号处理堆栈交换上进行双重检查，确保我的建议不会误导您。 - ramzeek

看第二张图时，另一个想法是，如果您查看噪声分布（大约在0.5左右），我认为这应该是正常的，您可以开始估计接近信号频率的值是噪声的可能性。 - ramzeek

Answer 2

这是一个信号处理问题，但您可能想要计算功率谱密度，可以使用 scipy.signal 轻松完成，然后设置合理的总功率阈值来分配周期性。

有一个类似问题（不是我的）的更长答案 here。

- eschibli

---

是的，补充eschibli的信息和相关帖子，PSD应该显示数据集中每个频率上有多少能量。理想情况下，您的数据集应足够长，以便使用类似于汉宁窗口的方法将其分成具有一定重叠的部分。这样，您就可以集合您的时间序列并在谱能量上获得误差条。 - ramzeek

Answer 3

你可以使用自相关函数或偏自相关函数。对于每个滞后期（周期），您将获得一个系数，该系数将测量信号与延迟版本的相似性。确保您使用大于某些周期（>5）的滞后期，并且如果所有系数都小于阈值（比如说0.5，但我没有找到任何理论方法来支持这个决定），那么您的信号就没有季节性。

- Guinther Kovalski

---

如果你愿意的话，我可以添加代码来自动完成这个任务，但前提是管理员同意不删除/关闭这个问题。 - Guinther Kovalski