窄带过程的可视化

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import rayleigh
from scipy.signal import hilbert

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# ============================================================
# 1. 参数设置
# ============================================================
fs = 2000.0
T = 20.0
N = int(fs * T)
t = np.linspace(0, T, N, endpoint=False)

f0 = 100.0
omega0 = 2 * np.pi * f0
freq = np.fft.rfftfreq(N, d=1/fs)
omega = 2 * np.pi * freq

# ============================================================
# 2. 生成三种不同带宽的窄带过程
# ============================================================
def make_spectrum(Delta_f):
    S = np.zeros_like(omega)
    mask = (omega > omega0 - np.pi*Delta_f) & (omega < omega0 + np.pi*Delta_f)
    S[mask] = 1.0
    return S

def generate_process(S_omega, N, fs):
    S_full = np.concatenate([S_omega, S_omega[-2:0:-1]])
    phase = np.random.uniform(0, 2*np.pi, len(S_full))
    X = np.sqrt(S_full) * np.exp(1j * phase)
    x = np.fft.ifft(X).real
    x = x * np.sqrt(fs * N)
    return x - np.mean(x)

def compute_moments(y, fs):
    Y = np.fft.rfft(y)
    S = np.abs(Y)**2 / (fs * len(y))
    df = fs / len(y)
    return (np.sum(S)*df, np.sum(omega**2 * S)*df, np.sum(omega**4 * S)*df)

def find_peaks(y):
    peaks = []
    for i in range(1, len(y)-1):
        if y[i] > y[i-1] and y[i] > y[i+1] and y[i] > 0:
            peaks.append(y[i])
    return np.array(peaks)

bandwidths = [0.5, 5.0, 20.0]
labels = ['理想窄带 Δf=0.5 Hz', '中等窄带 Δf=5 Hz', '较宽窄带 Δf=20 Hz']
colors = ['r', 'g', 'b']

np.random.seed(42)
results = []
for Delta_f, label, color in zip(bandwidths, labels, colors):
    S = make_spectrum(Delta_f)
    y = generate_process(S, N, fs)
    m0, m2, m4 = compute_moments(y, fs)
    nu0 = np.sqrt(m2/m0) / (2*np.pi)
    nup = np.sqrt(m4/m2) / (2*np.pi)
    gamma = nu0 / nup
    env = np.abs(hilbert(y))
    peaks = find_peaks(y)
    results.append({
        'Delta_f': Delta_f, 'label': label, 'color': color,
        'y': y, 'env': env, 'peaks': peaks,
        'nu0': nu0, 'nup': nup, 'gamma': gamma, 'sigma': np.sqrt(m0)
    })
    print(f"{label}: ν0={nu0:.2f} Hz, νp={nup:.2f} Hz, γ={gamma:.4f}, 峰值数={len(peaks)}")

# ============================================================
# 3. 单张大图：2行 x 3列 布局
# ============================================================
fig = plt.figure(figsize=(18, 11))
fig.suptitle('窄带高斯过程：从波形到峰值分布', fontsize=16, fontweight='bold')

# ---- 第一行：波形 + 包络 ----
for i, res in enumerate(results):
    ax = fig.add_subplot(2, 3, i+1)
    t_show = t < 2.0
    ax.plot(t[t_show], res['y'][t_show], color=res['color'], linewidth=0.6, alpha=0.7, label='y(t)')
    ax.plot(t[t_show], res['env'][t_show], 'k-', linewidth=1.5, label='包络 A(t)')
    ax.plot(t[t_show], -res['env'][t_show], 'k-', linewidth=1.5)
    ax.set_title(f"{res['label']}\nγ = {res['gamma']:.4f}   ν0 = {res['nu0']:.1f} Hz   νp = {res['nup']:.1f} Hz")
    ax.set_xlabel('时间 (s)')
    ax.set_ylabel('y(t)')
    ax.legend(loc='upper right', fontsize=8)
    ax.grid(True, alpha=0.3)

# ---- 第二行：峰值分布 ----
for i, res in enumerate(results):
    ax = fig.add_subplot(2, 3, i+4)
    peaks = res['peaks']
    sigma = res['sigma']
    ax.hist(peaks, bins=40, density=True, alpha=0.6, color=res['color'], label='峰值直方图')
    a_range = np.linspace(0, 4*sigma, 200)
    ax.plot(a_range, rayleigh.pdf(a_range, scale=sigma), 'k--', linewidth=2, label='瑞利分布')
    ax.set_title(f"峰值分布\nγ = {res['gamma']:.4f}")
    ax.set_xlabel('峰值 a')
    ax.set_ylabel('p_p(a)')
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)

plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()

# ============================================================
# 4. 第二张大图：拆解过程（谱密度 → 波形 → 包络 → 峰值）
# ============================================================
fig2, axes2 = plt.subplots(4, 3, figsize=(16, 14))
fig2.suptitle('拆解过程：谱密度 → 波形 → 包络 → 峰值分布', fontsize=16, fontweight='bold')

for i, res in enumerate(results):
    y = res['y']
    env = res['env']
    peaks = res['peaks']
    sigma = res['sigma']
    color = res['color']
    
    # 4.1 谱密度
    ax = axes2[0, i]
    S = make_spectrum(res['Delta_f'])
    ax.plot(freq, S, color=color, linewidth=1.5)
    ax.set_title(f"① 谱密度 S(ω)\nΔf = {res['Delta_f']} Hz")
    ax.set_xlabel('频率 (Hz)')
    ax.set_ylabel('S(ω)')
    ax.set_xlim(0, 250)
    ax.grid(True, alpha=0.3)
    
    # 4.2 时间波形（1秒）
    ax = axes2[1, i]
    t_show = t < 1.0
    ax.plot(t[t_show], y[t_show], color=color, linewidth=0.7)
    ax.set_title(f"② 时间波形 y(t)\nν0 = {res['nu0']:.1f} Hz")
    ax.set_xlabel('时间 (s)')
    ax.set_ylabel('y(t)')
    ax.grid(True, alpha=0.3)
    
    # 4.3 包络
    ax = axes2[2, i]
    ax.plot(t[t_show], env[t_show], color=color, linewidth=1.5, label='包络 A(t)')
    ax.plot(t[t_show], -env[t_show], color=color, linewidth=1.5)
    ax.fill_between(t[t_show], -env[t_show], env[t_show], color=color, alpha=0.15)
    ax.set_title(f"③ 包络 A(t)\n纺锤状")
    ax.set_xlabel('时间 (s)')
    ax.set_ylabel('A(t)')
    ax.grid(True, alpha=0.3)
    
    # 4.4 峰值分布
    ax = axes2[3, i]
    ax.hist(peaks, bins=40, density=True, alpha=0.6, color=color, label='峰值直方图')
    a_range = np.linspace(0, 4*sigma, 200)
    ax.plot(a_range, rayleigh.pdf(a_range, scale=sigma), 'k--', linewidth=2, label='瑞利分布')
    ax.set_title(f"④ 峰值分布\nγ = {res['gamma']:.4f}")
    ax.set_xlabel('峰值 a')
    ax.set_ylabel('p_p(a)')
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)

plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()