不同频率的穿越率

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

# 设置中文字体，避免中文显示为方块
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# ============================================================
# 1. 参数设置
# ============================================================
fs = 4000.0                  # 采样频率 (Hz)
T = 15.0                     # 总时长 (s)
N = int(fs * T)              # 采样点数
t = np.linspace(0, T, N, endpoint=False)  # 时间轴

# 频率轴（只取正频率部分，rfft 对应）
freq = np.fft.rfftfreq(N, d=1/fs)         # 频率 (Hz)
omega = 2 * np.pi * freq                  # 角频率 (rad/s)

# ============================================================
# 2. 高斯型谱密度（不同中心频率）
# ============================================================
def spectrum_gaussian(f0, Delta_f):
    """
    生成高斯型窄带谱密度。
    参数:
        f0      : 中心频率 (Hz)
        Delta_f : 带宽 (Hz)
    返回:
        S       : 谱密度数组，长度与 omega 相同
    """
    omega0 = 2 * np.pi * f0                # 中心角频率
    sigma_w = np.pi * Delta_f / 3          # 高斯谱的标准差（角频率）
    # 高斯函数：S(ω) = exp(-(ω-ω0)^2 / (2σ_w^2))
    S = np.exp(-(omega - omega0)**2 / (2*sigma_w**2))
    S[omega < 0] = 0                       # 负频率部分置零
    return S

def generate_process(S_omega, N, fs):
    """
    根据给定的单边谱密度 S_omega 生成平稳高斯过程。
    方法：在频域赋予随机相位，再做逆傅里叶变换。
    """
    # 构造双边谱：正频率 + 负频率（对称）
    S_full = np.concatenate([S_omega, S_omega[-2:0:-1]])
    # 随机相位，均匀分布在 [0, 2π)
    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.mean(x)

def count_crossings(y, fs, level=0.0):
    """
    统计正穿越 level 的次数，并返回频率 (Hz)。
    正穿越定义：y[i] < level 且 y[i+1] >= level
    """
    idx = np.where((y[:-1] < level) & (y[1:] >= level))[0]
    return len(idx) / (len(y) / fs)

def compute_theory(S_omega, fs, N):
    """
    由谱密度直接计算理论穿零频率。
    ν0 = (1/2π) * sqrt(m2 / m0)
    其中 m0 = ∫S(ω)dω, m2 = ∫ω^2 S(ω)dω
    """
    m0 = np.sum(S_omega) * (fs / N)        # 零阶谱矩
    m2 = np.sum(omega**2 * S_omega) * (fs / N)  # 二阶谱矩
    return (1/(2*np.pi)) * np.sqrt(m2 / m0)

np.random.seed(42)            # 固定随机种子，保证结果可复现

# ============================================================
# 3. 不同中心频率，同一高斯谱型
# ============================================================
f0_list = [30, 50, 80, 120, 180, 250]   # 6 个中心频率 (Hz)
Delta_f = 10.0                           # 带宽固定为 10 Hz
colors = plt.cm.viridis(np.linspace(0, 1, len(f0_list)))  # 颜色渐变

results = []                             # 存储每个频率的结果
for f0, color in zip(f0_list, colors):
    S = spectrum_gaussian(f0, Delta_f)   # 生成谱密度
    y = generate_process(S, N, fs)       # 生成时域过程
    nu0_emp = count_crossings(y, fs, 0.0)  # 实测穿零频率
    nu0_theory = compute_theory(S, fs, N)  # 理论穿零频率
    results.append({
        'f0': f0, 'S': S, 'y': y, 'color': color,
        'nu0_emp': nu0_emp, 'nu0_theory': nu0_theory
    })
    # 打印对比结果
    print(f"f0={f0:3d} Hz: ν0实测={nu0_emp:6.2f} Hz, ν0理论={nu0_theory:6.2f} Hz")

# ============================================================
# 4. 单张大图：每个频率一个子图
# ============================================================
fig = plt.figure(figsize=(20, 14))
fig.suptitle('不同中心频率的高斯谱：穿越率变化', fontsize=18, fontweight='bold', y=0.98)

