宽带高斯过程，窄带高斯过程
宽带高斯过程：时间轴上“几乎没有固定节奏”，相邻过零点间隔非常随机、非常短。
窄带高斯过程：时间轴上“有近似固定的节奏”，相邻过零点间隔虽然仍有随机性，但平均值接近一个固定周期
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, rayleigh
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# ============================================================
# 1. 参数设置
# ============================================================
fs = 1000.0 # 采样频率 (Hz)
T = 20.0 # 总时长 (s)
N = int(fs * T) # 采样点数
t = np.linspace(0, T, N, endpoint=False)
# 频率轴
freq = np.fft.rfftfreq(N, d=1/fs) # 0 ~ fs/2
omega = 2 * np.pi * freq
# ============================================================
# 2. 定义宽带和窄带谱密度形状
# ============================================================
# 宽带：白噪声，谱密度在整个频带内为常数
S_wide = np.ones_like(omega)
# 窄带：中心频率附近很窄的矩形，Delta_omega << omega0
f0 = 50.0 # 中心频率 (Hz)
omega0 = 2 * np.pi * f0
Delta_f = 2.0 # 带宽 (Hz)
Delta_omega = 2 * np.pi * Delta_f
S_narrow = np.zeros_like(omega)
mask = (omega > omega0 - Delta_omega/2) & (omega < omega0 + Delta_omega/2)
S_narrow[mask] = 1.0
# ============================================================
# 3. 用频域法生成高斯随机过程
# ============================================================
def generate_gaussian_process(S_omega, N, fs):
"""
根据给定的单边谱密度 S_omega 生成平稳高斯过程。
方法：在频域赋予随机相位，再做逆傅里叶变换。
"""
# 双边谱（因为 rfft 只取正频率，需要构造完整的双边谱）
S_full = np.concatenate([S_omega, S_omega[-2:0:-1]])
# 随机相位
phase = np.random.uniform(0, 2*np.pi, len(S_full))
# 频域幅值：sqrt(S) * exp(i*phase)
X = np.sqrt(S_full) * np.exp(1j * phase)
# 逆傅里叶变换得到时域信号
x = np.fft.ifft(X).real
x = x * np.sqrt(fs * N) # 归一化，使方差与谱密度积分一致
return x
np.random.seed(42)
y_wide = generate_gaussian_process(S_wide, N, fs)
y_narrow = generate_gaussian_process(S_narrow, N, fs)
# 去除均值，确保零均值
y_wide -= np.mean(y_wide)
y_narrow -= np.mean(y_narrow)
# ============================================================
# 4. 计算统计量
# ============================================================
def compute_stats(y, fs):
"""计算均值、方差、导数方差、穿零频率、峰值"""
m = np.mean(y)
sigma2 = np.var(y)
# 数值求导
dy = np.gradient(y, 1/fs)
m_dy = np.mean(dy)
sigma2_dy = np.var(dy)
# 正穿越零的次数
zero_crossings = np.where((y[:-1] < 0) & (y[1:] >= 0))[0]
nu0 = len(zero_crossings) / (len(y) / fs)
# 峰值（局部极大值）
peaks = []
for i in range(1, len(y)-1):
if y[i] > y[i-1] and y[i] > y[i+1]:
peaks.append(y[i])
peaks = np.array(peaks)
return {
'm': m,
'sigma2': sigma2,
'sigma': np.sqrt(sigma2),
'm_dy': m_dy,
'sigma2_dy': sigma2_dy,
'sigma_dy': np.sqrt(sigma2_dy),
'nu0': nu0,
'peaks': peaks
}
stats_wide = compute_stats(y_wide, fs)
stats_narrow = compute_stats(y_narrow, fs)
print("=" * 60)
print("宽带高斯过程统计量")
print("=" * 60)
for k, v in stats_wide.items():
if k != 'peaks':
print(f"{k:12s} = {v:.6f}")
print("\n" + "=" * 60)
print("窄带高斯过程统计量")
print("=" * 60)
for k, v in stats_narrow.items():
if k != 'peaks':
print(f"{k:12s} = {v:.6f}")
# ============================================================
# 5. 可视化对比
# ============================================================
fig, axes = plt.subplots(3, 2, figsize=(14, 12))
# 5.1 时间波形（前2秒）
t_show = t < 2.0
axes[0, 0].plot(t[t_show], y_wide[t_show], 'b-', linewidth=0.8)
axes[0, 0].set_title('宽带高斯过程：时间波形')
axes[0, 0].set_xlabel('时间 (s)')
axes[0, 0].set_ylabel('y(t)')
axes[0, 0].grid(True, alpha=0.3)
axes[0, 1].plot(t[t_show], y_narrow[t_show], 'r-', linewidth=0.8)
axes[0, 1].set_title('窄带高斯过程：时间波形')
axes[0, 1].set_xlabel('时间 (s)')
axes[0, 1].set_ylabel('y(t)')
axes[0, 1].grid(True, alpha=0.3)
# 5.2 谱密度
axes[1, 0].plot(freq, S_wide, 'b-', linewidth=1.5)
axes[1, 0].set_title('宽带谱密度 S(ω)')
axes[1, 0].set_xlabel('频率 (Hz)')
axes[1, 0].set_ylabel('S(ω)')
axes[1, 0].set_xlim(0, 200)
axes[1, 0].grid(True, alpha=0.3)
axes[1, 1].plot(freq, S_narrow, 'r-', linewidth=1.5)
axes[1, 1].set_title('窄带谱密度 S(ω)')
axes[1, 1].set_xlabel('频率 (Hz)')
axes[1, 1].set_ylabel('S(ω)')
axes[1, 1].set_xlim(0, 200)
axes[1, 1].grid(True, alpha=0.3)
# 5.3 一维概率密度
y_range = np.linspace(-4, 4, 200)
axes[2, 0].hist(y_wide, bins=50, density=True, alpha=0.6, color='b', label='宽带直方图')
axes[2, 0].plot(y_range, norm.pdf(y_range, 0, stats_wide['sigma']), 'b-', linewidth=2, label='高斯拟合')
axes[2, 0].set_title('宽带过程：一维概率密度 p(y)')
axes[2, 0].set_xlabel('y')
axes[2, 0].set_ylabel('p(y)')
axes[2, 0].legend()
axes[2, 0].grid(True, alpha=0.3)
axes[2, 1].hist(y_narrow, bins=50, density=True, alpha=0.6, color='r', label='窄带直方图')
axes[2, 1].plot(y_range, norm.pdf(y_range, 0, stats_narrow['sigma']), 'r-', linewidth=2, label='高斯拟合')
axes[2, 1].set_title('窄带过程：一维概率密度 p(y)')
axes[2, 1].set_xlabel('y')
axes[2, 1].set_ylabel('p(y)')
axes[2, 1].legend()
axes[2, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# ============================================================
# 6. 峰值分布对比：宽带 vs 窄带 vs 瑞利
# ============================================================
fig2, axes2 = plt.subplots(1, 2, figsize=(14, 5))
# 宽带峰值分布
peaks_wide = stats_wide['peaks']
peaks_wide = peaks_wide[peaks_wide > 0] # 只看正峰
axes2[0].hist(peaks_wide, bins=40, density=True, alpha=0.6, color='b', label='宽带峰值直方图')
a_range = np.linspace(0, 4*stats_wide['sigma'], 200)
axes2[0].plot(a_range, rayleigh.pdf(a_range, scale=stats_wide['sigma']), 'b--', linewidth=2, label='瑞利分布拟合')
axes2[0].set_title('宽带过程：峰值分布（不服从瑞利分布）')
axes2[0].set_xlabel('峰值 a')
axes2[0].set_ylabel('p_p(a)')
axes2[0].legend()
axes2[0].grid(True, alpha=0.3)
# 窄带峰值分布
peaks_narrow = stats_narrow['peaks']
peaks_narrow = peaks_narrow[peaks_narrow > 0]
axes2[1].hist(peaks_narrow, bins=40, density=True, alpha=0.6, color='r', label='窄带峰值直方图')
a_range = np.linspace(0, 4*stats_narrow['sigma'], 200)
axes2[1].plot(a_range, rayleigh.pdf(a_range, scale=stats_narrow['sigma']), 'r--', linewidth=2, label='瑞利分布拟合')
axes2[1].set_title('窄带过程：峰值分布（服从瑞利分布）')
axes2[1].set_xlabel('峰值 a')
axes2[1].set_ylabel('p_p(a)')
axes2[1].legend()
axes2[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
浙公网安备 33010602011771号