用Python实战解析随机信号:从统计特征到平稳性检验
在传感器监测、金融时间序列分析或语音处理中,我们常遇到看似毫无规律的数据波动。传统教学中复杂的数学公式往往让人望而生畏,而实际工程更需要的是能快速上手的代码实现。本文将以Jupyter Notebook为战场,用NumPy和Matplotlib这些Python利器,带您直击随机信号分析的核心实战技能。
1. 理解随机信号的DNA:统计特征
随机信号就像一群自由奔放的野马,虽然单个轨迹难以预测,但群体行为却存在规律。我们首先需要掌握描述这些规律的三大核心指标:
import numpy as np import matplotlib.pyplot as plt # 生成模拟随机信号 np.random.seed(42) stationary_signal = np.cumsum(np.random.randn(1000)) * 0.1 non_stationary = np.concatenate([ np.random.normal(0, 1, 300), np.random.normal(2, 1.5, 400), np.random.normal(-1, 0.8, 300) ])1.1 均值与方差:信号的身份证
均值反映信号的基准线,方差则体现波动强度。对于上述模拟信号:
def print_stats(signal, name): print(f"{name}信号统计:") print(f"均值: {np.mean(signal):.4f}") print(f"方差: {np.var(signal):.4f}") print(f"标准差: {np.std(signal):.4f}\n") print_stats(stationary_signal, "平稳") print_stats(non_stationary, "非平稳")注意:实际工程中建议使用pandas.DataFrame.rolling()计算滑动统计量,能更直观观察局部特征变化。
1.2 概率密度分布:信号的指纹图谱
直方图和核密度估计(KDE)比单一统计量更能揭示信号本质:
from scipy.stats import gaussian_kde plt.figure(figsize=(12,5)) plt.subplot(121) plt.hist(stationary_signal, bins=30, density=True, alpha=0.6) kde = gaussian_kde(stationary_signal) xvals = np.linspace(min(stationary_signal), max(stationary_signal), 200) plt.plot(xvals, kde(xvals), 'r-') plt.title("平稳信号分布") plt.subplot(122) plt.hist(non_stationary, bins=30, density=True, alpha=0.6) kde = gaussian_kde(non_stationary) xvals = np.linspace(min(non_stationary), max(non_stationary), 200) plt.plot(xvals, kde(xvals), 'r-') plt.title("非平稳信号分布") plt.show()2. 平稳性检验:信号稳定性的照妖镜
平稳性是许多信号处理算法的前提假设,下面介绍两种实用的检验方法。
2.1 滑动窗口检验法
通过观察统计量随时间的变化来判断平稳性:
def rolling_stat_test(signal, window_size=100): means = [np.mean(signal[i:i+window_size]) for i in range(0, len(signal)-window_size)] stds = [np.std(signal[i:i+window_size]) for i in range(0, len(signal)-window_size)] plt.figure(figsize=(12,4)) plt.subplot(121) plt.plot(means) plt.title("滑动均值变化") plt.subplot(122) plt.plot(stds) plt.title("滑动标准差变化") plt.show() rolling_stat_test(stationary_signal) rolling_stat_test(non_stationary)2.2 ADF单位根检验
Augmented Dickey-Fuller检验是时间序列分析中的经典方法:
from statsmodels.tsa.stattools import adfuller def adf_test(signal): result = adfuller(signal) print('ADF统计量: %f' % result[0]) print('p值: %f' % result[1]) print('临界值:') for key, value in result[4].items(): print('\t%s: %.3f' % (key, value)) if result[1] < 0.05: print("拒绝原假设,信号可能平稳") else: print("无法拒绝原假设,信号可能非平稳") print("平稳信号检验结果:") adf_test(stationary_signal) print("\n非平稳信号检验结果:") adf_test(non_stationary)3. 实战案例:EEG信号分析
让我们用真实的EEG数据(来自MNE库)演示完整分析流程:
# 安装MNE库: pip install mne import mne from mne.datasets import sample data_path = sample.data_path() raw_fname = data_path / 'MEG/sample/sample_audvis_filt-0-40_raw.fif' raw = mne.io.read_raw_fif(raw_fname, preload=True) # 选取EEG通道 eeg_data = raw.pick_types(eeg=True).get_data()[0, :10000]3.1 特征提取与可视化
plt.figure(figsize=(15,5)) plt.plot(eeg_data[:1000]) plt.title("EEG信号片段") plt.xlabel("采样点") plt.ylabel("幅值(μV)") plt.show() # 计算统计特征 eeg_features = { "均值": np.mean(eeg_data), "方差": np.var(eeg_data), "偏度": mne.stats.stats.skew(eeg_data), "峰度": mne.stats.stats.kurtosis(eeg_data) } print("EEG信号特征:") for k, v in eeg_features.items(): print(f"{k}: {v:.4f}")3.2 平稳性分析
# 分段平稳性检验 segment_length = 2000 for i in range(0, len(eeg_data), segment_length): segment = eeg_data[i:i+segment_length] print(f"\n段{i//segment_length}检验结果:") adf_test(segment)4. 进阶技巧:处理非平稳信号
当检测到非平稳性时,常用的处理方法包括:
| 方法 | 原理 | 适用场景 | Python实现 |
|---|---|---|---|
| 差分法 | 计算相邻数据点的差值 | 趋势性非平稳 | np.diff() |
| 对数变换 | 压缩数据尺度 | 方差变化型 | np.log1p() |
| 分段处理 | 将信号划分为准平稳段 | 突变型非平稳 | 自定义分段 |
| 小波变换 | 时频域联合分析 | 复杂非平稳 | pywt库 |
差分法示例:
diff_signal = np.diff(non_stationary) plt.figure(figsize=(12,4)) plt.subplot(121) plt.plot(non_stationary) plt.title("原始信号") plt.subplot(122) plt.plot(diff_signal) plt.title("一阶差分信号") plt.show() print("差分后平稳性检验:") adf_test(diff_signal)在金融时间序列分析中,我经常发现简单的对数差分组合就能使大多数价格序列变得平稳。但要注意过度差分会导致信息损失,实践中需要通过ACF/PACF图辅助判断最优差分阶数。