08-EFT.WP.Core.Sea v1.0 | 第5章 频谱分析与特征

目录 ／ 文档-技术白皮书 ／ 08-EFT.WP.Core.Sea v1.0

第5章 频谱分析与特征

I. 目标与范围

• 给出从离散时域 x[n] 到功率谱密度 S_xx(f) 的估计方法族（periodogram/Welch/multitaper），统一单位与归一化。
• 定义特征族（峰值、带能量、谱心、谱平坦度、谐波簇、相干与相位）及稳健性评估，与第2章采样与量化、第4章 H(f) 调理、第3章 tau_mono/ts 时间基准保持一致。
• 输出 P85-* 公设、S85-* 最小方程与 Mx-5 频谱分析流程；与 I80-4 fft/psd/feature_extract 接口对齐。

II. 基本对象与单位约定

1. DFT 与频率轴
   • X[k] = ∑_{n=0}^{N-1} x[n] * w[n] * exp( -j * 2 * pi * k * n / N )。
   • f[k] = k * fs / N，Delta_f = fs / N，k ∈ {0,...,N-1}。
   • 窗函数功率归一化量 U_w = ( 1 / N ) * ∑_{n=0}^{N-1} w[n]^2；相干增益 G_c = ( 1 / N ) * ∑_{n=0}^{N-1} w[n]。
2. 谱的单位与一/双边
   • 单边 PSD：S_xx^1(f) ≈ 2 * S_xx^2(f)（除 f=0 与 f=fs/2）；单位保持为 unit(x)^2 / Hz。
   • 振幅谱密度 ASD(f) = sqrt( S_xx(f) )，单位 unit(x) / sqrt(Hz)。
   • Parseval 一致性：∑_{n} |x[n]|^2 / N ≈ ∑_{k} S_xx(f[k]) * Delta_f（在正确归一化下成立）。

III. 公设 P85-*（频谱估计一致性）

• P85-1 单位与归一化公设
  任何 PSD 报告必须显式给出 {one_sided|two_sided}、U_w 与是否乘以 2 的处理；特征计算基于相同口径。
• P85-2 泄漏与窗函数公设
  使用窗函数时，必须给出主瓣宽与旁瓣抑制 A_sl_dB；目标带内峰值测量需保证峰落在主瓣内或采用插值修正。
• P85-3 统计独立性公设
  Welch 或多 taper 平均时，必须报告有效独立平均数 K_eff 或等效自由度 nu_eff；所有方差界限以此为准。
• P85-4 时频一致公设
  采用 STFT 时，需声明窗长 L 与步长 H，并记录 Delta_t ≈ L / fs、Delta_f ≈ fs / L；时间对齐使用 tau_mono，发布用 ts。
• P85-5 到达时与相位公设
  涉及相位或群时延的特征，必须扣除第4章 H(f) 的相位响应；否则标注为未补偿并给出可能偏差。

IV. 最小方程 S85-*（估计与界限）

1. S85-1 Periodogram 归一化
   P_xx[k] = | X[k] |^2 / ( N * fs * U_w )（双边）；单边修正为 P_xx^1[k] = 2 * P_xx[k]，k ∈ (0, N/2)。
2. S85-2 Welch 平均与方差
   将序列分为 K 段、每段长度 N_s、重叠 o ∈ [0,1)：
   • S_hat_xx(f[k]) = ( 1 / K ) * ∑_{i=1}^{K} | X_i[k] |^2 / ( N_s * fs * U_w )。
   • 等效自由度（独立段）nu_ind ≈ 2 * K * ( ( ∑ w^2 )^2 / ∑ w^4 )；考虑重叠相关性 nu_eff = nu_ind / c_ovlp，c_ovlp >= 1。
   • 方差近似：Var[ S_hat_xx(f) ] ≈ 2 * S_xx(f)^2 / nu_eff。
3. S85-3 Multitaper 估计
   以 P 个 DPSS taper：S_hat_xx(f) = ( 1 / P ) * ∑_{p=1}^{P} | FFT( v_p * x ) |^2 / ( N * fs * U_{v_p} )，nu_eff ≈ 2 * P。
