Modbus-rtu/HostComputer/bin_to_time_freq.py

# -*- coding: utf-8 -*-
"""
.bin 波形文件 -> 时频曲线查看/出图工具（PLSR 项目用）

数据格式: 逻辑分析仪导出的单通道数字采样，每字节 1 个采样点；
          电平 >= 阈值(默认128) 判为高（脉冲），否则为低。
频率定义: 用连续 N 个脉冲的上升沿间隔估计频率；
          频率点放在 N 周期测量窗口中间；
          频率 = N * 采样率 /（第 N 个后续上升沿 - 当前上升沿）；
时间定义: 频率点时间为测量窗口中点；高电平宽度取中心脉冲。

用法:
    python bin_to_time_freq.py <bin文件> [采样率Hz] [选项]
    未填写采样率时默认使用 3125000 Hz。
示例:
    python bin_to_time_freq.py "Document/PLSR_document/波形/10段.bin" 3125000
    python bin_to_time_freq.py xxx.bin 3125000 --ymax=12000
    python bin_to_time_freq.py xxx.bin 3125000 --save=out.png
    python bin_to_time_freq.py xxx.bin 3125000 --selftest
    python bin_to_time_freq.py xxx.bin 3125000 --reverse --cycles=8
    python bin_to_time_freq.py xxx.bin 3125000 --d1 --d2
    python bin_to_time_freq.py xxx.bin 3125000 --d1 --dwindow-ms=10
    python bin_to_time_freq.py xxx.bin 3125000 --d1 --d2 --d1-window-ms=5 --d2-window-ms=15

选项:
    --thresh=128   电平阈值（默认128）
    --cycles=8     用连续多少个周期估计频率（默认8；1为逐周期）
    --reverse      反相电平后再检测脉冲（低电平作为脉冲高电平）
    --ymax=12000   频率轴上限（默认自适应）
    --d1           另开一张图：频率对时间的一阶导数 df/dt（红）
    --d2           另开一张图：频率对时间的二阶导数 d²f/dt²（黑）
    --dwindow-ms=5 一、二阶导数的公共拟合窗口，单位 ms（默认5）
    --d1-window-ms 原函数 -> 一阶导数的拟合窗口（默认继承 --dwindow-ms）
    --d2-window-ms 一阶导数 -> 二阶导数的拟合窗口（默认继承 --dwindow-ms）
    --save=路径    保存 PNG 后退出（不弹窗）
    --selftest     仅打印统计，不画图
    --help, -h     显示这份中文帮助后退出

交互（弹窗模式）:
    滚轮           缩放 X 轴（向上放大/向下缩小，以鼠标位置为中心）
    Ctrl+滚轮      缩放 Y 轴（以鼠标 Y 位置为中心）
    左键拖拽       双向平移（上下左右跟随鼠标）
    双击           复位到全图
    鼠标悬停       高亮最近的点并显示其测量点序号、时间与频率

依赖: pip install numpy matplotlib
"""

import sys
import os

import numpy as np

try:
    import matplotlib
    matplotlib.rcParams['font.sans-serif'] = ['Microsoft YaHei', 'SimHei']
    matplotlib.rcParams['axes.unicode_minus'] = False
except Exception:
    pass


