- 实验目的
1、掌握血管壁搏动位移的超声检测方法;
2、初步掌握信号相似性的评价方法;
3、初步掌握信号滤波的方法。
- 实验内容
1、简述基于超声射频信号的管壁搏动位移及脉搏波速超声测量过程。
2、 基于数据sim_bmd6_45,使用归一化相关系数为管壁追踪的代价函数,测量13线空间位置的管壁搏动位移曲线。
3、 对第2题中获得的管壁搏动位移曲线进行滤波,消除台阶状的量化误差。
4、 使用第3题中已滤波的13线管壁搏动位移曲线估计脉搏波速(PWV)。
必要的实验参数:
1)13线空间位置平均分布在38mm的颈动脉上。
2)每线空间位置上,有1333个射频信号。
3)1333个射频信号的采样时间为1秒。
4)射频信号的采样频率为100MHz。
- 方法步骤
1、超声探头对颈部血管位置扫描时发射13个声束(如图1所示),这时我们取其中一个位置来分析。
图 1
图 2
图2中显示的是第一个位置(声束)的成像,可以看到这个声束中不同扫描深度处的成像。由于扫描时对每一个声束进行了1333次超声扫描,所以可以将1333次扫描看作是时间轴。我们用第一次扫描结果中的血管壁位置(我选的是左边的血管壁)与其它 1332次扫描结果进行相关系数测量,然后找到相关系数最大的位置,用这个位置减起始位置得到像素位移(以像素个数为单位),再将像素位移转换为实际位移(以mm为单位),得到的结果即为这个位置的管壁搏动位移。同理,计算得到13个位置的管壁搏动位移,观察它们相同幅度在时间上的延时dt,然后再计算得到13个测量位置dx,计算dx/dt即可计算出脉搏波速pwv。
1. 管壁搏动位移曲线测量
- 1参数选择
up=3101;
down=3320;
sampling_width=down-up;%3100~3320
number_k=1333;
pixed_width=1540/10^8*1000/2;
这里我选取的测量位置是图2中左边血管壁的位置,在3100—3320处。sampling_width是采样宽度,即血管壁宽度,pixed_width是每个像素实际表示多少mm的宽度,这里除2是因为超声成像是回声成像。
- 2数据读取
for f=1:13
url="line_noise"+num2str(f);
load(url,'-mat')
data=RFline_noise;
这里用循环读入13组数据。
- 3位移
for m=1:number_k
center=data(up:down,1);
search=data(up-300:down+300,m);
for i=1:300*2-sampling_width
ncc(i)=corr2(center,search(i:i+sampling_width,:));
end
l=300-find(ncc==max(ncc))+1;
distanse(m,f)=l*pixed_width;
end
这里外循环用于计算13个位置中每一个位置的位移,内循环用于在一个位置中计算第一次扫描结果与剩下的1332次扫描结果的相关系数。具体实现步骤是:以第一次扫描结果为中心块,中心块左右300的位置为搜索块,在搜索块中找到与中心块最匹配的位置,然后计算像素位移再转换为实际位移distanse,这里用300减最匹配位置是因为脉搏跳动的时候血管壁要向up一侧扩张。
- 4数据保存
由于上述过程使用了三个循环嵌套,所以shij2复杂度非常高,为了后续计算速度能更快,所以要保存上述的数据。
- 5绘图
x_axis=(0:1332).*ones(13,1);
x_axis=x_axis';
y_range=(0:12)*38/12;
y_axis=y_range.*ones(1333,1);
figure
plot3(x_axis,y_axis,distanse,'Color',[0 0 0]),hold on,grid on
ylabel('Position(mm)'),xlabel('Time(ms)'),zlabel('Distension(mm)')
在这一部分我计算了13个位置的x坐标,即时间坐标,y坐标。然后用黑色的线条绘制了13条管壁搏动位移曲线。
2. 滤波、消除台阶状的量化误差
wp=6/1333;
ws=10/1333;%设置阻带频率
Rp=0.5;%设置通带波纹系数
Rs=1.5;%设置阻带波纹系数
这里我设计了低通滤波器对13个位置的脉搏位移曲线进行了滤波,这里通带频率为6Hz,阻带频率为10Hz,通带波纹系数为0.5,阻带衰减系数为1.5。其中对通带频率和阻带频率进行了归一化处理。
[N,Wn]=buttord(wp,ws,Rp,Rs);
[b,a]=butter(N,Wn);
filter_datas(:,f)=filter(b,a,distanse(:,f));
使用巴特沃斯滤波器对每一条曲线进行滤波。
3. 估计脉搏波速(PWV)
s0=0.06;
for i=1:13
dt0=find(filter_datas(:,i)>s0);
dts=dt0(1);
dt(i)=dts;
end
dt(2:13)=dt(2:13)-dt(1);
dt(1)=0;
ds=0:38/12:38;
cftool(dt,ds)
这里我找了大于0.06的第一个位置,然后在13条曲线中找到这个位置对应的时间即为dt,对dt和ds进行数据拟合,斜率即为脉搏波速PWV。
- 实验结果
对第一个位置进行管壁搏动位移曲线测量得到的曲线如图3所示,可以看到搏动最大位移为0.4mm,曲线存在阶梯状的量化误差。
图 3
图 4
通过计算得到13条管壁搏动位移曲线如图4所示,图中的曲线存在阶梯状的量化误差,为了消除这个误差我对曲线进行了滤波处理。
图 5
通过滤波得到13条管壁搏动位移曲线如图5所示,对比曲线的趋势可以初步判断这是中年组的管壁搏动位移曲线,所以测量得到的PWV应该为6.0m/s左右。
图 6
拟合ds和dt得到的结果如图6所示,观察拟合曲线的斜率可以得到PWV为6.378m/s,这个结论和上面的预期一致,是中年组的PWV。
完整代码:
%%
clear
clc
up=3101;
down=3320;
sampling_width=down-up;%3100~3320
number_k=1333;
wp=6/1333;
ws=10/1333;
Rp=0.5;
Rs=1.5;
pixed_width=1540/10^8*1000/2;
for f=1:13
url="line_noise"+num2str(f);
load(url,'-mat')
data=RFline_noise;
for m=1:number_k
center=data(up:down,1);
search=data(up-300:down+300,m);
for i=1:300*2-sampling_width
ncc(i)=corr2(center,search(i:i+sampling_width,:));
end
l=300-find(ncc==max(ncc))+1;
distanse(m,f)=l*pixed_width;
end
[N,Wn]=buttord(wp,ws,Rp,Rs);
[b,a]=butter(N,Wn);
filter_datas(:,f)=filter(b,a,distanse(:,f));
end
save distanse
save filter_datas
save distanse
%%
close all
load filter_datas.mat
load distanse.mat
x_axis=(0:1332).*ones(13,1);
x_axis=x_axis';
y_range=(0:12)*38/12;
y_axis=y_range.*ones(1333,1);
figure
plot3(x_axis,y_axis,distanse,'Color',[0 0 0]),hold on,grid on
ylabel('Position(mm)'),xlabel('Time(ms)'),zlabel('Distension(mm)')
figure
plot3(x_axis,y_axis,filter_datas,'Color',[0 0 0]),hold on,grid on
ylabel('Position(mm)'),xlabel('Time(ms)'),zlabel('Distension(mm)')
s0=0.06;
for i=1:13
dt0=find(filter_datas(:,i)>s0);
dts=dt0(1);
dt(i)=dts;
end
dt(2:13)=dt(2:13)-dt(1);
dt(1)=0;
ds=0:38/12:38;
cftool(dt,ds)
%view([0 0])%正视图
%%