# GridSpec: 3 行 6 列，hspace/wspace 控制子图间距
gs = fig.add_gridspec(3, 6, hspace=0.5, wspace=0.4,
                      left=0.05, right=0.98, top=0.93, bottom=0.05)

# ---- 第一行：6个谱密度，各自一个子图 ----
for i, res in enumerate(results):
    ax = fig.add_subplot(gs[0, i])
    ax.plot(freq, res['S'], color=res['color'], linewidth=1.8)     # 画谱密度曲线
    ax.fill_between(freq, 0, res['S'], color=res['color'], alpha=0.2)  # 填充下方区域
    ax.axvline(res['f0'], color='k', linestyle='--', linewidth=1.0, alpha=0.5)  # 标注中心频率
    ax.set_title(f'谱密度 f0={res["f0"]} Hz', fontsize=10)
    ax.set_xlabel('频率 (Hz)')
    ax.set_ylabel('S(ω)')
    ax.set_xlim(0, 300)                  # 限制频率范围
    ax.grid(True, alpha=0.3)

# ---- 第二行：6个时间波形，各自一个子图 ----
for i, res in enumerate(results):
    ax = fig.add_subplot(gs[1, i])
    t_show = t < 0.3                     # 只显示前 0.3 秒，便于观察波形
    ax.plot(t[t_show], res['y'][t_show], color=res['color'], linewidth=0.8)
    ax.set_title(f'波形 f0={res["f0"]} Hz\nν0={res["nu0_emp"]:.1f} Hz', fontsize=10)
    ax.set_xlabel('时间 (s)')
    ax.set_ylabel('y(t)')
    ax.grid(True, alpha=0.3)

# ---- 第三行：穿越率对比 ----
# 3.1 穿越率柱状图（占 2 列）
ax = fig.add_subplot(gs[2, 0:2])
f0_arr = [r['f0'] for r in results]
nu0_emp_arr = [r['nu0_emp'] for r in results]
nu0_theory_arr = [r['nu0_theory'] for r in results]
x = np.arange(len(f0_arr))
width = 0.35
ax.bar(x - width/2, nu0_emp_arr, width, label='实测', color='steelblue')
ax.bar(x + width/2, nu0_theory_arr, width, label='理论', color='coral')
ax.set_xticks(x)
ax.set_xticklabels([f'{f}' for f in f0_arr])
ax.plot(x, f0_arr, 'k--o', markersize=6, label='ν0 = f0')   # 理想线
ax.set_title('③ 穿越率对比', fontsize=12)
ax.set_xlabel('中心频率 f0 (Hz)')
ax.set_ylabel('ν0⁺ (Hz)')
ax.legend()
ax.grid(True, alpha=0.3)

# 3.2 穿越率 vs 中心频率（散点 + 拟合线，占 2 列）
ax = fig.add_subplot(gs[2, 2:4])
ax.plot(f0_arr, nu0_emp_arr, 'bo-', markersize=10, label='实测')
ax.plot(f0_arr, nu0_theory_arr, 'r--', linewidth=2, label='理论')
ax.plot(f0_arr, f0_arr, 'k:', linewidth=2, label='ν0 = f0')  # 理想线
ax.set_title('④ 穿越率 vs 中心频率', fontsize=12)
ax.set_xlabel('中心频率 f0 (Hz)')
ax.set_ylabel('ν0⁺ (Hz)')
ax.legend()
ax.grid(True, alpha=0.3)

# 3.3 归一化穿越率（占 2 列）
ax = fig.add_subplot(gs[2, 4:6])
ratio = np.array(nu0_emp_arr) / np.array(f0_arr)   # ν0 / f0
ax.plot(f0_arr, ratio, 'bo-', markersize=10)
ax.axhline(1.0, color='r', linestyle='--', linewidth=2, label='ν0/f0 = 1')  # 理想值
ax.set_title('⑤ 归一化穿越率 ν0/f0', fontsize=12)
ax.set_xlabel('中心频率 f0 (Hz)')
ax.set_ylabel('ν0⁺ / f0')
ax.set_ylim(0.9, 1.1)                  # 限制纵轴范围，突出偏离程度
ax.legend()
ax.grid(True, alpha=0.3)

plt.show()