简介：光声峰峰值成像利用光吸收产生的超声信号重建组织内部光吸收分布，在生物医学光学成像与病变识别中具有实用价值。面向光声成像研究者、生物医学工程相关专业学生以及需要快速构建成像算法的开发者，这份MATLAB资源提供了一套完整的峰峰值幅值成像代码包。压缩包内共有4个m文件，整体仅3KB，分别对应信号预处理与滤波、主流程控制、三维曲面可视化以及数据翻转平滑等模块，代码精简、分工明确，便于直接阅读和二次修改。已有433人学习下载。通过研读和运行这些脚本，读者能够掌握光声信号去噪、峰峰值特征提取、三维图像重建等关键操作，理解从原始声波信号到光吸收分布图的完整处理链条，还可根据自身数据调整参数以优化成像质量，适合作为光声成像入门与算法验证的参考资料。
1. 光声峰峰值成像：为什么用峰峰值而不是最大值或均值
做光声成像的同学应该遇到过这样一个选择：从采集到的射频（RF）信号重建图像时，到底取信号的哪个特征来代表这个像素？有人直接取最大值，有人计算均值，还有人在做“幅值成像”时下意识地用了信号的峰峰值。区别在于，光声信号本质上是激光激发产生的超声波瞬态波，它通常是双极性振荡，包含着正负极值。直接取最大值只关注了正向半波，漏掉了负向半波的信息，而取均值几乎不能反映信号的瞬态强度。
峰峰值（Peak-to-Peak, PP）定义为信号在某一时间窗口内最大值与最小值之差。在光声峰峰值成像中，这个值表达了该空间位置处光声信号的绝对振荡幅度，对直流偏移不敏感，也天然抑制了缓慢变化的背景漂移。相比最大值、均值或均方根，峰峰值能在不依赖绝对基线的情况下刻画吸收体的浓度和尺寸，而且计算简单、无需拟合。这篇文章就围绕这个主题展开：从信号模型到MATLAB实现，再到参数调优和三维扩展。适合正在写光声信号处理脚本、做图像重建或调试采集系统的工程师阅读，哪怕是刚接触MATLAB光声数据分析的人，也可以照着例子跑通。
2. 光声信号模型与峰峰值特征提取原理
2.1 光声信号的产生与探测
光声成像的原理并不复杂：短脉冲激光照射吸收体（如血管中的血红蛋白），吸收体快速热膨胀产生压力波，这个压力波被超声换能器接收，形成时域信号，我们称之为一个A-line。典型的光声信号类似于带阻尼的振荡，正半周和负半周幅度可能不对称，但整体包络反映了吸收体的光吸收分布。由于换能器带宽有限，信号并不是理想的狄拉克脉冲，而是一个有限频宽的瞬态波。
这种信号的特征参数很多。对于一个采集到的RF数据矩阵，行代表时间采样，列代表空间位置（一维扫描时是一个A-line索引，二维B-scan时是扫描角度或x位置，三维体积数据则是x和y两个方向的排列）。光声峰峰值成像要做的，就是沿着时间轴（每一列）计算信号的峰值到峰值的幅度，再把这个幅度映射到对应的空间坐标上，形成强度图像。
2.2 峰峰值（Peak-to-Peak）的定义与数学表达
设光声信号为 ( s(t) )，其采样序列为 ( s[n], n=1,\dots,N )。对应一个空间点的峰峰值定义为：
[ PP = \max_{n \in W} s[n] - \min_{n \in W} s[n] ]
其中 ( W ) 是所选的时间窗口，可以是整个A-line长度，也可以是激光发射后的一个特定时间区间。窗口的选择直接决定了成像的深度范围：窗口起点对应探测起始深度，窗口长度对应成像深度厚度。
在MATLAB中，这个计算可以用最朴素的循环完成，但更好的做法是使用max和min函数在指定维度上作用。对于二维矩阵data(Nt, Nx)，一行代码就能完成：
pp = max(data, [], 1) - min(data, [], 1);这里的max(data, [], 1)是沿第一个维度（时间维度）求最大值，返回一个1×Nx的行向量。min同理。两个向量相减得到每个A-line的峰峰值。这个操作没有显式循环，效率很高。
不过要注意，直接对原始RF信号做max和min有一个陷阱：信号中的随机噪声会产生额外的尖峰。如果噪声是零均值高斯白噪声，那么取最大和最小会受到噪声极值的影响，且采样点越多，这种影响越大。所以在实际工作中，我更倾向于先对信号进行带通滤波，滤除带外噪声，再计算峰峰值。但滤波本身会引入振铃，可能改变信号的极值，这一点在后面的章节会详细讨论。
2.3 与最大值、均方根等幅值特征对比
如果需要向同事解释为什么要用峰峰值，最直接的办法是把几种常用的幅值特征放在一起对比。下表列出了常见特征的定义、适用场景和局限性：
| 特征 | 定义 | 对直流偏移敏感性 | 抗噪性 | 光声适用性 |
|---|---|---|---|---|
| 最大值 | (\max_{n \in W} s[n]) | 敏感 | 差 | 只反映正向半波，容易受基线漂移影响 |
| 最小值 | (\min_{n \in W} s[n]) | 敏感 | 差 | 只反映负向半波，很少单独使用 |
| 峰峰值 | (\max - \min) | 不敏感 | 中等 | 完整反映振荡幅度，适合双极性信号 |
| 均值 | (\frac{1}{N}\sum s[n]) | 不敏感 | 好 | 正负抵消后接近零，无法表征光声强度 |
| 均方根RMS | (\sqrt{\frac{1}{N}\sum s[n]^2}) | 不敏感 | 中等 | 反映能量，但计算量大且受直流偏置影响 |
| 包络（希尔伯特） | (| \mathcal{H}(s) |) | 不敏感 | 较好 | 反映瞬时幅度，但需要额外变换，计算量较大 |
在光声成像中，信号是一个近似零均值的衰减振荡，所以均值几乎等于零，没有区分度。最大值只捕捉了正向峰值，如果一个吸收体产生的信号正向半波很小而负向半波很大（这通常发生在信号相位反转时），最大值就会错误地低估信号强度。均方根虽然能反映能量，但对直流偏置没有峰峰值那种“自动抵消”的能力——不过光声信号直流偏置通常已被前期去基线处理。因此，峰峰值在这种场景下是最简单、最稳健的幅值特征。另外，很多超声成像里的“幅值”指的是包络检测后的峰值，但那是经过了Hilbert变换或正交解调之后的事情；这里我们讨论的是RF域直接计算的峰峰值，两者概念不同，写代码时别混淆。
3. 用MATLAB实现峰峰值成像的完整流程
3.1 数据准备：光声RF数据的读取与预处理
3.1.1 读入二维或三维数据矩阵
光声数据常见的存储格式有.mat、.h5、.dat或厂商自定义格式。MATLAB处理这些并不复杂。假设你已经把数据读取到一个变量rf_data中，它的大小是[Nt, Nx]或[Nt, Nx, Ny]。作为示例，我们先生成一段模拟光声信号，以便演示完整流程：
% 模拟光声RF数据：一个类似阻尼振荡的信号 fs = 50e6; % 采样率 50 MHz t = (0:999)/fs; % 时间向量，20微秒 f0 = 5e6; % 中心频率 5 MHz t0 = 5e-6; % 信号到达时间 % 双极性衰减振荡 signal = sin(2*pi*f0*(t-t0)) .* exp(-abs(t-t0)/2e-6); % 构造两个空间位置信号，一个强一个弱 rf_data = zeros(1000, 4); rf_data(:,1) = signal * 1.0 + 0.01*randn(1000,1); rf_data(:,2) = signal * 0.5 + 0.01*randn(1000,1); rf_data(:,3) = signal * 0.2 + 0.01*randn(1000,1); rf_data(:,4) = 0.01*randn(1000,1); % 只有噪声的区域这里rf_data有1000个时间采样点，4个空间位置。真实的系统里，每一列就是一个A-line的原始采样数据。注意我们加了高斯噪声，模拟真实采集中的电子噪声。
读取真实数据时，如果用的是.h5，可以用h5read('file.h5', '/dataset')；如果是.dat且知道采样格式，用fread加reshape。这里不展开，但有一个通用原则：先确认数据的采样率和时间轴长度，因为它决定了后续频率滤波的参数。另外，多通道数据可能按通道交织存储，读取后要重排成[通道数, 采样点数, 帧数]，再转置成我们习惯的[采样点数, 通道数]。
3.1.2 去基线与带通滤波的取舍
原始RF信号往往带有直流偏置，这来自换能器放大电路或模拟前端的偏压。直流偏置如果不去除，在计算峰峰值时会被max和min共同抵消掉，因为峰峰值是差值，所以实际上直流偏置对峰峰值没有影响。这是峰峰值的一个天然优势，因此在做峰峰值成像时，我一般不做专门去基线。
带通滤波则是一个需要认真思考的步骤。光声信号的频谱集中在换能器带宽内，典型范围从几百kHz到几十MHz。带外噪声（特别是低频漂移和高频白噪声）会影响max和min的取值。建议先做一个带通滤波，但要用零相位滤波以避免相位失真导致峰值偏移：
% 带通滤波参数 f_low = 1e6; % 低频截止 1 MHz f_high = 15e6; % 高频截止 15 MHz % 设计巴特沃斯滤波器 [b, a] = butter(4, [f_low, f_high]/(fs/2), 'bandpass'); % 零相位滤波 rf_filtered = filtfilt(b, a, rf_data);这里butter的阶数为4，通带范围1~15 MHz。为什么用4阶而不是更高？阶数越高，频带边缘越陡峭，但群延迟越明显，filtfilt可以减轻相位失真，但高阶滤波器对瞬态信号容易产生振铃。光声信号本身是宽带瞬态，过高阶数的滤波器会在信号前后沿处产生明显振荡，反而干扰峰峰值的计算。我的经验是2~4阶巴特沃斯就足够，最常用的反而是2阶。
如果不想滤波，也可以直接计算峰峰值，然后在成像前对峰峰值结果做平滑。但那样噪声尖峰的影响可能会残留在图像中。所以，建议顺序是：滤波 -> 计算峰峰值。注意滤波后信号的正负峰值幅度会略微改变，但不影响相对对比度。
3.2 峰峰值计算的核心操作
3.2.1 沿时间轴提取每个A-line的峰峰值
有了预处理后的rf_filtered，计算峰峰值非常简单。对于二维数据：
pp_image_1d = max(rf_filtered, [], 1) - min(rf_filtered, [], 1);得到的pp_image_1d是一个长度为通道数的向量，每个值就是该A-line的峰峰值。对于三维数据（例如 B-scan 或多个角度采集的体积数据），rf_data大小是[Nt, Nx, Ny]，我们需要沿第一个维度计算，保留另外两个维度：
pp_image_2d = max(rf_filtered, [], 1) - min(rf_filtered, [], 1); % 结果维度为 [1, Nx, Ny]，squeeze后变为 [Nx, Ny] pp_image_2d = squeeze(pp_image_2d);这里的max(x, [], 1)在本版本MATLAB中沿第一维返回1×Nx×Ny数组。如果不写squeeze，后面看图像时会有一个单一维度，不易操作。这是新手常踩的坑。
3.2.2 参数说明：窗口选择、噪声阈值与输出映射
计算峰峰值时，不一定总要用整个时间长度。有时信号只出现在某个深度范围，而其他时间段都是纯噪声。如果整段计算，噪声会导致每个位置都有一个基础峰峰值，降低了对比度。此时应该选择信号所在的时间窗口。最常见的做法是用激光触发延迟和换能器焦距来估计时间窗口。例如，光声信号在 t=5μs 到达，持续约5μs，那么可以只取 t=3μs 到 t=10μs 之间的采样点：
win_start = 3e-6; % 窗口起始时间（秒） win_end = 10e-6; % 窗口结束时间 idx_start = round(win_start * fs) + 1; idx_end = round(win_end * fs); windowed_segment = rf_filtered(idx_start:idx_end, :, :); pp_windowed = squeeze(max(windowed_segment, [], 1) - min(windowed_segment, [], 1));之后可设置一个噪声阈值：把峰峰值低于某个值的像素置为背景值，通常这个阈值取无信号区域的峰峰值平均加若干倍标准差。也可以直接用> noise_level作为掩膜。噪声水平可以通过测量采集一段无激光时的数据得到，或者取每个A-line前几个采样点（还没有接收到光声信号）的峰峰值来估计。
成像映射这一步，对于二维B-scan，直接imagesc(pp_image)；对于一维线扫描，用plot(pp_image)或image逐列填充。但要注意，实际坐标轴映射还需要知道换能器的扫描位置与时间-深度转换关系。基本关系是深度 (d = c \cdot t / 2)，声速 (c) 在生物组织中约为1540 m/s，所以时间轴可以换算为深度轴。如果imagesc时手动指定坐标：
depth = t * 1540 / 2 * 1e3; % 单位 mm imagesc(pp_image_1d); colormap('hot'); xlabel('通道位置');这里只是示意，真实数据还需要根据采集几何定义x轴和y轴。
3.3 用MATLAB把峰峰值结果显示为灰度图
现在我们把上面的步骤整合到一个完整的、可运行的示例中，并添加必要的注释：
% 完整峰峰值成像示例（模拟数据） fs = 50e6; t = (0:999)/fs; % 生成两个不同强度的光声信号，加上噪声 s = sin(2*pi*5e6*(t-5e-6)) .* exp(-abs(t-5e-6)/1.5e-6); rf_data = [s; 0.6*s; 0.3*s; zeros(size(s))]' + 0.02 * randn(1000,4); % 滤波 [b, a] = butter(4, [1e6 15e6]/(fs/2), 'bandpass'); rf_f = filtfilt(b, a, rf_data); % 窗口选择：2~12 us w_start = 2e-6; w_end = 12e-6; [~, i1] = min(abs(t - w_start)); [~, i2] = min(abs(t - w_end)); rf_w = rf_f(i1:i2, :); % 峰峰值 pp = max(rf_w, [], 1) - min(rf_w, [], 1); % 显示 figure; subplot(2,1,1); imagesc(rf_w); title('滤波后RF数据（每列一个通道）'); subplot(2,1,2); plot(pp, 'o-'); xlabel('通道序号'); ylabel('峰峰值'); grid on;代码中[~, i1] = min(abs(t - w_start))是一种常见的按时间找索引方法：求时间向量与目标值差的绝对值，最小值所在位置就是最近的索引。这种写法比round(w_start*fs)+1更直观，尤其在时间向量有非均匀间隔时。运行后可以看到，前三个通道的峰峰值明显高于只有噪声的通道，而且近似呈1:0.6:0.3的关系，符合我们设置信号的幅度比例。这说明峰峰值成像能定量反应光声信号幅值。
如果数据是三维的，只需要把rf_f设为[Nt, Nx, Ny]，然后max和min沿第一维计算，最后squeeze得到二维图像。后续显示时用imagesc即可。
4. 峰峰值成像参数调优：窗口、阈值与对数压缩
4.1 时域窗口长度对峰峰值的影响
上面提到，窗口的选择直接决定了参与计算的信号长度。窗口太短，可能只截取到半个振荡周期，导致max和min无法捕捉完整的正负峰，峰峰值偏小。窗口太长，则会引入更多的噪声采样点，噪声极值变大，峰峰值升高。那么如何定量选择窗口？
最稳妥的方法是基于已知信号的带宽：光声信号的振荡周期约为中心频率的倒数。以5 MHz中心频率为例，一个周期是200 ns。5个周期的持续时间约为1μs。如果系统采集的信号长度有5μs，那么信号振荡大约25个周期。为了捕捉完整的正负峰值，窗口至少需要覆盖信号的主要能量段，通常取2~3倍信号持续时间。但实际中我们总希望窗口尽量窄，以提高轴向分辨率。一个折中做法是：先观察几个A-line的波形，手动估计信号从开始到衰减至10%的时间长度，然后取这个时间窗口。
窗口起点也很关键。如果起点在信号到达之前，那么这段纯噪声会被纳入计算，每个位置的峰峰值都会被噪声抬高。我一般先取整段信号计算一个粗略峰峰值图像，找到信号覆盖的时间区间，然后再用这个区间重新计算。这个过程可以做成半自动：用包络检波（Hilbert变换）找到信号能量超过噪声阈值的区域，然后用该区域作为窗口。
4.2 噪声阈值与动态范围
峰峰值成像的后处理中，阈值设定的意义不仅在于美观，更在于避免低幅值区域被噪声主导。光声信号幅度差异较大，例如强吸收体（血管）与弱吸收体（脂肪）可能相差两个数量级。如果直接显示线性灰度图，低幅值区域的噪声会非常明显。常见的做法是：
- 计算每个A-line在无信号时间段的峰峰值，统计其均值和标准差，得到噪声水平 (N_0)。
- 设阈值为 (T = N_0 + k \cdot \sigma_N)，k 通常取3~5。
- 将小于 (T) 的峰峰值置为 (T) 或直接置为背景值。
在MATLAB中实现如下：
% 假设已经得到峰峰值矩阵 pp_image (MxN) % 估计噪声：取每个A-line最前面100个采样点 noise_seg = rf_data(1:100, :); noise_pp = max(noise_seg, [], 1) - min(noise_seg, [], 1); threshold = mean(noise_pp) + 3 * std(noise_pp); pp_thresholded = pp_image; pp_thresholded(pp_thresholded < threshold) = 0;注意这种阈值处理是对每个A-line独立估计的，可以有效应对通道间噪声水平不一致的情况。对于三维数据，可以先把噪声段改成一个三维子块，然后对时间维计算峰峰值，得到一张二维噪声图，再扩展为整个体积使用的阈值分布。
阈值设得太高，小吸收体会丢失；太低则噪声背景明显。我一般在调试时使用滑块交互式调整：
figure; im = imagesc(pp_image); thr = mean(noise_pp) + 3*std(noise_pp); caxis([thr max(pp_image(:))]);这虽然不改变数据，但通过调节显示范围caxis可以快速评估阈值。后续再将同样的裁剪逻辑应用到最终保存的数据。
4.3 对数压缩与动态范围显示
光声图像往往需要显示很大的灰度范围。峰峰值从几十到几万，直接imagesc会压缩低幅值信息。通常采用对数压缩，将数据映射到人眼可感知的动态范围。MATLAB中可以直接用log(pp_image + 1)或db函数。其中加1是为了避免log(0)产生无穷大。常见的动态范围压缩公式是：
[ I_{\text{display}} = \frac{\log_{10}(1 + \alpha \cdot I)}{\log_{10}(1 + \alpha \cdot I_{\max})} ]
其中 (\alpha) 是压缩系数，控制压缩强度。(\alpha) 越大，低幅值区域被分配更多灰度级。实现为：
alpha = 100; I_disp = log10(1 + alpha * (pp_image / max(pp_image(:)))); imagesc(I_disp);也可以使用商用的显示方案：先去掉底噪，再取对数。很多人会把这个处理直接放进行业成熟的重建脚本里。但要注意，对数压缩会改变图像之间的相对比例，如果后续要做定量分析（如血氧饱和度），应该在原始峰峰值图像上进行，而不是压缩后的显示数据。
4.4 性能优化：向量化与并行计算
光声数据集往往很大，例如三维体积数据有1000个时间点 × 256 × 256 的空间点，数据量约6550万个采样点，用双层循环计算峰峰值会非常慢。上面的max(..., [], 1)方法已经向量化，但对三维数据需要小心处理。一个常见的错误是直接使用max(rf, [], 1)，结果维度变成[1, Nx, Ny]，然后错误地使用max(result, [], 2)等，导致计算混乱。正确做法见前面代码。
如果体积数据太大，内存无法同时容纳，可以采用分块处理。例如沿沿x方向每次读取一块256×16×Nt的子体积，计算峰峰值后再拼接。MATLAB的tall数组和imageDatastore也能处理，但更直接的还是for循环配合matfile对象。此外，如果有Parallel Computing Toolbox，可以对通道维使用parfor：
% 将数据按通道分给不同并行worker pp_all = zeros(Nx, Ny); parfor i = 1:Nx pp_all(i, :) = max(rf_time(:, i, :), [], 1) - min(rf_time(:, i, :), [], 1); end注意这里rf_time是一个[Nt, Nx, Ny]的三维数组，parfor切片时MATLAB会自动传输每个i需要的数据。不过并行开销大，只有当数据维度较长时使用才划算。对于单次二维图像，向量化操作已经足够快，不需要并行。
5. 验证与进阶：三维峰峰值投影、多波长比值与滤波顺序技巧
5.1 三维体积数据上的峰峰值投影
当你有三维光声数据（B-scan 堆叠成体积）时，除了生成每一层的二维图像，往往还需要做一次最大强度投影（MIP）或平均强度投影展示整体结构。峰峰值成像可以在三维上直接计算体积数据，得到一个三维幅值矩阵，然后再沿深度方向做投影。假设rf_volume维度是[Nt, Nx, Ny]，先求峰峰值体数据：
pp_volume = max(rf_volume, [], 1) - min(rf_volume, [], 1); pp_volume = squeeze(pp_volume); % 变成 Nx×Ny这样每个(x,y)位置都是一个深度上最大振荡幅度的投影。但这个投影丢失了深度信息。更常用的是最大强度投影沿深度轴（时间轴）做，而不是先算峰峰值。两种做法目的不同：如果你想要的是“这个位置不管多深，都存在强吸收体”，应该对原始RF信号取深度方向的峰峰值；如果你想要的是“某个深度切面”，应该对每个时间片分别做峰峰值。前者把三维体积压缩成二维地图，适合观察整体血管形态；后者用于逐层分析。
MATLAB中实现最大强度投影的方式是：
mip_image = max(pp_volume, [], 3); % 如果 pp_volume 是 Nx×Ny×Nz如果你的原始三维数据维度是[Nt, Nx, Ny]，那么pp_volume恰好是[Nx, Ny]，无法再沿深度投影。这种情况下，你需要把体积按时间窗口分割成多个子体积，每个子体积计算峰峰值，得到多幅二维图像，然后再沿深度轴（即子体积序号）取最大。这个过程中，窗口的划分会影响投影结果的轴向分辨率。小窗口带来高分辨率，但噪声变大；大窗口则相反。实践中我通常将时间窗口设为信号振荡周期的3~5倍。
5.2 多波长光声中的峰峰值比值
光声成像的一个重要应用是多波长成像，通过不同波长下吸收体的吸收系数差异来计算血氧饱和度等参数。峰峰值成像在这种场景下相当方便：因为峰峰值与光声信号幅值成正比，而信号幅值又正比于吸收系数和光通量。如果我们采集了两个波长 (\lambda_1) 和 (\lambda_2) 的峰峰值图像 (PP_1) 和 (PP_2)，那么比值 (R = PP_1 / PP_2) 可以抵消掉光通量不均和表面衰减的影响，保留与血氧饱和度相关的吸收比信息。
计算时要注意分母不能为零，而且要使用同一时间窗口。MATLAB代码可以写成：
PP1 = max(data_lambda1, [], 1) - min(data_lambda1, [], 1); PP2 = max(data_lambda2, [], 1) - min(data_lambda2, [], 1); ratio = (PP1 + eps) ./ (PP2 + eps);其中eps是为了防止分母为零。比值图像一般不需要再做大动态范围压缩，直接显示即可。不过多波长成像更精细的做法是先用上述峰峰值方法得到强度图像，再在空间上做一个中值滤波平滑掉噪声，然后再算比值。经过平滑之后，比值图会稳定很多。
5.3 一个被忽略的细节：滤波顺序对峰峰值的影响
最后分享一个我踩过的坑。很多人会先对RF信号做带通滤波，再计算峰峰值。这个顺序很自然，但滤波器的振铃效应可能让信号产生额外的正负极值，特别是在信号边缘处。filtfilt虽说是零相位，但零相位不代表零振铃。对于短促的光声脉冲，一个通带截止频率接近中心频率的滤波器，会在脉冲前后形成对称的过冲，这个过冲幅度有时甚至比原始信号峰值还大。
解决的办法是：先计算峰峰值，然后再对峰峰值图像进行空间域平滑。或者，如果一定要在时间域滤波，那就选用 Butterworth 2阶且通带带宽足够宽的滤波器，并且在计算峰峰值之前截去信号两端各几微秒的数据，避开振铃的起始段。我个人的策略是：原始RF信号先只做去除直流（减均值），用max和min计算峰峰值；然后对峰峰值图像做3×3中值滤波和低通空间滤波。这样既保留了峰峰值对双极性振荡的真实度量，又避免了滤波振铃的干扰。在实际的大光声数据集上，这个策略对比度稳定、伪影少。这也是我在处理光声峰峰值成像时最常推荐的做法。
本文还有配套的精品资源，点击获取