4. S85-4 谱峰插值（抛物插值）
   设功率谱在 k-1,k,k+1 处为 P_- , P_0 , P_+：
   • delta = ( P_+ - P_- ) / ( 2 * ( 2 * P_0 - P_- - P_+ ) )，
   • f_peak ≈ ( k + delta ) * fs / N；幅值修正 A_peak ≈ P_0 - 0.25 * ( P_- - P_+ ) * delta。
5. S85-5 带能量与带宽
   • BandPower( f1,f2 ) = ∑_{k: f1<=f[k]<=f2} S_hat_xx(f[k]) * Delta_f。
   • f_centroid = ( ∑ f[k] * S_hat_xx[k] ) / ( ∑ S_hat_xx[k] )，
     SpecBandwidth = sqrt( ( ∑ ( f[k] - f_centroid )^2 * S_hat_xx[k] ) / ( ∑ S_hat_xx[k] ) )。
6. S85-6 谱平坦度与调谐度
   • SF = exp( mean( ln S_hat_xx[k] ) ) / mean( S_hat_xx[k] )。
   • Tonalness ≈ 1 - SF（启发式，范围 [0,1]）。
7. S85-7 交叉谱与相干
   • S_xy(f) = E[ X(f) * Y(f)^* ]，Coh^2(f) = | S_xy(f) |^2 / ( S_xx(f) * S_yy(f) ) ∈ [0,1]。
   • phi_xy(f) = arg S_xy(f)；在 nu_eff 自由度下，Coh^2 的零假设显著性阈值可用 alpha 近似：Coh_alpha^2 ≈ 1 - alpha^( 1 / ( K_eff - 1 ) )（独立段近似）。

V. 窗函数与泄漏管理

1. 主瓣宽与旁瓣抑制
   • Hamming/Hann：较低旁瓣，主瓣稍宽；Blackman：旁瓣更低；Kaiser：以 beta 调控 A_sl_dB 与主瓣宽折中。
   • 选择规则：峰检优先低旁瓣窗；能量守恒与 PSD 平滑优先适中主瓣与较小 ENBW。
2. ENBW 与分辨率
   • ENBW = fs * ( ∑ w[n]^2 ) / ( ( ∑ w[n] )^2 )；等效分辨率带宽 RBW ≈ ENBW。
   • 峰间可分辨性约束：| f_a - f_b | >= RBW。

VI. 特征族定义（F85-*）

1. F85-1 峰值与谐波
   • f_peak、A_peak 按 S85-4；基础频 f0 通过扫描最大化 ∑_{h=1}^{H} S_hat_xx( h * f0 ) 得到。
   • HarmonicRatio = ( ∑_{h=2}^{H} S_hat_xx( h * f0 ) ) / S_hat_xx( f0 )。
2. F85-2 带能量与信噪比
   SNR_band = 10 * log10( BandPower(signal_band) / BandPower(noise_band) )；噪声带不含已检测峰附近 ±RBW。
3. F85-3 谱心与带宽
   f_centroid 与 SpecBandwidth 按 S85-5；可在多带上分别计算以作指纹。
4. F85-4 平坦度与粗糙度
   SF 按 S85-6；Roughness ≈ Var( diff( log S_hat_xx[k] ) )（离散近似）。
5. F85-5 相干与相位差
   Coh_band = mean_{f ∈ band}( Coh^2(f) )，Phase_band = atan2( ∑ Im S_xy , ∑ Re S_xy )。Coh^2(f)、phi_xy(f)；可在感兴趣频带取分位数或带内平均：
6. F85-6 容错与缺失
   若存在缺失掩码 m[n] ∈ {0,1}，对每段应用相同窗并以有效样本数归一化 U_w，并在清单记录缺失率。

VII. 时频分析与跟踪

