一种自适应窗口长度与角度优化的时频脊线追踪算法（Python）

算法特点

自适应：算法根据信号特性自动调整分析参数

窗口长度与角度优化：同时优化时频分析的窗口长度和调频变化率角度

时频脊线追踪：核心目标是提取信号在时频平面上的能量脊线

多参数协同优化：通过对比度指标协同优化多个分析参数

算法流程

开始
│
├─ 信号预处理
│   ├─ 计算信号长度Ns和总时长tend
│   ├─ 确定Nyquist频率fend = fs/2
│   └─ 对信号进行零扩展：前后各补Ns个零
│
├─ 参数网格生成
│   ├─ 角度网格Ana：[-π/2, π/2]区间Nna等分
│   ├─ 窗口长度网格Awl：[tend/(Nwl+1), tend]区间Nwl等分
│   ├─ 频率网格Af：[fend/Nf, fend]区间Nf等分
│   ├─ 时间采样点At：时间轴上1/fs间隔点
│   └─ 脊线时间点Atao：[tend/Ntao, tend]区间Ntao等分
│
├─ 四维时频张量初始化 (Ntao×Nf×Nwl×Nna)
│
├─ 主计算循环 (遍历所有窗口长度和角度)
│   ├─ 计算当前窗口参数
│   │   ├─ 半窗长Nhalfwin = round(当前窗口长度/2 * fs)
│   │   └─ 全窗长Nwin = 2*Nhalfwin + 1
│   │
│   ├─ 信号截取
│   │   └─ 在扩展信号上截取以脊线点为中心的Nwin长度段
│   │
│   ├─ 构造高斯窗函数
│   │   ├─ 窗函数：exp(-0.5*(归一化时间偏移/Gc)^2)
│   │   └─ 窗函数归一化
│   │
│   ├─ 计算窗内时间偏移向量
│   │
│   ├─ 角度有效性检查
│   │   └─ 若|tan(当前角度)| * 窗口长度 > 总时长 → 标记为NaN
│   │
│   ├─ 有效角度处理
│   │   ├─ 计算线性调频核kernela：(1 + 角度因子 * 相对时间)
│   │   ├─ 计算时间项矩阵：t_center + t_rel + (角度因子 * t_rel²)
│   │   ├─ 计算相位补偿核kernelb：exp(-j2π * 时间项矩阵 * 频率)
│   │   └─ 计算时频表示：截取信号 * 窗 * kernela 与 kernelb 的加权和
│   │
│   └─ 结果存入四维张量
│
├─ 计算对比度指标CI
│   └─ CI = E[|TFR|⁴] / (E[|TFR|²])² (沿频率轴计算)
│
├─ 确定最优参数
│   ├─ 对每个时间点：寻找使CI最大的(窗口长度,角度)组合
│   ├─ 提取最优窗口长度dwl
│   └─ 提取最优角度dna
│
├─ 提取最终时频表示TFR
│   └─ 从四维张量中选取dwl和dna对应的二维切片
│
└─ 返回(TFR, dwl, dna, Atao, Af)
结束

此算法适用于分析具有复杂频率变化规律的非平稳信号，如机械振动、生物医学信号和通信信号等。通过动态优化局部时频分析的窗口参数，显著提升对非线性调频信号的时频聚集性。

for iwl in range(Nwl):
        ipwl = Awl[iwl]  # 当前窗长
        Nhalfwin = int(round(ipwl * fs / 2))  # 半窗长点数
        Nwin = 2 * Nhalfwin + 1               # 全窗长点数

        # 构造分析窗 (高斯窗)
        win_x = np.linspace(-1, 1, Nwin) * (2 * Gc)
        window = np.exp(-0.5 * win_x**2)
        window /= np.sum(window)  # 归一化

        # 截取信号段 (所有分析时间点)
        start_idx = Ns + Itao - Nhalfwin  # 在扩展信号中的起始索引
        idx_matrix = start_idx[:, None] + np.arange(Nwin)  # 索引矩阵
        struc = sextend[idx_matrix]  # 截取信号段 (Ntao×Nwin)

        # 窗内时间向量 (中心点为0)
        t_win = (np.arange(Nwin) - Nhalfwin) / fs

        for ina in range(Nna):
            tanna = np.tan(Ana[ina])  # 当前角度正切值

            # 检查角度是否有效
            if np.abs(tanna) * ipwl > tend:
                continue  # 超出范围跳过

            # 幅度校正因子
            kernela = 1 + (tanna / tend) * 2 * t_win
            kernela /= np.mean(kernela)  # 归一化

            # 相位校正因子 (核函数)
            # 修正关键行：确保相位参数与频率向量维度匹配
            phase_arg = t_win + (tanna / tend) * t_win**2
            kernelb = np.exp(-1j * 2 * np.pi * np.outer(phase_arg, Af))

            # 计算子时频表示
            # 修正关键行：确保矩阵乘法维度匹配
            weighted_sig = struc * window * kernela
            subTFRs[:, :, iwl, ina] = weighted_sig @ kernelb

完整代码通过知乎学术咨询获得（哥廷根数学学派）：

工学博士，担任《Mechanical System and Signal Processing》审稿专家，《中国电机工程学报》，《控制与决策》，《系统工程与电子技术》，《电力系统保护与控制》，《宇航学报》等EI期刊审稿专家。擅长领域：现代信号处理，机器学习，深度学习，数字孪生，时间序列分析，设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。