非线性奇异谱分解算法：精细化处理时间序列数据，提取CSV文件信号特征，生成希尔伯特谱分析报告

SSD–fft–hht，奇异谱分解算法，是对原始小波分解的一种改进，对小波分解中的高频部分进行二次分解，提高分辨率。 一种非线性时间序列分解方法，可用于处理各种复杂数据，包括金融，气候，生物医学等等。 主要功能: 1. 读取 CSV 文件中的数据并进行预处理,包括设置采样频率、时间间隔等参数。 2. 绘制原始信号的时域图像。 3. 使用 SSD (Sparse Signal Decomposition) 算法对信号进行分解,得到多个 SSD 分量。 4. 绘制 SSD 分量的时域图像、包络线、包络谱等。 5. 计算并绘制 SSD 分量的中心频率、峰度值、能量熵和包络熵。 6. 绘制 SSD 分量特征值的柱状图和峭度值柱状图。 7. 绘制 SSD 分量的希尔伯特谱。 8. 使用前 3 个 SSD 分量重构信号,并比较重构信号与原始信号。 9. 绘制 SSD 分量的功率谱密度图。 10. 计算重构后的信号与原始信号的相关系数。 总的来说,这个代码主要用于对时域信号进行 SSD 分解,并分析各个 SSD 分量的特征,最终重构信号并与原始信号进行比较。 这种方法在信号处理和分析领域有广泛的应用。

最近在信号处理领域折腾了个有意思的工具——基于SSD-FFT-HHT的奇异谱分解。这玩意儿算是小波分解的升级版，尤其擅长对付高频信号分辨率不足的老大难问题。咱们今天直接上代码实战，手把手拆解它的十八般武艺。

先整点测试数据最实在。假设咱们手头有个混了3个频率的仿真信号：

import numpy as np import matplotlib.pyplot as plt fs = 1000 # 采样率 t = np.arange(0, 1, 1/fs) # 1秒时长 f1, f2, f3 = 20, 100, 300 # 三个特征频率 signal = 1.5*np.sin(2*np.pi*f1*t) + 0.8*np.cos(2*np.pi*f2*t) + 0.3*np.sin(2*np.pi*f3*t)

原始信号长这样：

plt.figure(figsize=(12,4)) plt.plot(t, signal) plt.title('原始信号时域图') plt.xlabel('时间（秒）') plt.ylabel('幅值') plt.tight_layout() plt.show()

!原始信号时域图示例

接下来是重头戏——SSD分解。这里用到了自定义的分解类（具体算法实现涉及专利，咱们用等效代码示意）：

from scipy.signal import hilbert class SSDExtractor: def __init__(self, signal, fs): self.signal = signal self.fs = fs def decompose(self, n_components=5): # 返回分量的字典列表 components = [] # 伪代码，实际需实现分解算法 for i in range(n_components): comp = {'data': np.random.randn(len(self.signal)) * 0.1, # 示例数据 'envelope': np.abs(hilbert(comp['data'])), 'freq': np.random.choice([f1, f2, f3])} components.append(comp) return components extractor = SSDExtractor(signal, fs) components = extractor.decompose(n_components=5)

分解后的分量分析才是精华。咱们搞个组合图展示前三分量：

fig, axs = plt.subplots(3, 2, figsize=(15, 10)) for i in range(3): # 时域波形+包络线 axs[i,0].plot(t, components[i]['data'], label='分量') axs[i,0].plot(t, components[i]['envelope'], 'r--', label='包络线') axs[i,0].set_title(f'分量{i+1}时域特征') # 包络谱 f = np.fft.rfftfreq(len(t), 1/fs) spectrum = np.abs(np.fft.rfft(components[i]['envelope'])) axs[i,1].plot(f, spectrum) axs[i,1].set_title(f'分量{i+1}包络谱') plt.tight_layout() plt.show()

!分量分析示例图

SSD–fft–hht，奇异谱分解算法，是对原始小波分解的一种改进，对小波分解中的高频部分进行二次分解，提高分辨率。 一种非线性时间序列分解方法，可用于处理各种复杂数据，包括金融，气候，生物医学等等。 主要功能: 1. 读取 CSV 文件中的数据并进行预处理,包括设置采样频率、时间间隔等参数。 2. 绘制原始信号的时域图像。 3. 使用 SSD (Sparse Signal Decomposition) 算法对信号进行分解,得到多个 SSD 分量。 4. 绘制 SSD 分量的时域图像、包络线、包络谱等。 5. 计算并绘制 SSD 分量的中心频率、峰度值、能量熵和包络熵。 6. 绘制 SSD 分量特征值的柱状图和峭度值柱状图。 7. 绘制 SSD 分量的希尔伯特谱。 8. 使用前 3 个 SSD 分量重构信号,并比较重构信号与原始信号。 9. 绘制 SSD 分量的功率谱密度图。 10. 计算重构后的信号与原始信号的相关系数。 总的来说,这个代码主要用于对时域信号进行 SSD 分解,并分析各个 SSD 分量的特征,最终重构信号并与原始信号进行比较。 这种方法在信号处理和分析领域有广泛的应用。

特征值计算这块，能量熵和峭度值最实用：

def calc_features(component): data = component['data'] # 能量熵 energy = np.sum(data**2) p = data**2 / energy entropy = -np.sum(p * np.log(p)) # 峭度 kurtosis = np.mean((data - np.mean(data))**4) / (np.std(data)**4) return {'能量熵': entropy, '峭度': kurtosis} for i, comp in enumerate(components): features = calc_features(comp) print(f'分量{i+1}：能量熵={features["能量熵"]:.2f}，峭度={features["峭度"]:.2f}')

信号重构与对比是检验分解效果的关键：

reconstructed = sum([comp['data'] for comp in components[:3]]) plt.figure(figsize=(12,4)) plt.plot(t, signal, alpha=0.5, label='原始信号') plt.plot(t, reconstructed, '--', label='重构信号') plt.legend() plt.title('重构信号对比') plt.show() # 计算相关系数 corr_coef = np.corrcoef(signal, reconstructed)[0,1] print(f'相关系数：{corr_coef:.3f}')

实际应用中遇到过个有意思的案例：某轴承振动信号里混着转速频率和故障冲击，传统小波分解总把冲击特征分散到多个分量里。换成SSD分解后，前两个分量直接锁定了轴频和故障冲击，包络谱里的特征频率清晰得像教科书似的。

不过要注意，这算法对采样率敏感。上次处理EEG信号时，采样率不足导致高频分量出现混叠，后来插值到2000Hz才解决。建议正式分析前先做抗混叠滤波，别像我一样踩坑。

最后放个功率谱密度图收尾，方便观察频率分布：

from scipy.signal import welch fig, axs = plt.subplots(2, 1, figsize=(12,6)) for i in range(2): f, Pxx = welch(components[i]['data'], fs) axs[i].semilogy(f, Pxx) axs[i].set_title(f'分量{i+1}功率谱密度') plt.tight_layout() plt.show()

这套方法在金融时间序列预测中也挺有用——把价格序列分解成趋势、周期、噪声分量后，分别用不同策略处理。不过要注意市场数据的非平稳特性，建议配合自适应窗长使用。