1. STFT 定义
   • STFT[m,k] = ∑_{n} x[n] * w[n - m*H] * exp( -j * 2 * pi * k * n / N_w )。
   • 谱图 Spec[m,k] = | STFT[m,k] |^2 / ( N_w * fs * U_w )。
2. 时频分辨率
   Delta_t ≈ N_w / fs，Delta_f ≈ fs / N_w，满足 Delta_t * Delta_f ≈ 1；随 H 控制时间采样。
3. 轨迹与漂移
   峰轨迹 f_peak[m] 由 S85-4 的时窗内插得到；以 drift = max | f_peak[m+1] - f_peak[m] | 作为漂移速率特征；与第9章 drift 监测一致。

VIII. 稳健性与不确定度

• 振幅偏置校正
  单音幅值估计需除以相干增益：A_tone ≈ ( 2 / N ) * | X[k0] | / G_c（单边）。
• 方差与置信区间
  CI_alpha(f) ≈ S_hat_xx(f) * [ nu_eff / χ^2_{1-α/2}(nu_eff) , nu_eff / χ^2_{α/2}(nu_eff) ]（卡方近似，独立段）。
• 交叉验证
  通过改变 N_s/overlap/window 计算 S_hat_xx 的敏感度；若差异超过门限，则提升 K 或改用 multitaper。

IX. 执行流程 Mx-5（频谱分析与特征）

• 设定目标：fs/Delta_f/RBW/K or P/one_sided 与感兴趣带集合 {band_i}；从第4章取得 H(f) 相位补偿方案。
• 选择估计器：method ∈ {"periodogram","welch","multitaper"}，设 N or N_s/overlap（或 P）。
• 窗与归一化：选择 w[n]，计算 U_w/G_c/ENBW；明确单边/双边口径与单位。
• 计算 PSD：得到 S_hat_xx(f)（及 S_hat_yy(f)、S_hat_xy(f) 如需）；估计 nu_eff 与方差。
• 特征提取：对每个 band 计算 BandPower/SNR_band/f_peak/f_centroid/SF/Coh^2/phi_xy 等。
• 稳健性评估：改变 N_s/overlap/window 或 P，检查特征相对变化率 <= epsilon_feat。
• 清单与发布：记录参数与版本，输出特征表与不确定度列；对与 T_arr 相关的相位特征，标注补偿状态。

X. 接口绑定与返回形制（I80-4）

1. fft(sig:any, window:str="hann") -> complex[]
   返回 X[k]；由调用方完成单位归一化与单边转换；建议同时返回 G_c/U_w。
2. psd(sig:any, method:str="welch", seg:int=8, overlap:float=0.5) -> dict
   约定返回键：{"f":float[], "Sxx":float[], "one_sided":bool, "U_w":float, "G_c":float, "nu_eff":float, "window":str, "overlap":float}。
3. feature_extract(sig:any, feats:list[str]) -> dict
   • 支持键："f_peak","A_peak","bandpower:{f1,f2}","snr:{sig,noise}","f_centroid","spec_bw","sf","coh:{y,band}","phi:{y,band}"。
   • 结果需携带 {"uncertainty":..., "nu_eff":..., "manifest_ref":...}。

XI. 清单（manifest）最小字段

• {method, window, N or N_s, overlap, one_sided, U_w, G_c, ENBW, nu_eff, RBW, bands:[{f1,f2}], H_phase_comp:{"enabled":bool,"ver":str}, ts_range, missing_rate}。
• 特征表需引用该清单并记录 epsilon_feat、CI_alpha 参数与 alpha。

XII. 互锁与跨卷引用

• 与第2章：ENOB/DR/Delta_f 的一致性检查；RBW 不得小于物理可分辨极限。
• 与第3章：tau_mono/ts 的时间基准在 STFT 索引中保持单调；抖动 J 的 SNR_jitter 对高频谱峰的影响按第2章与第4章合并评估。
• 与第4章：任何相位与群时延相关特征必须使用 H(f) 补偿。
• 与第8章：到达时 T_arr 的频域估计需在同一 gamma(ell) 与 c_ref/n_eff 口径下报告。