库库插值是良方！（高频流数据版）
最近写了不少程序，记录一个软件上面的设计，目前是将解析好的数据通过 USB发送给上位机，其中数据打包的思路是把结构化的多通道数据转换为一维字节流，并加入“帧头和帧尾”以确保上位机能够准确断帧和对齐频道。 ## **固定帧长数据包（Packit）** 现在的数据属于连续的高频流数据，如果直接发送 ASCII 字符串（如使用 `printf`），会因为字符转换效率低、数据量翻倍而导致卡顿。因此，**必须发送原始二进制字节（Binary Data）**。 现在设计的数据包（共 38 个字节）： |字节偏移|字段名称|数据类型|说明| | --- | --- | --- | --- | |0 ~ 1|帧头 (Header)|uint8_t|固定为 0xAA, 0x55（用于上位机对齐查找包头）| |2|包计数器 (Counter)|uint8_t|0 ~ 255 循环，用于上位机检测是否丢包| |3 ~ 34|通道数据 (Payload)|int32_t × 8|8 个通道的 32 位脑电数据（8 × 4 = 32 字节）| |35|校验和 (Checksum)|uint8_t|前 35 字节的累加和，用于检测数据传输错误| |36 ~ 37|帧尾 (Tail)|uint8_t|固定为 0x0D, 0x0A (\r\n，符合多数串口软件换行习惯)| 利用 C 语言结构体的内存连续性，可以直接打包： ``` // 强制编译器不进行内存对齐优化，确保结构体大小严格等于 38 字节 #pragma pack(1) typedef struct { uint8_t header[2]; // 帧头 0xAA, 0x55 uint8_t packet_cnt; // 包计数器 int32_t channels[8]; // 8通道数据 (32字节) uint8_t checksum; // 校验和 uint8_t tail[2]; // 帧尾 0x0D, 0x0A } EEG_Packet_t; #pragma pack() // 全局变量 EEG_Packet_t TxPacket; uint8_t global_packet_counter = 0; ``` 这段是数据的容器，要装在这里，接下来看数据的装载： ``` void Pack_And_Send_EEG(int32_t *parsed_channels) { // 装载帧头和帧尾 TxPacket.header[0] = 0xAA; TxPacket.header[1] = 0x55; TxPacket.tail[0] = 0x0D; TxPacket.tail[1] = 0x0A; // 装载计数器并自增 TxPacket.packet_cnt = global_packet_counter++; // 复制解析好的 8 通道 int32_t 数据 (共 32 字节) memcpy(TxPacket.channels, parsed_channels, sizeof(TxPacket.channels)); // 计算校验和 (从 header 到 channels 的所有字节累加) uint8_t sum = 0; uint8_t *ptr = (uint8_t *)&TxPacket; for(int i = 0; i 丢包补偿只能恢复时间位置，不能恢复真实信号，插值数据必须始终带有“无效/估计”标记，不能伪装成真实采样值。 ## **包计数器**  ### **跨越255时仍然成立** 例如上一包为254，当前包为1： 因此丢失2包，即255和0。 Python中可以写成： ``` distance = (current_seq - previous_seq) & 0xFF lost_packets = distance - 1 ``` 这里的 `& 0xFF` 相当于对256取模。 上面写的是： > 每当单片机成功组装并发送完一个数据包，计数器加1。 这句话不够严谨，更严谨的说法是： > 每次采集系统生成一个应当发送的数据包时，就为它分配唯一的序号。 例如： ``` packet.seq = packet_counter++; enqueue_tx_packet(&packet); ``` 如果 MCU 内部发生以下问题：ADC采集完成，但发送队列已经满了；DMA发送任务来不及处理；某个数据包在 MCU 内部被覆盖；BLE发送接口返回忙，程序直接放弃该包；如果计数器只在“实际发送成功”后递增，上位机看到的序号仍然连续，反而无法发现 MCU 内部已经丢失过采样数据。 ### **8位序号存在歧义** 8位序号最多只能表示256种状态。如果一次断开时间太长，丢失超过255包，上位机无法知道究竟循环了多少次。 > 例如上一包为5，当前包仍为5，可能表示：收到重复包；恰好丢失了255包；设备重新启动；丢失了511包，有条件使用最低16位序号；更推荐32位 `sample_index`。 32位计数器即使采样率为1 kSPS，也可以连续运行约： 天 如果每个数据包包含多个采样点，持续时间还会更长。 ## **丢包为什么会影响时间轴** 假设采样率为： 采样周期： 真实时间应该由采样序号决定： 假设采样序列是： ``` 100、101、102、105、106 ``` 说明103、104缺失；如果绘图程序只是把收到的数据依次放在相邻位置： ``` 显示位置：0、1、2、3、4 真实位置：100、101、102、105、106 ``` 那么时间轴被压缩了8 ms。后续所有数据都会向前错位，因此真正需要恢复的是时间位置： ``` 100、101、102、缺失、缺失、105、106 ``` 这里不一定必须填写一个数值，但必须保留两个缺失位置。 ## **对FFT影响的一点讨论** ### **情况一：丢失位置被删除** 如果丢失采样点直接从序列里消失，后面的时间轴被压缩。此时数据已经不再是严格均匀采样： 而普通FFT默认： 因此直接进行FFT，频率和相位分析都会产生误差。 ### **情况二：缺失位置补0** 补0可以恢复样本数量和时间位置，但在信号中引入了突然下降到0的边沿。这会产生宽带频谱泄漏，并激励IIR滤波器的瞬态响应。 ### **情况三：前值保持** 前值保持通常不会产生立刻跳到0的巨大阶跃，但会形成一段水平信号，并在真实数据恢复时产生跳变。 ### **情况四：线性插值** 短缺口下通常比较平滑，但会人为生成一段近似直线，相当于局部低通处理，高频成分会被压低。 ## **补偿方法** 既然可以找到哪里缺了数据，那就可以补回来，这里讨论三个方案 ### **前值保持** 公式为： 适合实时显示；缺口极短；信号变化相对缓慢；不能等待后一个样本到达；缺点是会产生平台；恢复时可能产生阶跃；高频信号误差较大；不能作为真实测量数据。 ### **线性插值** 若两端真实值为 和 ，中间缺失 个点，则第 个插值为： 例如缺失两个点： 则插值分别为： 适合：短缺口；平滑显示；低频趋势恢复；离线处理或允许等待下一包的实时显示；但插值值是估算，不是真实测量；这个缺口越长，可信度越低；肌电等高频、随机性较强的信号不适合长区间插值；可能人为减少功率谱中的高频能量；可能制造不存在的斜坡，插值结果最好保留浮点数，避免额外量化，只有在最后输出到整数格式时才转换。 ### **NaN加质量标记** 这是保存和科研分析中最可靠的基础方法： ``` data[missing_position] = np.nan valid[missing_position] = False ``` 明确： > 这里没有测量数据。 有了 NaN，不能简单把NaN直接送进IIR滤波器，如果把 NaN 直接送入普通 IIR 滤波器： ``` y = filter.process(np.nan) ``` 由于IIR滤波器有反馈状态，NaN可能进入内部状态，之后所有输出都变成NaN。 IIR滤波器一般形式： 一旦历史项中出现NaN，后续计算会持续受到污染。 ## **可以采用分级补偿策略** 不应该对所有长度的缺口都使用同一种插值，假设缺失时间为： 可以根据任务设置策略，如： |缺口|实时显示|原始保存|频谱/特征分析| | --- | --- | --- | --- | |1～2点|线性插值|NaN+无效标记|可酌情使用，但保留标记| |3～10点|插值但标色|NaN+无效标记|通常排除该区域| |较长缺口|显示断线|NaN+无效标记|分段分析| |连接中断|新建数据段|记录事件|不跨段滤波/FFT| ## **详细模拟** 我们使用Python 模拟插值，以及研究缺失数据在频域会有什么影响，展示的时候使用分辨率高的 FFT。  ### **仿真设置** ``` 采样点数 : 92,160 缺失样本数/比例 : 792 / 0.8594% 记录长度理论频率间隔 : 1/T = 0.005556 Hz Hann 窗 ENBW : 1.5/T = 0.008333 Hz 零填充 FFT 栅格间距 : fs/NFFT = 0.000488 Hz 注意：更细的 FFT 栅格不是更强的真实分辨能力；真实分辨率主要靠延长 T。 10 Hz附近主要谱峰： Reference : [(10.123046875, -0.005393100133749951), (10.1728515625, -5.18754105913096), (10.10986328125, -31.439126135762386), (10.13623046875, -31.494161173770728)] Zero fill : [(10.123046875, -0.06699910014441023), (10.1728515625, -5.216832638381734), (10.10986328125, -31.412882691863118), (10.13623046875, -32.08808645213892)] ZOH : [(10.123046875, -0.03559907459472774), (10.1728515625, -5.205575543568598), (10.10986328125, -31.199437020777808), (10.13623046875, -31.503423386907578)] Linear : [(10.123046875, -0.03330143431458907), (10.1728515625, -5.193685957757536), (10.10986328125, -31.188468711646827), (10.13623046875, -31.69564251215458)] Deleted : [(10.126953125, -5.586233149229062), (10.18359375, -6.546209263501059), (10.154296875, -14.058907252136022), (10.24755859375, -14.16637288270421)] 采样点数 : 92,160 缺失样本数/比例 : 792 / 0.8594% 记录长度理论频率间隔 : 1/T = 0.005556 Hz Hann 窗 ENBW : 1.5/T = 0.008333 Hz 零填充 FFT 栅格间距 : fs/NFFT = 0.000488 Hz 注意：更细的 FFT 栅格不是更强的真实分辨能力；真实分辨率主要靠延长 T。 10 Hz附近主要谱峰： Reference : [(10.123046875, -0.005393100133749951), (10.1728515625, -5.18754105913096), (10.10986328125, -31.439126135762386), (10.13623046875, -31.494161173770728)] Zero fill : [(10.123046875, -0.06699910014441023), (10.1728515625, -5.216832638381734), (10.10986328125, -31.412882691863118), (10.13623046875, -32.08808645213892)] ZOH : [(10.123046875, -0.03559907459472774), (10.1728515625, -5.205575543568598), (10.10986328125, -31.199437020777808), (10.13623046875, -31.503423386907578)] Linear : [(10.123046875, -0.03330143431458907), (10.1728515625, -5.193685957757536), (10.10986328125, -31.188468711646827), (10.13623046875, -31.69564251215458)] Deleted : [(10.126953125, -5.586233149229062), (10.18359375, -6.546209263501059), (10.154296875, -14.058907252136022), (10.24755859375, -14.16637288270421)] ``` 采用： 记录时间： 实际采样点数： 信号包含： 其中前两个信号只相差： > 丢包按照每包8个采样点产生，包括：随机整包丢失；连续3包丢失；连续10包丢失；连续25包丢失。 最终缺失： 点 丢失率： > 比较五种情况：完整参考数据；缺失位置补0；前值保持；线性插值；直接删除缺失点。  ### **“高分辨率”FFT** 程序使用： 因此FFT频率栅格间距为： > 这个数很小，但它主要来自零填充，只代表频谱曲线计算得更密。 真正由采集时间决定的基本频率间隔为： 使用Hann窗后，其等效噪声带宽约为： 因此： > `0.000488 Hz`：频谱显示栅格； `0.00556 Hz`：记录长度对应的基本分辨尺度； `0.00833 Hz`：Hann窗的ENBW。 不能把大点数零填充误认为真正提高了频率分辨率。真正提高分辨率，主要依靠延长采集时间。  ## **两个近邻单音能否分开？** 完整信号检测到： |设置信号|FFT检测谱峰| | --- | --- | |10.123 Hz|10.12305 Hz| |10.173 Hz|10.17285 Hz| 频率误差只有FFT栅格量级，两个相隔0.05 Hz的信号能够清楚分开，这说明180秒记录足以分辨它们，因为： 如果要分辨相隔0.001 Hz的两个信号，仅仅把FFT改成几百万点不够。记录时间至少应达到约： 实际考虑窗函数，通常还需要更长。 ## **缺失数据为什么产生频谱泄漏？** 把完整信号写成 ，有效性掩码写成 ： 数据有效数据缺失 补零后的信号实际是： 时域相乘对应频域卷积： 如果没有丢包， 恒等于1，它的频谱主要集中在直流位置，不会明显扩散原信号频谱。 > 发生丢包后，掩码中出现很多缺口， 会包含宽带频率成分。原本集中的单音谱线与它卷积后，会出现：主谱峰幅度下降；谱峰附近出现裙边；宽带噪声底升高；产生不存在的伪频率分量；连续突发缺口产生类似 sinc 的频谱结构。 因此，即使丢失比例不到1%，高动态范围频谱也能清楚看到影响。 ### **各种处理方法的时域误差** 缺失区域RMSE结果： |方法|缺失区域RMSE|全记录RMSE| | --- | --- | --- | |补零|0.8172|0.07576| |前值保持|0.7470|0.06925| |线性插值|0.5280|0.04895| 在本次信号条件下，线性插值最好，因为它利用了缺口前后的两个真实值。 但这不意味着线性插值能恢复真实高频信号。它只是在两个端点之间构造直线： 其中 为缺失点数；缺口越长，真实信号越可能在缺口内发生多次振荡，直线插值与真实值的差别就越大。 ### **补零的频域影响** 补零相当于直接使用缺失掩码： 主要问题是：信号突然跳到0；缺口结束时又突然跳回真实值；两个边沿产生宽带频谱；单音幅度被略微压低；噪声底升高。 本次仿真中，10.123 Hz幅度从： 下降到： 大约下降： 虽然主谱峰变化不大，但在参考信号较弱、动态范围较高时，新增的宽带泄漏会遮盖微弱分量。 ## **前值保持的频域影响** 前值保持为： 它不会像补零一样立即跳到0，但会形成一段水平线;缺口结束时，保持值与当前真实值之间仍可能产生阶跃。对低频或缓慢变化信号比较适合；能用于实时处理，因为不需要等待未来数据；对高频信号产生幅度和相位误差；长缺口会形成明显平台；缺口恢复边沿会带来高频泄漏。 > 仿真中，前值保持比补零的时域RMSE稍低，但频谱误差仍很明显。 ### **线性插值的频域影响** 线性插值比补零和前值保持平滑，因此通常减少突变导致的宽带泄漏，仿真中，10.123 Hz幅度变为： 相对完整参考值的误差约： 10.173 Hz幅度几乎保持不变： 但是线性插值本质上具有低通倾向：缺失区被替换成直线；高频细节被消除；高频功率通常被低估；相位信息不再真实；长缺口可能生成不存在的慢斜坡。线性插值适合：实时波形显示；极短缺口；低频趋势观察；不要求缺失区真实频域特征的场景。 ### **直接删除最危险** 原始数据长度为： 删除缺失点后变为： 如果仍按512 Hz采样率解释，程序会认为记录长度为： 但真实过程持续了180秒。 相当于把180秒发生的信号压缩到了178.45秒，平均时间压缩比例为： 所以频率会整体向上偏移。仿真中10.123 Hz附近的主峰移动到约： 同时出现：10.1543 Hz；10.1836 Hz；10.2476 Hz；等多个伪峰或严重变形的峰。  这里没有严格按照平均压缩比例整体平移，是因为缺失点不是均匀删除，而是在不同时间发生随机和突发删除，所以时间轴是分段跳变，而不是均匀缩放。 结论是： > 绝对不要删除缺失点后，再假定剩余样本仍然是等间隔采样。 ## **小结** 上位机算力大，可以任意插值，另外高点数只能让FFT失真看得更细，但不能把插值数据变回真实测量数据。