def load_time_freq(path, fs, threshold=128, cycles=8, reverse=False):
    """读取 bin，返回 (t_sec, freq_hz, hi_width_sec, meta)。
    t_sec 以第一个脉冲高电平中点为 0 时刻；hi_width 为中心脉冲高电平宽度。"""
    data = np.fromfile(path, dtype=np.uint8)
    if data.size == 0:
        raise ValueError("文件为空: %s" % path)
    if not 0 <= threshold <= 255:
        raise ValueError("阈值必须在0到255之间")
    samples = (data >= threshold).astype(np.int8)
    if reverse:
        samples = 1 - samples

    # run 级压缩
    changes = np.flatnonzero(np.diff(samples) != 0) + 1
    starts = np.concatenate(([0], changes))
    ends = np.concatenate((changes, [len(samples)]))
    runs = np.column_stack((starts, ends, samples[starts]))

    # 每个高电平 run 的起点都是一个可用于测周期的上升沿。
    # 若采集从高电平开始，第一个高电平起点不是被观测到的上升沿，丢弃它。
    hi = np.flatnonzero(runs[:, 2] == 1)
    if hi.size and runs[hi[0], 0] == 0:
        hi = hi[1:]

    if hi.size == 0:
        raise ValueError("未检测到上升沿脉冲: %s" % path)

    # 保留浮点上升沿位置：在原始电平相邻采样点之间线性插值，
    # 再配合多周期测量，显著降低整数采样点带来的量化抖动。
    rise_idx = runs[hi, 0].astype(np.int64)
    t_start = rise_idx.astype(np.float64)
    valid = rise_idx > 0
    if np.any(valid):
        y0 = data[rise_idx[valid] - 1].astype(np.float64)
        y1 = data[rise_idx[valid]].astype(np.float64)
        delta = y1 - y0
        frac = np.zeros_like(y0)
        np.divide(float(threshold) - y0, delta, out=frac, where=delta != 0)
        t_start[valid] = rise_idx[valid] - 1.0 + np.clip(frac, 0.0, 1.0)

    hi_width = (runs[hi, 1] - runs[hi, 0]).astype(np.float64)

    if fs <= 0:
        raise ValueError("采样率必须大于0")
    try:
        cycles = int(cycles)
    except (TypeError, ValueError):
        raise ValueError("cycles 必须是正整数")
    if cycles < 1:
        raise ValueError("cycles 必须是正整数")
    if len(t_start) < 2:
        raise ValueError("脉冲不足两个，无法计算频率")

    # 用下降沿的阈值 crossing 计算高电平宽度，时间单位仍为采样点。
    # 文件末尾若停在高电平，频率仍可用该上升沿作为前一窗口终点；
    # 该脉冲的宽度只取到文件末尾，不用于猜测频率。
    has_fall = np.zeros(len(hi), dtype=bool)
    adjacent = (hi + 1) < len(runs)
    has_fall[adjacent] = runs[hi[adjacent] + 1, 2] == 0
    fall_idx = np.where(has_fall, runs[hi, 1], len(data)).astype(np.int64)
    fall_edge = fall_idx.astype(np.float64)
    valid = has_fall & (fall_idx > 0) & (fall_idx < len(data))
    if np.any(valid):
        y0 = data[fall_idx[valid] - 1].astype(np.float64)
        y1 = data[fall_idx[valid]].astype(np.float64)
        delta = y1 - y0
        frac = np.zeros_like(y0)
        np.divide(float(threshold) - y0, delta, out=frac, where=delta != 0)
        fall_edge[valid] = fall_idx[valid] - 1.0 + np.clip(frac, 0.0, 1.0)
    hi_width = np.maximum(fall_edge - t_start, 0.0)
    pulse_mid = t_start + hi_width / 2.0

    # 单周期只看到整数采样点，周期接近半个采样点时会严重跳变。
    # 跨 cycles 个周期计算总时长，再换算回单周期频率，可显著抑制量化抖动。
    cycles = min(cycles, len(t_start) - 1)
    left = np.arange(len(t_start) - cycles, dtype=np.int64)
    right = left + cycles
    period = t_start[right] - t_start[left]
    freq = cycles * float(fs) / np.maximum(period, np.finfo(np.float64).eps)

    # 频率点放在测量窗口的时间中心；宽度取中心附近脉冲，仅用于兼容绘图接口。
    t_center = (t_start[left] + t_start[right]) / 2.0
    center_idx = np.minimum(left + cycles // 2, len(hi_width) - 1)
    hi_width_out = hi_width[center_idx]
    # 时间原点仍为第一个脉冲的高电平中点；多周期测量点自然位于窗口中心。
    t_sec = (t_center - pulse_mid[0]) / fs

    meta = {
        'total_samples': int(len(samples)),
        'duration': len(samples) / fs,
        'pulses': int(len(freq)),
        'detected_pulses': int(len(t_start)),
        'cycles': int(cycles),
        'offset_ms': t_start[0] / fs * 1000,
        'activity_duration': max(0.0, (fall_edge[-1] - pulse_mid[0]) / fs),
    }
    return t_sec, freq, hi_width_out / fs, meta


def selftest(path, fs, thresh, cycles, reverse=False):
    import statistics
    t, freq, hi_w, meta = load_time_freq(path, fs, thresh, cycles, reverse)
    print("文件: %s" % path)
    print("总采样: %d  (%.3f s @ %.2f MS/s)" % (meta['total_samples'], meta['duration'], fs / 1e6))
    print("起始采集偏移: %.1f ms" % meta['offset_ms'])
    print("检测到脉冲: %d，频率测量点: %d" % (
        meta['detected_pulses'], meta['pulses']))
    if meta['pulses']:
        print("频率 min=%.0f max=%.0f 中位=%.0f Hz" % (
            freq.min(), freq.max(), statistics.median(freq)))
        print("末测量点: t=%.3f ms, 高电平 %.3f ms, %.0f Hz" % (
            t[-1] * 1000, hi_w[-1] * 1000, freq[-1]))
        print("估计周期窗口: %d 个周期" % meta['cycles'])
        print("波形活动时长: %.1f ms" % (meta['activity_duration'] * 1000))


def local_linear_derivative(t_sec, values, window_s=5e-3, min_points=7):
    """对非等间隔时间点做局部线性拟合，返回 d(values)/dt。

    直接对相邻频率点差分会放大采样量化抖动；局部最小二乘斜率
    使用真实时间间隔，并对这种抖动做平均。
    """
    t = np.asarray(t_sec, dtype=np.float64)
    y = np.asarray(values, dtype=np.float64)
    if t.ndim != 1 or y.ndim != 1 or len(t) != len(y):
        raise ValueError("t_sec 和 values 必须是等长一维数组")
    if len(t) < 2:
        raise ValueError("至少需要两个点才能计算导数")
    if window_s <= 0:
        raise ValueError("导数拟合窗口必须大于0")
    if np.any(np.diff(t) <= 0):
        raise ValueError("时间点必须严格递增")

    n = len(t)
    min_points = min(n, max(2, int(min_points)))
    half_window = window_s / 2.0
    left = np.searchsorted(t, t - half_window, side='left')
    right = np.searchsorted(t, t + half_window, side='right')

    # 低频段在固定时间窗口内可能点数太少，至少补足 min_points 个点。
    idx = np.arange(n)
    fallback_left = np.clip(idx - min_points // 2, 0, n - min_points)
    fallback_right = fallback_left + min_points
    too_few = (right - left) < min_points
    left[too_few] = fallback_left[too_few]
    right[too_few] = fallback_right[too_few]

    # 移动时间原点减少长时采集时前缀和相减的精度损失。
    x = t - (t[0] + t[-1]) / 2.0
    sx = np.concatenate(([0.0], np.cumsum(x)))
    sy = np.concatenate(([0.0], np.cumsum(y)))
    sxx = np.concatenate(([0.0], np.cumsum(x * x)))
    sxy = np.concatenate(([0.0], np.cumsum(x * y)))
    count = (right - left).astype(np.float64)
    sum_x = sx[right] - sx[left]
    sum_y = sy[right] - sy[left]
    sum_xx = sxx[right] - sxx[left]
    sum_xy = sxy[right] - sxy[left]
    denominator = count * sum_xx - sum_x * sum_x
    numerator = count * sum_xy - sum_x * sum_y

    derivative = np.empty(n, dtype=np.float64)
    good = np.abs(denominator) > np.finfo(np.float64).eps
    derivative[good] = numerator[good] / denominator[good]
    derivative[~good] = np.gradient(y, t)[~good]
    return derivative


def plot_show(path, fs, thresh, ymax, cycles, reverse, save_path=None,
              draw_d1=False, draw_d2=False, d1_window_ms=5.0,
              d2_window_ms=None):
    t, freq, hi_w, meta = load_time_freq(path, fs, thresh, cycles, reverse)

    import matplotlib.pyplot as plt
    from matplotlib.ticker import MaxNLocator

    fig, ax = plt.subplots(figsize=(15, 7))
    fig.subplots_adjust(bottom=0.10, top=0.92)

    # 全脉冲散点：复用 artist，重绘时按可见范围抽稀（大文件流畅）
    MAX_VISIBLE_POINTS = 5000   # 可见范围内最多绘制的点数，超过则等间隔抽稀
    (scatter,) = ax.plot([], [], '.', ms=1.5, color='C0', alpha=0.5, zorder=1)
    tx_ms = t * 1000
    x_full = tx_ms
    y_full = freq
    ax.set_xlim(0.0, meta['activity_duration'] * 1000 * 1.02)
    ax.set_ylim(0, ymax if ymax > 0 else freq.max() * 1.08)
    ax.set_title("%s 时频曲线（检测到脉冲 %d，频率测量点 %d，%d 周期平均，活动时长 %.0f ms）" % (
        os.path.basename(path), meta['detected_pulses'], meta['pulses'], meta['cycles'],
        meta['activity_duration'] * 1000))
    ax.set_xlabel("时间 (ms)")
    ax.set_ylabel("频率 (Hz)")
    # 更密的刻度
    ax.xaxis.set_major_locator(MaxNLocator(nbins=20))
    ax.yaxis.set_major_locator(MaxNLocator(nbins=15))
    ax.grid(True, which='major', alpha=0.35)
    ax.grid(True, which='minor', alpha=0.15)
    ax.minorticks_on()

    fig_der = None
    draw_d1 = draw_d1 and (len(freq) >= 2)
    draw_d2 = draw_d2 and (len(freq) >= 3)
    if draw_d1 or draw_d2:
        f64 = freq.astype(np.float64)
        if d2_window_ms is None:
            d2_window_ms = d1_window_ms
        d1_full = local_linear_derivative(t, f64, d1_window_ms * 1e-3)
        d2_full = (local_linear_derivative(t, d1_full, d2_window_ms * 1e-3)
                   if draw_d2 else None)
        plot_count = int(draw_d1) + int(draw_d2)
        fig_der, axes_der = plt.subplots(plot_count, 1, figsize=(15, 4.5 * plot_count),
                                         sharex=True, squeeze=False)
        fig_der.subplots_adjust(bottom=0.14, top=0.90)
        row = 0
        if draw_d1:
            ax_der = axes_der[row, 0]
            ax_der.plot(x_full, d1_full, '-', lw=0.9, color='red', label='df/dt')
            ax_der.set_ylabel("一阶导数 (Hz/s)")
            ax_der.legend(loc='upper right')
            row += 1
        if draw_d2:
            ax_der = axes_der[row, 0]
            ax_der.plot(x_full, d2_full, '-', lw=0.9, color='black', label='d²f/dt²')
            ax_der.set_ylabel("二阶导数 (Hz/s²)")
            ax_der.legend(loc='upper right')
        for ax_der in axes_der[:, 0]:
            ax_der.set_xlim(0.0, meta['activity_duration'] * 1000 * 1.02)
            ax_der.grid(True, which='major', alpha=0.35)
            ax_der.minorticks_on()
        axes_der[-1, 0].set_xlabel("时间 (ms)")
        axes_der[0, 0].set_title(
            "时频曲线导数（一阶窗口 %.3g ms / 二阶窗口 %.3g ms）" %
            (d1_window_ms, d2_window_ms))

    if save_path:
        scatter.set_data(x_full, y_full)
        root, ext = os.path.splitext(save_path)
        fig.savefig(save_path, dpi=130)
        print("已保存: %s" % save_path)
        if fig_der is not None:
            der_path = root + "_d" + ext
            fig_der.savefig(der_path, dpi=130)
            print("已保存: %s" % der_path)
        plt.close(fig)
        if fig_der is not None:
            plt.close(fig_der)
        return

    # 交互：滚轮缩放 X（Ctrl+滚轮缩放 Y）/ 左键双向拖拽平移 / 双击复位 / 悬停高亮
    state = {'press_x': None, 'press_y': None,
             'press_xlim': None, 'press_ylim': None}

    # 悬停高亮：一个红点标记 + 一个带框文本
    (hl_marker,) = ax.plot([], [], 'o', ms=9, mfc='red', mec='white',
                           mew=1.0, zorder=5, visible=False)
    hl_text = ax.text(0, 0, '', fontsize=10, color='black',
                      bbox=dict(boxstyle='round,pad=0.3', fc='yellow', ec='red', alpha=0.9),
                      zorder=6, visible=False)
    hover_last = {'idx': -1, 'visible': False}

    # 可见范围内抽稀绘制散点（大文件性能优化）
    def update_scatter():
        x0, x1 = ax.get_xlim()
        mask = (x_full >= x0) & (x_full <= x1)
        n_vis = int(np.count_nonzero(mask))
        if n_vis > MAX_VISIBLE_POINTS:
            # 等间隔抽稀：取 n_vis 中的 MAX_VISIBLE_POINTS 个
            step = (n_vis + MAX_VISIBLE_POINTS - 1) // MAX_VISIBLE_POINTS
            idx = np.flatnonzero(mask)[::step]
        else:
            idx = np.flatnonzero(mask)
        scatter.set_data(x_full[idx], y_full[idx])

    def on_scroll(event):
        if event.inaxes is not ax or event.xdata is None:
            return
        zoom_in = event.button == 'up'
        factor = 1.0 / 1.5 if zoom_in else 1.5
        if event.key in ('control', 'ctrl'):
            # Ctrl+滚轮：缩放 Y 轴，以鼠标 Y 位置为锚点（锚点数据点不动）
            y0, y1 = ax.get_ylim()
            cy = event.ydata
            n0 = cy - (cy - y0) * factor
            n1 = cy + (y1 - cy) * factor
            if n1 - n0 < 1.0:
                return
            ax.set_ylim(n0, n1)
        else:
            # 普通滚轮：缩放 X 轴，以鼠标 X 位置为锚点（锚点数据点不动）
            x0, x1 = ax.get_xlim()
            n0 = event.xdata - (event.xdata - x0) * factor
            n1 = event.xdata + (x1 - event.xdata) * factor
            if n1 - n0 < 1e-6:
                return
            ax.set_xlim(n0, n1)
        update_scatter()
        fig.canvas.draw_idle()

    def on_press(event):
        if event.inaxes is ax and event.button == 1:
            state['press_x'] = event.xdata
            state['press_y'] = event.ydata
            state['press_xlim'] = ax.get_xlim()
            state['press_ylim'] = ax.get_ylim()

    def on_motion(event):
        if event.inaxes is not ax or event.xdata is None:
            return
        if state['press_x'] is not None:
            # 拖拽平移：X/Y 双向跟随鼠标
            x0, x1 = state['press_xlim']
            y0, y1 = state['press_ylim']
            dx = event.xdata - state['press_x']
            dy = event.ydata - state['press_y']
            ax.set_xlim(x0 - dx, x1 - dx)
            ax.set_ylim(y0 - dy, y1 - dy)
            update_scatter()
            fig.canvas.draw_idle()
            return
        # 悬停：找可见范围内距鼠标最近的脉冲点（屏幕像素距离 < 20px 才高亮）
        x0, x1 = ax.get_xlim()
        mask = (x_full >= x0) & (x_full <= x1)
        if not np.any(mask):
            if hover_last['visible']:
                hl_marker.set_visible(False)
                hl_text.set_visible(False)
                fig.canvas.draw_idle()
            hover_last['visible'] = False
            return
        px, py = ax.transData.transform(np.column_stack([x_full[mask], y_full[mask]])).T
        dist = np.hypot(px - event.x, py - event.y)
        k = int(np.argmin(dist))
        if dist[k] <= 20.0:
            idx = int(np.flatnonzero(mask)[k])
            hl_marker.set_data([x_full[idx]], [y_full[idx]])
            hl_text.set_text("测量点 #%d\nt = %.3f ms\nf = %.0f Hz"
                             % (idx + 1, x_full[idx], y_full[idx]))
            hl_text.set_position((x_full[idx] + (x1 - x0) * 0.01,
                                  y_full[idx] + (ax.get_ylim()[1] - ax.get_ylim()[0]) * 0.02))
            hl_marker.set_visible(True)
            hl_text.set_visible(True)
            if hover_last['idx'] != idx or not hover_last['visible']:
                fig.canvas.draw_idle()
            hover_last['idx'] = idx
            hover_last['visible'] = True
        else:
            if hover_last['visible']:
                hl_marker.set_visible(False)
                hl_text.set_visible(False)
                fig.canvas.draw_idle()
            hover_last['visible'] = False

    def on_release(event):
        state['press_x'] = None

    def on_double(event):
        if event.dblclick:
            ax.set_xlim(0.0, meta['activity_duration'] * 1000 * 1.02)
            ax.set_ylim(0, ymax if ymax > 0 else freq.max() * 1.08)
            update_scatter()
            fig.canvas.draw_idle()

    update_scatter()
    fig.canvas.mpl_connect('scroll_event', on_scroll)
    fig.canvas.mpl_connect('button_press_event', on_press)
    fig.canvas.mpl_connect('button_release_event', on_release)
    fig.canvas.mpl_connect('motion_notify_event', on_motion)
    fig.canvas.mpl_connect('button_press_event', on_double)

    print("打开窗口：滚轮缩放X / Ctrl+滚轮缩放Y / 左键双向拖拽平移 / 双击复位")
    print("检测到脉冲 %d，频率测量点 %d，使用 %d 周期平均，活动时长 %.1f ms" % (
        meta['detected_pulses'], meta['pulses'], meta['cycles'],
        meta['activity_duration'] * 1000))
    plt.show()


def main():
    args = [a for a in sys.argv[1:] if not a.startswith('--')]
    opts = {a.split('=', 1)[0]: (a.split('=', 1)[1] if '=' in a else True)
            for a in sys.argv[1:] if a.startswith('--')}

    if len(args) < 1 or '--help' in opts or '-h' in sys.argv[1:]:
        print(__doc__)
        sys.exit(0)

    path = args[0]
    fs = float(args[1]) if len(args) > 1 else 3.125e6
    thresh = int(opts.get('--thresh', 128))
    cycles = int(opts.get('--cycles', 8))
    reverse = '--reverse' in opts
    ymax = float(opts.get('--ymax', 0))
    save_path = opts.get('--save', '')
    draw_d1 = '--d1' in opts
    draw_d2 = '--d2' in opts
    derivative_window_ms = float(opts.get('--dwindow-ms', 5.0))
    d1_window_ms = float(opts.get('--d1-window-ms', derivative_window_ms))
    d2_window_ms = float(opts.get('--d2-window-ms', derivative_window_ms))

    if not os.path.isfile(path):
        print("错误：文件不存在 - %s" % path)
        sys.exit(1)

    if '--selftest' in opts:
        selftest(path, fs, thresh, cycles, reverse)
    else:
        plot_show(path, fs, thresh, ymax, cycles, reverse, save_path or None,
                  draw_d1, draw_d2, d1_window_ms, d2_window_ms)


if __name__ == '__main__':
    main()