使用Periodogram查找开关系统时间序列周期k的问题排查

我有一个开关系统的时间序列，序列的重复模式长度为k，模式包含k-1个连续0和1个1。尝试用Periodogram找出模式长度k，但代码运行结果有误：针对长度为7的模式，代码返回了3，而非正确的7。

原始代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import periodogram

def find_repetitive_pattern_length(time_series):
    f, Pxx = periodogram(time_series)
    max_power_index = np.argmax(Pxx)
    dominant_frequency = f[max_power_index]
    period_length = int(1 / dominant_frequency)

    
    plt.figure(figsize=(10, 6))
    plt.plot(f, Pxx, label='Periodogram')
    plt.scatter(f[max_power_index], Pxx[max_power_index], color='red', label=f'Dominant Frequency: {dominant_frequency:.2f} Hz')
    plt.title('Periodogram with Dominant Frequency')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Power/Frequency Density')
    plt.legend()
    plt.show()

    return period_length


time_series4 = np.array([0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1])
time_series7 = np.array([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1])
period_length = find_repetitive_pattern_length(time_series4)
print(f"The length of the first repetitive pattern is: {period_length}")
period_length = find_repetitive_pattern_length(time_series7)
print(f"The length of the second repetitive pattern is: {period_length}")

运行结果 #

The length of the first repetitive pattern is: 4 (correct)
The length of the second repetitive pattern is: 3 (incorrect)

错误原因 #

1. 未排除直流分量：序列中0的数量远多于1，直流分量（频率0）的功率通常很高，若误取该频率会导致周期计算错误；即使未取直流，脉冲串类信号的高次谐波功率可能与基频相当，np.argmax可能选中谐波对应的频率点。
2. 频率离散性与频谱泄漏：FFT的频率点是离散的（f = n/N，N为序列长度），当基频不在这些离散点上时，能量会泄漏到相邻频率点，导致最大功率点并非真实基频。
3. 仅依赖单个频率点：只取最大功率的单个频率点，未考虑信号的谐波特性——脉冲串信号的频谱包含基频（1/k）及其整数倍谐波，谐波对应的周期是k/m（m为整数），会得到更小的错误周期。

解决办法 #

优化思路 #

• 排除直流分量，从非零频率中寻找候选频率；
• 选取多个高功率频率点，计算对应的周期，取其中最大的周期（基频对应最长周期）；
• 结合自相关函数验证，自相关函数在延迟k、2k等处会出现明显峰值，可辅助确认模式长度。

修改后的代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import periodogram

def find_repetitive_pattern_length(time_series):
    fs = 1  # 采样频率，每个样本间隔1单位时间
    f, Pxx = periodogram(time_series, fs=fs, window='boxcar', scaling='density')
    
    # 排除直流分量（频率0），避免浮点误差用1e-6做阈值
    non_dc_mask = f > 1e-6
    f_non_dc = f[non_dc_mask]
    Pxx_non_dc = Pxx[non_dc_mask]
    
    if len(f_non_dc) == 0:
        return 0  # 全0序列，无模式
    
    # 取前5个功率最高的频率点，覆盖基频及主要谐波
    top_n = 5
    top_indices = np.argsort(Pxx_non_dc)[-top_n:][::-1]
    top_frequencies = f_non_dc[top_indices]
    # 用round取整，避免int截断导致的误差
    top_periods = np.round(1 / top_frequencies).astype(int)
    
    # 基频对应最长周期，取最大的候选周期
    period_length = max(top_periods)
    
    # 绘制周期图，标记所有高功率频率点并标注周期
    plt.figure(figsize=(10, 6))
    plt.plot(f, Pxx, label='Periodogram')
    for freq, power in zip(top_frequencies, Pxx_non_dc[top_indices]):
        plt.scatter(freq, power, color='red', s=50, zorder=5)
        plt.text(freq + 0.01, power, f'T={round(1/freq)}', fontsize=10)
    plt.title('Periodogram with Top High-Power Frequencies')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Power/Frequency Density')
    plt.legend()
    plt.grid(True)
    plt.show()
    
    # 自相关函数验证模式长度
    autocorr = np.correlate(time_series, time_series, mode='full')
    autocorr = autocorr[len(autocorr)//2:]  # 只保留正延迟部分
    # 排除延迟0的自相关峰值（自身匹配）
    autocorr_peaks = np.where(autocorr[1:] > 0.5 * np.max(autocorr[1:]))[0] + 1
    if len(autocorr_peaks) > 0:
        print(f"自相关验证的模式长度：{autocorr_peaks[0]}")
    
    return period_length


time_series4 = np.array([0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1])
time_series7 = np.array([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1])

period_length = find_repetitive_pattern_length(time_series4)
print(f"第一个序列的模式长度：{period_length}")
period_length = find_repetitive_pattern_length(time_series7)
print(f"第二个序列的模式长度：{period_length}")

代码说明 #

1. 排除直流分量：通过阈值过滤掉接近0的频率，避免误选直流分量；
2. 多频率点候选：选取前5个高功率频率点，覆盖基频及谐波，通过取最大周期得到真实模式长度；
3. 自相关验证：利用自相关函数的峰值特性，进一步确认模式长度，提高结果可靠性；
4. 绘图优化：标记所有高功率频率点并标注对应周期，便于直观观察频谱特性。

内容的提问来源于stack exchange，提问作者Lorenzo Cutrupi