ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

RF超声时间序列读取与MATLAB成像处理:从原始数据到B超图像

RF超声时间序列读取与MATLAB成像处理:从原始数据到B超图像 简介RFRadio Frequency数据是超声成像中最原始的信号形态记录着声波在组织界面反射后的完整信息其质量直接决定图像分辨率与诊断价值。这份资源针对超声成像学习者与信号处理研究者提供读取与预处理 RF 数据的 MATLAB 脚本覆盖时间序列分析、滤波去噪、FFT 时频转换、动态范围压缩等关键环节帮助用户从原始信号到可视化超声图像快速上手。压缩包内含 1 个 m 文件体积仅 2KB轻量实用。已有 148 人学习下载。通过该脚本读者可以掌握 RF 数据加载、预处理及时间序列可视化的完整思路结合描述中的方法可进一步探索声速估计、多普勒血流分析与组织弹性成像等高级应用适合作为超声信号处理入门与科研验证的参考工具脚本结构简洁便于二次修改与功能扩展。1. RF 超声时间序列为什么 ReadRFdata.m 是读懂超声原始信号的第一关RF 超声时间序列是超声成像里最接近原始物理量的数据形态。探头接收到的回声电压没有经过任何显示处理直接落盘就是它这份 ReadRFdata.zip 里的 ReadRFdata.m 就是专门把这类数据读进 MATLAB 的脚本。很多拿到 .raw 文件的人第一反应是 imshow结果看出一片噪声其实不是数据坏了而是缺了从时间序列到图像的那几步。它适合做超声成像课题、医学信号处理作业的工程师和研究生。别指望一个脚本能直接输出 B 超图它只是把数据从黑匣子里取出来后面的滤波、包络检波、动态范围压缩还得自己接。2. RF 数据与显示图像的差距从探头电压到 MATLAB 矩阵的完整链路2.1 RF 信号、IQ 数据和 B 超图像三者差在了相位做超声成像的人常把 RF、IQ、B 模式图像混着叫但三者的信息量差一个层级。RF 是探头收到的原始电压序列中心频率通常在 3.5 MHz 到 10 MHz 之间采样率按奈奎斯特条件一般取中心频率的 4 到 8 倍。IQ 数据是 RF 经过正交解调后的基带复数形式采样率可以降下来但相位信息完整保留。B 超图像则是把 RF 或 IQ 做包络检测、对数压缩、扫描转换之后的结果相位信息已经丢光了。所以当你手里只有 B 超图时多普勒血流方向、组织弹性这类需要相位信息的分析基本没法做。RF 时间序列的价值就在这它在固定扫描位置连续采集多帧每一帧都是一个采样点 × 扫描线的矩阵第三维是时间。ReadRFdata.m 这类脚本的实质工作就是把这个三维矩阵从二进制文件里还原出来。这里要强调一个常见误区不要把 RF 时间序列和普通的生理时间序列混为一谈。心电、脑电是一维标量随采样时刻变化而 RF 时间序列是「每个采样深度、每条扫描线都各自有一条随帧号变化的时间曲线」。后面做 LSTM 输入时这个结构直接决定你怎么 reshape。2.2 读取 RF 文件前的四个假设尺寸、字节序、交错方式和头部长度在跑脚本之前必须先确认文件的原始存储约定。常见做法是先用一个小工具把文件头 256 字节 dump 成十六进制看有没有 ASCII 关键字段没有的话就得按采集设备的参数反推。四个关键假设是每帧由 samplesPerLine × numLines 个采样点组成采样点一般是 int16 或 float32数据可能按帧连续存放也可能按扫描线交错存放文件开头可能有自定义头长度从几十字节到几千字节不等。ReadRFdata.m 能帮上忙的前提是你已经知道前两个值。一般会先看采集软件导出的配置或者用下面这段 MATLAB 试探文件长度和采样点的关系% 先用 dir 拿到字节数再反推帧数 info dir(rf_data.raw); bytes info.bytes; samplesPerLine 4096; % 每条线采样点由 AD 时钟和成像深度决定 numLines 128; % B 模式一帧的扫描线数 headerBytes 0; % 未知时先置 0有余数再改 frameBytes samplesPerLine * numLines * 2; % int16 每样本 2 字节 numFrames floor((bytes - headerBytes) / frameBytes);这段代码的逻辑是从总字节数里减去可能的头部再除以一帧的字节数得到合法帧数。参数说明里headerBytes 需要你从文件头分析得到置 0 时如果算出的 numFrames 带小数、或者乘回总字节数对不上就说明文件确实有自定义头或者 samplesPerLine 不是你以为的值。算出来的 numFrames 如果能与采集时设置的帧数互相印证脚本的读法才可靠。2.3 从单帧到时间序列第三维的物理含义很多超声设备支持「高帧率序列采集」。心脏这类运动组织帧率通常是 30 到 100 fps剪切波弹性成像则要求数千 fps。第三维的物理含义就是慢时间slow time它对应组织随时间的运动而第一维是快时间fast time对应声波在深度方向往返的时间差。读入之后rf(: , :, 1)是第一帧、rf(:, 64, 3)是第三帧第 64 条扫描线。想要观察某个深度点的波形随帧号变化直接squeeze(rf(depthIndex, lineIndex, :))。但要提醒一句慢时间的采样率是帧率不是 AD 采样率。做多普勒 FFT 时用的窗口是沿第三维切出来的千万不要和第一维混在一起。3. 把 RF 时间序列读进 MATLABReadRFdata.m 的 reshape 逻辑与参数校验3.1 三个核心参数采样率、扫描线数、成像深度怎么联动读取脚本能不能一次跑对取决于填入的三个参数采样率 fs、扫描线数 numLines、每条扫描线采样点数 samplesPerLine。它们之间的关系是参数典型取值错误后果fs 采样率40 MHz取 20 MHz 时奈奎斯特边界不到 10 MHz5 MHz 载波混叠fc 中心频率5 MHz带通滤波器中心设错会滤掉有效信号samplesPerLine4096reshape 尺寸错包络图整片斜纹numLines128一帧宽度不对图像比例失真numFrames100时间轴长度不对多帧数据错位samplesPerLine 的估算公式是depth 0.10; % 成像深度 10 cm c 1540; % 软组织声速 m/s fs 40e6; % AD 采样率 Hz samplesPerLine ceil(2 * depth / c * fs);这个公式来自声波往返时间。深度 10 cm、声速 1540 m/s 时往返时间约 130 微秒按 40 MHz 采样就是约 5200 个点。实际操作中设备往往按最大深度固定一个长度常见的 4096、5120 都可能出现。如果脚本里写死 4096而你用的是 5120读进来第一帧就会串入好几条线。3.2 一次标准的读取fopen、fread、reshape 三步ReadRFdata.m 这类脚本最核心的部分可以收敛成一个函数function rf ReadRFdata(filename, samplesPerLine, numLines, numFrames, headerBytes) % 读取原始 RF 二进制文件返回 (samplesPerLine, numLines, numFrames) fid fopen(filename, rb); if fid -1 error(打不开文件: %s, filename); end fseek(fid, headerBytes, bof); % 跳过自定义文件头 rf fread(fid, [samplesPerLine, numLines * numFrames], *int16); fclose(fid); rf reshape(rf, samplesPerLine, numLines, numFrames); end这里的关键点在fread的尺寸是用[samplesPerLine, numLines * numFrames]表示的。文件连续存放时第一条扫描线的采样点连续排列第二条紧接着开始MATLAB 的 fread 按列填矩阵正好一列是一条扫描线。先把二维读出来再用reshape恢复成三维顺序是「列优先」也就是第一维采样点、第二维扫描线、第三维帧。如果采集设备是 float32把*int16改成*float32但帧字节数也要从 2 改成 4。另外建议不要在脚本里直接覆盖原始变量。读取后立刻执行rf_raw rf;这样后面处理坏了还有后悔药可吃。MATLAB 里被覆盖的变量可没有撤销键临时文件一旦被清理就得重新解压读取。3.3 快速校验一帧包络图能看出 reshape 对不对读对了不代表万事大吉还得校验一次。最直接的办法是取第一帧沿深度方向做包络检测再看图像是否符合超声解剖结构frameEnv abs(hilbert(double(rf(:, :, 1)), [], 1)); imagesc(20*log10(frameEnv / max(frameEnv(:)) eps)); axis image; colormap gray; colorbar;RF 是 int16 时直接传给 hilbert 会丢失精度而且直流偏置会影响包络所以先double(rf(:, :, 1))。hilbert(..., [], 1)里的第三个参数表示沿第一维操作也就是沿深度方向做解析信号变换。如果 reshape 正确这张图应该有明显的深浅层次组织边界是连续曲线如果看到斜向条纹、整块错位多半是 samplesPerLine 或 numLines 填反了。没有 Signal Processing Toolbox 时可以退而求其次做滑动 RMS 近似包络但相位精度比不上希尔伯特后面做多普勒就别省这个工具箱。4. 从滤波到灰度映射RF 时间序列变成超声图像的四个处理环节4.1 带通滤波去掉带外噪声再谈信号原始 RF 信号里除了组织回声还有探头振铃、电路热噪声和低频漂移。直接做包络检测噪声会被一起放大图像看起来脏。标准做法是先做一个以中心频率 fc 为中心的带通滤波器带宽一般取 fc 的 40% 到 80%fs 40e6; % AD 采样率 fc 5e6; % 探头中心频率 fbw 0.6 * fc; % 带宽 3 MHz [b, a] butter(4, [(fc - fbw/2) / (fs/2), (fc fbw/2) / (fs/2)], bandpass); rfLine double(rf(:, 64, 1)); rfLineF filtfilt(b, a, rfLine); % 零相位滤波避免时间延迟butter(4, ...)设计的是 4 阶巴特沃斯滤波器括号里的截止频率必须归一化到奈奎斯特频率 fs/2 的 0 到 1 区间。计算带外噪声的幅度时用filtfilt比filter稳。filter有相位延迟会在包络峰值位置上偏移几十个采样点filtfilt双向滤波消除相位偏差代价是计算量翻倍处理大矩阵时要注意内存。如果发现滤波后的信号幅度整体变小先检查归一化截止频率是不是写错了。比如 fc10 MHz、fs20 MHz 时上限已经顶到奈奎斯特频率带通滤波器会变得很尖甚至数值不稳定。4.2 包络检测为什么要用希尔伯特而不是直接取绝对值RF 信号是高频载波乘以组织反射系数直接看波形正负交替非常快。取绝对值只能得到半个周期的脉冲仍然带着载波振荡放到图像上就是条纹状伪影。希尔伯特变换构造解析信号后取模得到的是瞬时幅度包络这才是 B 超需要的低频幅度信息envLine abs(hilbert(rfLineF)); plot((0:length(envLine)-1)/fs*1e6, envLine); xlabel(时间/us); ylabel(包络幅度);包络峰值对应组织界面的强反射宽度由脉冲长度决定而不是由载波周期决定。想看中心频率是否漂移可以在滤波后做一次 FFTf (0:1023) / 1024 * fs; plot(f, abs(fft(rfLineF, 1024))); xline(fc, --);这里补零到 1024 点只是为了画图平滑不是提高频率分辨率。如果频谱峰值明显偏离 fc可能探头中心频率和标称值不一致也可能滤波器的中心频率设错了。4.3 对数压缩与灰度映射动态范围参数怎么调组织回声强度跨度很大最强反射和最弱散射之间可以相差 60 到 80 dB线性灰度早就饱和了。所以要把包络幅度转成对数尺度再做灰度映射env abs(hilbert(double(rf(:, :, 1)), [], 1)); dB 20*log10(env / max(env(:)) eps); dB(dB max(dB(:)) - 45) max(dB(:)) - 45; % 45 dB 显示窗 img uint8(255 * (dB - min(dB(:))) / (max(dB(:)) - min(dB(:)))); imshow(img);这段代码里20*log10 把幅度比换算成 dBeps为了防止取对数时出现零。45 这个参数是显示动态范围它决定你保留多少弱的散射信号。调到 60图像里能看到大量斑点噪声但微弱组织的边界也更清楚调到 30图像干净但浅表弱回声可能直接消失。实际调试时我会先设 50再根据浅表区域和深部区域的对比度微调 5 到 10 dB。5. 避坑指南RF 时间序列读取与成像的五个常见翻车点5.1 reshape 出斜条纹帧交错方式搞反了现象按前面脚本读出来单帧包络图是一堆斜向条纹组织边界是断裂的甚至图上能看到明显接缝。原因fread 按物理文件顺序填矩阵但采集设备有时不是「帧内逐线连续」存储而是把多帧数据按「帧间交错」存放。先存每一帧的第一条扫描线再存第二条。这样直接用 [samplesPerLine, numLines * numFrames] 读每一列里混了不同帧的数据。解决先看采集软件有没有「interleaved」这个选项有的话读入后需要按帧做一次 permute。常见做法是先读成二维矩阵再重排三维顺序tmp fread(fid, [samplesPerLine, numLines * numFrames], *int16); rf reshape(tmp, samplesPerLine, numLines, numFrames); % 若出现斜纹尝试换顺序 rf permute(reshape(tmp, [numFrames, numLines, samplesPerLine]), [3 2 1]);不要死记某一种顺序要以第 3.3 节的包络图能否出现连续光滑边界为准。每换一种组合就重画一次哪个维度组合能出正常结构就用哪一种。5.2 满屏噪声或全 0字节序和采样位数不对现象读出来的矩阵最大值只有个位数或者全是 ±1 的假脉冲包络图看不出任何组织层次。原因文件可能是 16 位小端格式脚本却用了*int8也可能是大端字节序fread 默认按本机小端解析于是每个采样点的两个字节被交换幅值变成随机数。解决明确指定字节序rf fread(fid, [samplesPerLine, numLines * numFrames], *int16, ieee-le);不确定时读入前先看文件前 8 个字节如果十六进制是5831 7A2B这种明显高低字节交换过的样子就把ieee-le改成ieee-be。换一次字节序包络峰值位置和幅度都会恢复正常。5.3 包络图横向拉丝希尔伯特沿错了维现象包络图看起来是一行一行横向条纹深度方向反而没有纹理边界全部模糊。原因hilbert(frame)不指定维度时默认沿第一个非单一维度操作。RF 矩阵第一维是深度采样点第二维是扫描线如果你写成hilbert(double(rf(:, :, 1)))MATLAB 沿第二维做解析信号变换等于把相邻扫描线的相位信息混进了包络。解决强制沿深度维envFrame abs(hilbert(double(rf(:, :, 1)), [], 1));多帧循环里同样要写全第三个参数。这一点对时间序列尤其重要因为每帧都要做一次包络循环体里很容易因为图省事漏掉维度参数结果整批帧都错。5.4 文件长度算不清自定义头和数据尾现象numFrames 算出来带小数或者读到最后一帧时 fread 报错说文件已经到末尾。原因采集端在前面写了配置头尾部还附带了时间戳或 CRC 校验字段。只按数据长度整除当然对不上。解决先用 dir 拿总字节数再尝试减去不同的 headerBytes直到整除bytesNoHeader info.bytes - headerBytes; framesExact bytesNoHeader / (samplesPerLine * numLines * 2);如果 framesExact 不是整数就拿floor的结果做帧数并接受最后一帧之后可能还有几个字节的尾部数据。更稳的做法是让采集软件直接导出带明文头的格式比如 Verasonics 的 .bin 就有完整配置头头解析不对时先用十六进制工具确认头长度而不是反复试 numFrames。5.5 处理完的临时变量把原始数据覆盖了现象调试到一半发现 rf 已经被滤波后的数据覆盖想重新看原始波形只能从头解压。原因脚本里写了rf filtfilt(b, a, rf)这种原地操作MATLAB 工作区又没有备份。解决读取后立刻留底rf_raw rf; % 备份原始数据之后的滤波、包络、压缩全部写到新变量名里比如rf_filt、envSeq。如果内存吃紧也可以只保存原始文件路径和读取参数需要原始数据时重跑一次读取函数这比重下整个压缩包快得多。6. 进阶验证用合成 RF 信号测试脚本再谈 LSTM 时间序列预测的入口6.1 没有实测数据时先造一个已知答案的 RF 序列新拿到一个读取脚本最怕的是连数据本身是不是坏的都不知道。我的习惯是先合成一条 RF 信号把峰值位置算出来再用脚本读真实的 .raw对照包络峰值是否符合同样的物理规律。fs 40e6; fc 5e6; N 4096; t (0:N-1) / fs; centerDelay 10e-6; rfSynth cos(2*pi*fc*t) .* exp(-((t - centerDelay).^2) / (2*(0.5e-6)^2)); envSynth abs(hilbert(rfSynth)); [~, idx] max(envSynth); fprintf(包络峰值在 %.2f us应约为 %.2f us\n, idx/fs*1e6, centerDelay*1e6);高斯包络的中心时刻就是回波延时合成信号跑通后再拿真实 RF 数据走同样的包络流程。如果真实数据的第一条线看不出明确的脉冲峰值就要回头检查 5.1 到 5.4 的任何一个环节。这个方法不挑设备所有 RF 数据都能用。6.2 RF 时间序列特征怎么给 LSTM 用ReadRFdata.m 读出来的是三维矩阵但 LSTM 的输入通常是序列长度 × 特征维度的二维矩阵。常见做法是把每帧的扫描线维度压缩成一个特征向量取固定一条线或者取所有线的平均包络。seqData zeros(numFrames, samplesPerLine); for k 1:numFrames envFrame abs(hilbert(double(rf(:, :, k)), [], 1)); seqData(k, :) mean(envFrame, 2); % 每帧一条深度包络特征 end这样组织出来的 seqData每一行是一帧列是深度位置直接可以接sequenceInputLayer和lstmLayer做时间序列预测或异常检测。注意一点平均包络已经丢掉了相位信息做血流多普勒、组织弹性这类任务不能用这个特征应该改用 IQ 数据或复数 RF 的相位项。从那以后我每次拿到新的 RF 数据集都会先跑一遍第 6.1 节里的合成对照再读真实文件确认 reshape 和包络峰值都对得上才开始下面的处理。希望帮到你。本文还有配套的精品资源点击获取
返回列表