一种超声频谱自适应噪声估计包络提取方法
【专利摘要】一种超声频谱自适应噪声估计包络提取方法，步骤为：预存储经过处理的t时刻声功率谱线数据Pt(i)；读取数据Pt(i)，计算其正向谱最大值标签和负向谱最大值标签，同时分类确定数据Pt(i)的噪声序列分布；统计计算数据Pt(i)噪声分布序列的平均噪声；计算数据Pt(i)的正向平均功率、负向平均功率、全局平均功率，通过数据筛选获得数据Pt(i)的最佳平均噪声；按时间顺序连续计算三根相邻谱线Pt?2(i)、Pt?1(i)、Pt(i)的最佳噪声值，对三组最佳噪声进行数据处理获得数据Pt(i)的最终平均噪声。本发明方法受信噪比影响较小，有利于正负最大频谱的计算从而精确估算血流的流速。
【专利说明】
一种超声频谱自适应噪声估计包络提取方法
技术领域
[0001 ]本发明涉及到超声技术，尤其涉及一种超声频谱自适应噪声估计包络提取方法。
【背景技术】
[0002] 超声技术被广泛应用于医学成像和测量。在脉冲波多普勒系统测量血流速度的过 程中，系统按照一定时间间隔重复发射脉冲信号，在连续两次发射之间的某一时刻接收散 射回波，再通过解调、滤波后，多普勒信号频谱从几兆赫兹搬移到以零频为中心，带宽几千 赫兹的音频范围内。对解调出的音频信号进行两路处理，一路声音播放，通过扬声器可以听 到多普勒信号的声音，一路进行傅里叶变化显示为以时间为横轴，频率为纵轴的功率谱图。 在临床上对血管疾病进行诊断时，经常会截取多普勒信号声谱图中的某一段进行分析，从 中提取一些参数，比如平均流速，最大流速，最小流速，收缩期、舒张末期最大流速比，阻尼 指数，脉动指数等。所有这些参数都是基于声谱图包络曲线或者平均频率曲线计算获得，因 此准确计算声谱图的包络信息在临床上有重大的意义。
[0003] 传统的声谱图包络提取方法包括百分比法、过阈值法、改进的阈值法、混合法、几 何法、改进的几何法、自适应阈值法等等。专利US5，287，753、US5，935，074采用自适应估计 阈值的方法，实际操作中为了增加最大频率估算的鲁棒性而设置较大的阈值，而阈值来自 于频谱的信噪比，当信噪比较小时，估算不是很准确；专利CN 1875887A采用改进的几何法 估算最大频率并利用之前时刻的噪声计算当前谱线的噪声，这种方法对混叠状态的处理以 及信噪比较低时包络的处理不尽人意。以上的方法基本都存在噪声计算不准的情况，另外 通过预定噪声阈值的方法来确定实时的噪声鲁棒性较低。
[0004] 图1为现有技术中的多普勒信号处理流程，其包络提取在处理压缩和增益之前进 行，因此在图像冻结状态不能进行包络的调节。
【发明内容】
[0005] 针对现有技术中的不足，本发明提出了一种超声频谱自适应噪声估计包络提取方 法，用于实时计算谱线噪声，准确提取声谱图包络，从而可提供能准确计算临床应用参数的 数据。
[0006] -种超声频谱自适应噪声估计包络提取方法，其特征在于，包括以下步骤：
[0007] 1 )、预存储经过解调、滤波、DFT处理的t时刻声功率谱线数据Pt(i)，同时进入步骤 2)与步骤4);
[0008] 2)、读取预存储的t时刻声功率谱线数据Pt(i)，计算其正向谱最大值标签和反向 谱最大值标签，同时分类确定t时刻声功率谱线数据P t(i)的噪声序列分布，其中分类确定 噪声序列的计算公式如下：
[0009]
[0010]其中，P表示正向谱最大值标签；
[0011] η表示负向谱最大值标签；
[0012] MaxIndex表示最大谱频率标签；
[0013] i表示t时亥Ij谱线上的点序号；
[0014] I表示逻辑"或"关系；
[0015] 0:η_1 表不0、1、2......n_l 序列；
[0016] p+1 :Max 表不 p+1、p+2、p+3......MaxIndex序列；
[0017] p+1 :n_p 表示 p+1、p+2、p+3......n_p 序列；
[0018] 噪声序列计算完毕，进入步骤3);
[0019] 3 )、统计计算t时刻声功率谱线数据Pt ( i )噪声分布序列Ptnciise3l ( i )的平均噪声 Ptnoise，其中平均噪声Ptnoisd十算步骤如下：
[0020] 31)把噪声序列?_^1(1)按照各噪声?_^1(1)的值的大小，按从小到大的顺序排 列，然后把排好序的Ptr^s^a)划分为N等分；
[0021 ] 32)对每一等分内的噪声数据求平均值，记为mean( j)，其中(X j<N;
[0022] 33)对噪声序列mean( j)去掉一个最大值、去掉一个最小值后，获得新序列meannew (j)，对meannew( j)序列求其平均值，记为meanall;
[0023] 34)新序列meannew(j)与阈值meanall进行比较判断获得t时刻声功率谱线数据Pt (i)的平均噪声MeanNoiset，具体公式如下：
[0024]
[0025]其中，N表示噪声功率Ptnciisel(i)序列的等分数；
[0026] k表示在meannew( j)序列中第一个大于阈值meanall的数据meanall (k)的标签值， 0^k<N-2；
[0027] 之后进入步骤4);
[0028] 4)、计算t时刻声功率谱线数据Pt( i)的正向平均功率Ppt(i)，负向平均功率Pnt(i)， 全局平均功率Pmt( i)，对ppt( i)、Pnt( i)、Pmt( i)、MeanNoiset进行数据筛选获得t时刻声功率 谱线数据的最佳平均噪声noiset，其中所述数据筛选具体的计算公式如下：
[0029] noiset=Min{MeanNoiset，Max(PPt(i)，Pnt(i)，Pmt(i))}
[0030] 其中MinU表示取最小值函数；
[0031] MaxO表示取最大值函数；
[0032] 5)、重复步骤1)到4)，按时间顺序连续计算卜2、卜1、丨三时刻声功率谱线数据卩*一2 (i)、Pt-i(i)、Pt(i)的最佳平均噪声，分别记为noiset-2，noiset-i，noiset，对三组噪声进行时 间方向数据处理获得t时刻声功率谱线数据的最终平均噪声，记为NoiseEnd t，所述时间方 向数据处理的具体计算公式如下：
[0033]
[0034] 其中：t-Ι表示t时刻的前一时刻；
[0035] t_2表示t-Ι时刻的前一时刻；
[0036] Pt-i(i)表示t-1时刻声功率谱线数据；
[0037] Pt-2( i)表示t-2时刻声功率谱线数据；
[0038] MaxO表示取最大值函数；
[0039] 6)、对当前谱线Pt(i)减去最终计算噪声并积分，以基线为边界重新确定新的正向 谱最大值标签和负向谱最大值标签，具体计算公式如下：
[0040]
[0041]
[0042] 其中，base表示谱基线位置，可取0、1、2……M-1;
[0043] f c表示量化后的壁滤波截止频率；
[0044] M表示量化后的谱频率；
[0045] k的取值范围为0、1、2……M-I;
[0046] ;表示在base+fc到M-I之间取最大值函数；
[0047] 表示在0到base-fc之间取最小值函数；
[0048] find(x，y)表示已知y，查找X的查询函数；
[0049] fp表示正向谱最大值标签；
[0050] f η表示负向谱最大值标签；
[00511 7)、按时间顺序依次连接所有时刻谱线的fp以及fn值，即可获得PWD谱的正负包络 值，存储正负包络值完成包络的提取。
[0052]本发明所提供的自适应噪声估计包络提取方法，无需外部预测任何参数作为辅助 参数，在信噪比较低时也能准确计算谱线的噪声同时准确提取谱图包络，此外在图像冻结 和非冻结状态下均进行包络提取。
【附图说明】：
[0053]图1为现有技术中的多普勒信号处理方法流程图；
[0054]图2为本发明自适应噪声估计包络提取方法的流程图；
[0055]图3为多普勒声功率谱线统计计算平均噪声的流程图；
[0056]图4为本发明方法处理过程中不同时刻谱图显示以及参数示意图；
[0057]图5为所选取时刻的多普勒声功率谱线示意图。
【具体实施方式】：
[0058]本发明对PW(脉冲多普勒)和CW(连续多普勒)成像均适用，下面以PW多普勒成像为 例，对本发明自适应噪声估计包络提取方法进行详细说明。
[0059]本发明的处理过程为:预先缓存PW谱数据到计算机硬盘，再从硬盘读取不同时刻 的PW谱数据，接着，先计算初始正、负向最大频率，分类获得当前谱线的噪声序列，对噪声序 列数据统计、筛选、滤波处理获得当前谱线较为准确的最终平均噪声，原始的谱线减去最终 平均噪声后重新计算新的更准确的正、负向最大频率;连接所有正、负向最大频率即可获得 PW频谱包络。
[0060]图2为本发明自适应噪声估计包络提取方法的流程图，具体步骤如下：
[0061 ] 1)预存储经过解调、滤波、离散傅里叶变换(DFT)处理的声功率谱数据Pt (i)到硬 盘空间并以文件形式保存，声功率谱数据Pt(i)按照时间顺序不断更新以保证数据的实时 性，之后，同时进入步骤2)与步骤4);
[0062] 2)读取硬盘数据文件获得当前时刻t的声功率谱线数据Pt (i )，对当前时刻t的声 功率谱线数据Pt(i)预计算初始正向谱最大值标签和反向谱最大值标签，同时分类确定当 前谱线Pt(i)的噪声序列分布，其中，预计算初始正向谱最大值标签和反向谱最大值标签可 采用过阈值法和改进几何法等技术，此部分内容为现有技术，在此不赘述;分类确定噪声序 列的计算公式如下：
[0063]
[0064] 其中如图4所示：
[0065] p表示正向谱最大值标签；
[0066] η表示负向谱最大值标签；
[0067] MaxIndex表示最大谱频率标签，本实例中MaxIndex为255;
[0068] 另如图5所示：
[0069] i表示t时刻谱线上的点序号，i取值范围在0~255;
[0070] 噪声序列计算完毕后，进入步骤3)对噪声数据进行统计处理；
[0071] 3)统计计算t时刻声功率谱线数据Pt(i)噪声分布序列的平均噪声，记为P tncilse3l,如 图3所示，具体包括如下步骤：
[0072] 31)把噪声序列?_^1(1)按照各噪声?_^1(1)的值的大小，按从小到大的顺序排 列，然后把排好序的Ptr^s^a)划分为N等分；
[0073] 32)对每一等分内的噪声数据求平均值，记为mean( j)，其中0< j<N;
[0074] 33)对噪声序列mean( j)去掉一个最大值、去掉一个最小值后，获得新序列meannew (j)，对meannew( j)序列求其平均值，记为meanall;
[0075] 34)新序列meannew(j)与阈值meanall进行比较判断获得t时刻声功率谱线数据Pt (i)的平均噪声MeanNoiset，具体公式如下：
[0076]
[0077]其中，N表示噪声功率Ptncilsel(i)序列的等分数；
[0078] k表不在meannew( j)序列中第一个大于阈值meanall的数据meanall (k)的标签值， 0^k<N-2；
[0079] 之后进入步骤4);
[0080] 4)如图4所示，按正、负方向计算当前时刻t谱线的正向平均功率Ppt(i)，负向平均 功率P nt (i )，全局平均功率Pmt (i )，求平均功率为现有技术，在此不再赘述，在计算完上述三 个平均功率后，对数据进行筛选以获得当前t时刻谱线pt(i)的最佳平均噪声，其中最佳平 均噪声具体的计算公式如下：
[0081]
[0082]
[0083]其中
表示对z函数在X到y范围内的积分函数；
[0084] f c可设为0到31之间整数值；
[0085] base取值范围为0到255之间；
[0086] MinU表示对中括号中的两个值进行比较并取其最小值；
[0087] MaxO表示对小括号中的两个值进行比较并取其最大值；
[0088] 5)重复步骤1)到4)，按时间顺序连续计算t-2、t_l、t三时刻声功率谱线数据Pt-2 (i)、Pt-i(i)、Pt(i)的最佳平均噪声，分别记为noiset-2，noiset-i，noiset，对三组噪声进行时 间方向数据处理获得t时刻声功率谱线数据的最终平均噪声，记为NoiseEnd t，所述时间方 向数据处理的具体计算公式如下：
[0089]
[0090] 其中：t-Ι表示t时刻的前一时刻；
[0091] t-2表示t-Ι时刻的前一时刻；
[0092] Pt-i(i)表示t-Ι时刻声功率谱线数据；
[0093] Pt-2(i)表示t-2时刻声功率谱线数据；
[0094] MaxO表示取最大值函数；
[0095] 6)对当前谱线Pt(i)减去最终平均噪声NoiseEndt，并积分，以基线为边界重新确定 新的正向谱最大值坐标和负向谱最大值坐标，具体计算公式如下：
[0096]
[0097]
[0098] 其中，base表示谱基线位置，可取0、1、2……255;
[0099] f c表示量化后的高通滤波器截止频率；
[0100] M表示量化后的谱频率；
[0101] k的取值范围为0、1、2……255;
[0102] 表示在base+f c至_-1之间取最大值函数；
[0103] g示在0到base-fc之间取最小值函数；
[0104] find(x，y)表示已知y，查找X的查询函数；
[0105] f p表示正向谱最大值标签；
[0106] f η表示负向谱最大值标签；
[0107] 7)按时间顺序依次连接所有谱线的f ρ值，以及按时间顺序依次连接所有谱线的fn 值，即可获得PWD谱的正、负包络值，存储正、负包络值，完成包络的提取。
[0108]本发明的上述实施例仅是为清楚地说明本发明所作的举例，而并非是对本发明的 实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其 它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而这些属于本发 明的精神所引申出的显而易见的变化或变动仍处于本发明的保护范围之中。
【主权项】
1. 一种超声频谱自适应噪声估计包络提取方法，其特征在于，包括W下步骤： 1) 、预存储经过解调、滤波、DFT处理的t时刻声功率谱线数据Pt(i)，同时进入步骤2)与 步骤4); 2) 、读取预存储的t时刻声功率谱线数据Pt(i)，计算其正向谱最大值标签和反向谱最大 值标签，同时分类确定t时刻声功率谱线数据Pt(i)的噪声序列分布，其中分类确定噪声序 列的计算公式如下：其中，P表示正向谱最大值标签； η表示负向谱最大值标签； Maxindex表示最大谱频率标签； i表示t时刻谱线上的点序号； 陵示逻辑"或'关系； 0:n-l表示0、1、2……n-1序列； P+1 :Max表不 P+1、p~i~2、p+3......Maxindex 序列； P+1 :n-p 表示 p+l、p+2、p+3......n_p 序列； 噪声序列计算完毕，进入步骤3); 3 )、统计计算t时刻声功率谱线数据Pt( i )噪声分布序列Ptnoisel( i )的平均噪声Ptnoise，其 中平均噪声Ptnnise计算步骤如下： 31) 把噪声序列Ptnnisel(i)按照各噪声PtnniselQ)的值的大小，按从小到大的顺序排列， 然后把排好序的PtnDisel( i )划分为N等分； 32) 对每一等分内的噪声数据求平均值，记为mean(j)，其中0《j<N; 33) 对噪声序列mean( j)去掉一个最大值、去掉一个最小值后，获得新序列meannew( j)， 对meannew( j)序列求其平均值，记为meanall; 34) 新序列meannew(j)与阔值meanall进行比较判断获得t时刻声功率谱线数据Pt(i)的 平均噪声MeanNoiset，具体公式如下：其中，N表示噪声功率Ptnnisel ( i )序列的等分数； k表示在meannew (j)序列中第一个大于阔值meanal 1的数据meanal K k)的标签值，0《k <N-2； 之后进入步骤4); 4)、计算t时刻声功率谱线数据Pt(i)的正向平均功率Ppt (i)，负向平均功率Pnt (i)，全局 平均功率Pmt( i)，对Ppt( i)、Pnt( i)、Pmt( i)、MeanNoiset进行数据筛选获得t时刻声功率谱线 数据的最佳平均噪声noiset，其中所述数据筛选具体的计算公式如下： noiset=Min{MeanNoiset，Max(Ppt(i) ,Pnt(i) ,Pmt(i))} 其中Min{}表示取最小值函数； Max()表示取最大值函数； 5) 、重复步骤1)到4)，按时间顺序连续计算t-2、t-1、tS时刻声功率谱线数据Pt-2(i)、 Pt-i(i)、Pt(i)的最佳平均噪声，分别记为noiset-2，noiset-i，noiset，对Ξ组噪声进行时间方 向数据处理获得t时刻声功率谱线数据的最终平均噪声，记为NoiseEndt，所述时间方向数 据处理的具体计算公式如下：其中：t-1表不t时刻的前一时刻； t-2表不t-1时刻的前一时刻； Pt-i(i)表示t-1时刻声功率谱线数据； Pt-2( i)表示t-2时刻声功率谱线数据； Max()表示取最大值函数； 6) 、对当前谱线Pt(i)减去最终计算噪声并积分，W基线为边界重新确定新的正向谱最 大值标签和负向谱最大值标签，具体计算公式如下：其中，base表示谱基线位置，可取0、1、2……M-1; fc表示量化后的壁滤波截止频率； Μ表示量化后的谱频率； k的取值范围为0、1、2……M-1;表示在base+fc到M-1之间取最大值函数；良示在0菌Jbase-fc之间取最小值函数； find(x，y)表示已知y，查找X的查询函数； fp表示正向谱最大值标签； f η表示负向谱最大值标签； 7) 、按时间顺序依次连接所有时刻谱线的fpW及fn值，即可获得PWD谱的正负包络值，存 储正负包络值完成包络的提取。
【文档编号】A61B8/00GK106073822SQ201610377416
【公开日】2016年11月9日
【申请日】2016年5月31日
【发明人】丁勇, 张晓明, 肖均文
【申请人】青岛惠尔医疗科技有限公司