CN105997148B - 脉冲多普勒超高谱分辨率成像处理方法及处理*** - Google Patents

脉冲多普勒超高谱分辨率成像处理方法及处理*** Download PDF

Info

Publication number
CN105997148B
CN105997148B CN201610356317.8A CN201610356317A CN105997148B CN 105997148 B CN105997148 B CN 105997148B CN 201610356317 A CN201610356317 A CN 201610356317A CN 105997148 B CN105997148 B CN 105997148B
Authority
CN
China
Prior art keywords
matrix
freq
velocity
signal
distribution matrix
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201610356317.8A
Other languages
English (en)
Other versions
CN105997148A (zh
Inventor
郭建军
陈惠人
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
FEIYINUO TECHNOLOGY (SUZHOU) CO LTD
Feiyinuo Technology Co ltd
Original Assignee
Vinno Technology Suzhou Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Vinno Technology Suzhou Co Ltd filed Critical Vinno Technology Suzhou Co Ltd
Priority to CN201610356317.8A priority Critical patent/CN105997148B/zh
Publication of CN105997148A publication Critical patent/CN105997148A/zh
Priority to PCT/CN2017/073774 priority patent/WO2017202068A1/zh
Priority to EP17801924.6A priority patent/EP3466343B1/en
Priority to US16/304,268 priority patent/US11197656B2/en
Application granted granted Critical
Publication of CN105997148B publication Critical patent/CN105997148B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/48—Diagnostic techniques
    • A61B8/488—Diagnostic techniques involving Doppler signals
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/52—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/5207—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of raw data to produce diagnostic data, e.g. for generating an image
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/06—Measuring blood flow
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S15/00—Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems
    • G01S15/88—Sonar systems specially adapted for specific applications
    • G01S15/89—Sonar systems specially adapted for specific applications for mapping or imaging
    • G01S15/8906—Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques
    • G01S15/8977—Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques using special techniques for image reconstruction, e.g. FFT, geometrical transformations, spatial deconvolution, time deconvolution
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • G01S7/52017—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00 particularly adapted to short-range imaging
    • G01S7/52046—Techniques for image enhancement involving transmitter or receiver
    • G01S7/52047—Techniques for image enhancement involving transmitter or receiver for elimination of side lobes or of grating lobes; for increasing resolving power
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • G01S7/52017—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00 particularly adapted to short-range imaging
    • G01S7/52053—Display arrangements
    • G01S7/52057—Cathode ray tube displays
    • G01S7/5206—Two-dimensional coordinated display of distance and direction; B-scan display
    • G01S7/52066—Time-position or time-motion displays
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10132—Ultrasound image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20048—Transform domain processing
    • G06T2207/20056—Discrete and fast Fourier transform, [DFT, FFT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Radiology & Medical Imaging (AREA)
  • General Physics & Mathematics (AREA)
  • Biomedical Technology (AREA)
  • Veterinary Medicine (AREA)
  • Biophysics (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • Pathology (AREA)
  • Public Health (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Acoustics & Sound (AREA)
  • Quality & Reliability (AREA)
  • Theoretical Computer Science (AREA)
  • Hematology (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

本发明提供一种脉冲多普勒超高谱分辨率成像处理方法及处理***，所述方法包括：在每根扫查线上均获取N个采样点分别对应的IQ信号，并在快时间方向上，对其进行壁滤波处理，以组成每根扫查线对应的IQ信号壁滤波序列；分别对IQ信号壁滤波序列中的每个采样点依次在慢时间方向和快时间方向上做FFT变换，以获得每个采样点在不同频移上的能量分布矩阵；获取与所述能量分布矩阵匹配的速度分布矩阵；调取用于显示的速度序列，以所述速度序列查询所述速度分布矩阵和与所述速度分布矩阵匹配的能量分布矩阵，获取所述速度序列对应的能量序列，以用于最终的频谱显示。本发明在发射信号带宽较宽或者散射子速度较快时，速度分辨率都得到了很大保证，极大提高了速度谱分辨率。

Description

脉冲多普勒超高谱分辨率成像处理方法及处理***
技术领域
本发明属于医疗超声技术领域，主要涉及一种脉冲多普勒超高谱分辨率成像处理方法及处理***。
背景技术
超声成像因其无创性、实时性、操作方便、价格便宜等诸多优势，使其成为临床上应用最为广泛的辅助诊断的手段之一。其中，脉冲多普勒成像能方便快捷地测出血流的具体流速，在临床诊断中成为某些病症的判断标准。
传统的脉冲多普勒处理技术，在慢时间方向对回波信号做一维傅里叶变换，得到频移分布和其能量，根据频移分布和发射信号的中心频率，计算相应的速度，传统的脉冲多普勒处理技术根据多普勒效应，散射子运动导致的频移大小正比于发射信号频率和散射子运动速度的乘积；如此，传统的脉冲多普勒处理技术由于速度谱的不均一性，速度快的散射子频移带宽较宽，频谱分辨率较差，而速度慢的散射子频移带宽较窄，相应的频谱分辨率较好；同时，多普勒频谱分辨率受发射信号带宽的影响很大，当采用门较小时，发射信号较短，其带宽较宽，多普勒的频谱分辨率较差；另外，不同速度的具有同样反射能量的散射子，快速散射子较早走出取样容积，从而在某段时间内(N个脉冲重复时间(PRT))，快速散射子的反射信号较少，而传统的将采样门内的信号简单平均或者差值法获取采样信号的方法，就使得高速血流的能量被减弱较多，从而最终速度谱的能量不均一。
发明内容
本发明的目的在于提供一种脉冲多普勒超高谱分辨率成像处理方法及处理***。
为了实现上述发明目的之一，本发明一实施方式的脉冲多普勒超高谱分辨率成像处理方法，所述方法包括以下步骤：
S1、在每根扫查线上均获取N个采样点分别对应的IQ信号，以组成每根扫查线对应的IQ信号序列；
在快时间方向上，分别对每根扫查线对应的IQ信号序列进行壁滤波处理，以组成每根扫查线对应的IQ信号壁滤波序列；
S2、分别对IQ信号壁滤波序列中的每个采样点依次在慢时间方向和快时间方向上做FFT变换，以获得每个采样点在不同频移上的能量分布矩阵；
S3、根据所述能量分布矩阵中各个采样点的频移坐标值、初始发射信号的重复扫查频率、初始发射信号的中心频率，获取每个采样点在不同频移上与所述能量分布矩阵匹配的速度分布矩阵；
S4、调取用于显示的速度序列，以所述速度序列查询所述速度分布矩阵和与所述速度分布矩阵匹配的能量分布矩阵，获取所述速度序列对应的能量序列，以用于最终的频谱显示。
作为本发明一实施方式的进一步改进，所述步骤S1还包括：
预设M个存储空间，每个存储空间存储一组IQ信号序列，采用先进先出的方式依次存储获得M组IQ信号序列；
对M组IQ信号序列进行壁滤波处理后，组成M*N数据矩阵以用于依次在慢时间方向和快时间方向上做FFT变换。
作为本发明一实施方式的进一步改进，所述步骤S2具体包括：
P1、将所述M*N数据矩阵依次进行加窗，2D-FFT变换，并获取M*N数据矩阵中每个数据点的复数的模平方，以形成M*N的初始能量分布矩阵；
P2、对所述M*N的初始能量分布矩阵做转置处理，形成N*M的能量分布矩阵。
作为本发明一实施方式的进一步改进，所述步骤S3具体包括：
获取由V_Matrix(x_freq,y_freq)组成的N*M的速度分布矩阵；
当fsig(y_freq)<B_iq且fsig(y_freq)>-B_iq时，
V_Matrix(x_freq,y_freq)＝fd(x_freq)*c/(2*(fsig(y_freq)+f0))，
否则V_Matrix(x_freq,y_freq)＝-Prf*c/(2*f0)，
其中，V_Matrix(x_freq,y_freq)表示频移为fd(x_freq)＝x_freq/(N-1)*Prf、对应的频率分量为fsig(y_freq)＝y_freq/(M-1)*fs–fs/2处的速度大小，B_iq表示IQ信号序列进行壁滤波处理过程中壁滤波器的带宽，Prf为初始发射信号的重复扫查频率，f0为初始发射信号的中心频率，c为超声在组织内的传播速度，fs为信号采样频率。
作为本发明一实施方式的进一步改进，所述步骤S4具体包括：
M1、将所述用于显示的速度序列中每个数据点对应的初始能量值及初始记样点数均初始化为0；
M2、以用于显示的速度序列中的每个数据点依次查询速度分布矩阵，根据所述用于显示的速度序列中的每个数据点的速度大小，所述速度分布矩阵中每个矩阵速度点的大小，每个矩阵速度点的坐标值，按照预定规则获取累加系数，以及在所述能量分布矩阵中获取与各个累加系数对应的矩阵能量值；
M3、根据用于显示的速度序列中的每个数据点对应获得的所述矩阵能量值以及累加系数，获取所述用于显示的速度序列中的每个数据点对应的总能量值，以组成所述速度序列对应的能量序列，用于最终的频谱显示。
为了实现上述发明目的之一，本发明一实施方式提供一种脉冲多普勒超高谱分辨率成像处理***，
门采样模块，用于每根扫查线上均获取N个采样点分别对应的IQ信号，以组成每根扫查线对应的IQ信号序列；
壁滤波处理模块，用于在快时间方向上，分别对每根扫查线对应的IQ信号序列进行壁滤波处理，以组成每根扫查线对应的IQ信号壁滤波序列；
FFT变换模块，用于分别对IQ信号壁滤波序列中的每个采样点依次在慢时间方向和快时间方向上做FFT变换，以获得每个采样点在不同频移上的能量分布矩阵；
矩阵匹配模块，用于根据所述能量分布矩阵中各个采样点的频移坐标值、初始发射信号的重复扫查频率、初始发射信号的中心频率，获取每个采样点在不同频移上与所述能量分布矩阵匹配的速度分布矩阵；
处理输出模块，用于调取用于显示的速度序列，以所述速度序列查询所述速度分布矩阵和与所述速度分布矩阵匹配的能量分布矩阵，获取所述速度序列对应的能量序列，以用于最终的频谱显示。
作为本发明一实施方式的进一步改进，预设M个存储空间，每个存储空间存储一组IQ信号序列，所述门采样模块采用先进先出的方式依次存储获得M组IQ信号序列；
所述壁滤波处理模块对M组IQ信号序列进行壁滤波处理后，组成M*N数据矩阵以用于依次在慢时间方向和快时间方向上做FFT变换。
作为本发明一实施方式的进一步改进，所述FFT变换模块具体用于：
将所述M*N数据矩阵依次进行加窗，2D-FFT变换，并获取M*N数据矩阵中每个数据点的复数的模平方，以形成M*N的初始能量分布矩阵；
对所述M*N的初始能量分布矩阵做转置处理，形成N*M的能量分布矩阵。
作为本发明一实施方式的进一步改进，所述矩阵匹配模块具体用于：
获取由V_Matrix(x_freq,y_freq)组成的N*M的速度分布矩阵；
当fsig(y_freq)<B_iq且fsig(y_freq)>-B_iq时，
V_Matrix(x_freq,y_freq)＝fd(x_freq)*c/(2*(fsig(y_freq)+f0))，
否则V_Matrix(x_freq,y_freq)＝-Prf*c/(2*f0)，
其中，V_Matrix(x_freq,y_freq)表示频移为fd(x_freq)＝x_freq/(N-1)*Prf、对应的频率分量为fsig(y_freq)＝y_freq/(M-1)*fs–fs/2处的速度大小，B_iq表示IQ信号序列进行壁滤波处理过程中壁滤波器的带宽，Prf为初始发射信号的重复扫查频率，f0为初始发射信号的中心频率，c为超声在组织内的传播速度，fs为信号采样频率。
作为本发明一实施方式的进一步改进，所述处理输出模块具体用于：
M1、将所述用于显示的速度序列中每个数据点对应的初始能量值及初始记样点数均初始化为0；
M2、以用于显示的速度序列中的每个数据点依次查询速度分布矩阵，根据所述用于显示的速度序列中的每个数据点的速度大小，所述速度分布矩阵中每个矩阵速度点的大小，每个矩阵速度点的坐标值，按照预定规则获取累加系数，以及在所述能量分布矩阵中获取与各个累加系数对应的矩阵能量值。
M3、根据用于显示的速度序列中的每个数据点对应获得的所述矩阵能量值以及累加系数，获取所述用于显示的速度序列中的每个数据点对应的总能量值，以组成所述速度序列对应的能量序列，用于最终的频谱显示。
与现有技术相比，本发明的脉冲多普勒超高谱分辨率成像处理方法及处理***，对IQ信号序列壁滤波处理后做二维快速傅里叶变换，然后在IQ信号的有效带宽内，根据频移和信号的快时间频率分量，获取相应的能量分布矩阵及速度分布矩阵；最后在根据要显示的速度标尺，查询能量分布矩阵及速度分布矩阵，获取相应的能量序列用于显示成像；本发明在发射信号带宽较宽或者散射子速度较快时，速度分辨率都得到了很大保证，极大提高了速度谱分辨率，而且性噪比也得到了提升，提高了超声成像设备的方便性和使用效率，提升了超声图像的质量。
附图说明
图1是传统的成像***的整体模块示意图；
图2是传统的脉冲多普勒超高谱分辨率成像处理***的模块示意图；
图3是本发明一实施方式中脉冲多普勒超高谱分辨率成像处理方法的流程示意图；
图4是本发明一实施方式中脉冲多普勒超高谱分辨率成像处理***的模块示意图；
图5A、图5B是本发明一具体示例中分别采用传统的脉冲多普勒超高谱分辨率成像处理方法与本发明的脉冲多普勒超高谱分辨率成像处理方法获得的速度谱图的对比示意图；
图5C、图5D是本发明另一具体示例中分别采用传统的脉冲多普勒超高谱分辨率成像处理方法与本发明的脉冲多普勒超高谱分辨率成像处理方法获得的速度谱图的对比示意图。
具体实施方式
以下将结合附图所示的各实施方式对本发明进行详细描述。但这些实施方式并不限制本发明，本领域的普通技术人员根据这些实施方式所做出的结构、方法、或功能上的变换均包含在本发明的保护范围内。
需要说明的是，本发明主要应用于超声设备，相应的，所述待测物可为待测组织，在此不做详细赘述。
结合图1所示，多普勒成像***的模块示意图；脉冲多普勒超高谱分辨率成像过程中；
通过探头向组织中发射脉冲信号，所述脉冲信号经组织中反射形成超声信号经由探头换能器的不同基元转变为电模拟信号，通过前放模块放大，再由A/D数模转换模块转换为数字信号；各个不同基元的数字信号经过波束合成模块，合成为射频信号；射频信号经过固定频率的正交解调后，将正交解调结果I/Q信号送入相应的处理模块。
上述脉冲多普勒超高谱分辨率成像过程中；发射脉冲信号的扫查重复频率PRF、发射中心频率f0为***预设数值，对数字信号进行波束合成过程中，需要通过低通滤波器做正交解调，所述低通滤波器的解调带宽B_iq也为***预设阈值。
结合图2所示，传统的脉冲多普勒(PW)模式成像需经过如下过程：
在门门采样模块中，将在采样门里的信号采样点截出，再将截出的所有采样点求平均，送入后面的壁滤波模块，壁滤波一般为高通滤波器，主要用来在慢时间方向上虑除信号中低速运动的组织信号，传统方法有低阶FIR型滤波器、IIR型滤波器和多项式回归滤波器。在这些滤波器中，投影初始化IIR滤波器和多项式回归滤波器的性能具有阻带衰减大、过渡带窄和不损失数据等优点，其中投影初始化IIR滤波器的阻带截止频率更为灵活多变。壁滤波后的信号，就送入后面的FFT能量计算模块和双声道分离模块。
FFT能量计算模块，对壁滤波后的信号做快速傅里叶变换(FFT)，获得每个频移分量的能量大小，再通过动态范围压缩模块进行对数压缩，最后进行频谱显示。
双声道分离模块主要将信号的正负频谱进行左右声道分离，再经过速率转化为***声卡所需的采样率，最后进行音频播放。
结合图3所示，图3为本发明一实施方式中脉冲多普勒超高谱分辨率成像处理方法的流程图，所述方法包括：
S1、在每根扫查线上均获取N个采样点分别对应的IQ信号，以组成每根扫查线对应的IQ信号序列；
在快时间方向上，分别对每根扫查线对应的IQ信号序列进行壁滤波处理，以组成每根扫查线对应的IQ信号壁滤波序列。
本发明一具体示例中，为了方便描述，将每条回波信号上N个采样点分别对应的坐标点以(x,y)表示；
通常情况下，根据采样门大小在相应的回波信号上进行采样，所述采样门大小以SV表示。
相应的，对于第i根扫查线，i为正整数，其对应的IQ信号壁滤波序列可表示为：
{IQ(1,i),IQ(2,i),IQ(3,i),….IQ(N,i)}
其中，N＝2*SV/c*fs，c为超声在组织内的传播速度，fs为信号采样频率。
进一步的，本发明一优选实施方式中，所述步骤S1还包括：
预设M个存储空间，每个存储空间存储一组IQ信号序列，采用先进先出的方式依次存储获得M组IQ信号序列；对M组IQ信号序列进行壁滤波处理后，组成M*N数据矩阵以用于依次在慢时间方向和快时间方向上做FFT变换。
当M个存储空间全部存储数据后，下一个IQ信号序列需等待最前端存储空间中的IQ信号序列输出或清除后，再进行存储。
本发明具体实施方式中，对M组IQ信号序列进行壁滤波处理后，组成M*N数据矩阵为：
{IQ(1,i),IQ(2,i),IQ(3,i),…,IQ(N,i)}
{IQ(1,i+1),IQ(2,i+1),IQ(3,i+1),…,IQ(N,i+1)}
{IQ(1,i+2),IQ(2,i+2),IQ(3,i+2),…,IQ(N,i+2)}
……
{IQ(1,i+M-1),IQ(2,i+M-1),IQ(3,i+M-1),…,IQ(N,i+M-1)}
需要说明的是，在本发明的其他实施方式中，也可以将获取的壁滤波信号序列全部用于后续计算，在此不做详细赘述。
本发明一实施方式中，所述方法还包括：
S2、分别对IQ信号壁滤波序列中的每个采样点依次在慢时间方向和快时间方向上做FFT变换，以获得每个采样点在不同频移上的能量分布矩阵。
本实施方式中，分别对IQ信号壁滤波序列中的每个采样点依次在慢时间方向上做FFT变换，以分别获取IQ信号壁滤波序列中的每个采样点对应的不同频移上的能量分布；进一步的，在快时间方向上做FFT变换，可以进一步获得不同频移分量上的能量分布。
本发明具体实施方式中，所述步骤S2具体包括：
P1、将所述M*N数据矩阵依次进行加窗，2D-FFT变换，并获取M*N数据矩阵中每个数据点的复数的模平方，以形成M*N的初始能量分布矩阵；
相应的，形成的M*N的初始能量分布矩阵为：
{P_2DFFT(1,i),P_2DFFT(2,i),P_2DFFT(3,i),…,P_2DFFT(N,i)}
{P_2DFFT(1,i+1),P_2DFFT(2,i+1),P_2DFFT(3,i+1),…,P_2DFFT(N,i+1)}
{P_2DFFT(1,i+2),P_2DFFT(2,i+2),P_2DFFT(3,i+2),…,P_2DFFT(N,i+2)}
……
{P_2DFFT(1,i+M-1),P_2DFFT(2,i+M-1),P_2DFFT(3,i+M-1),…,P_2DFFT(N,i+M-1)}
进一步的，所述步骤S2还包括：
P2、对所述M*N的初始能量分布矩阵做转置处理，形成N*M的能量分布矩阵。
{P_2DFFT(1,i),P_2DFFT(1,i+1),P_2DFFT(1,i+2),…,P_2DFFT(1,i+M-1)}
{P_2DFFT(2,i),P_2DFFT(2,i+1),P_2DFFT(2,i+2),…,P_2DFFT(2,i+M-1)}
{P_2DFFT(3,i),P_2DFFT(3,i+1),P_2DFFT(3,i+2),…,P_2DFFT(3,i+M-1)}
……
{P_2DFFT(N,i),P_2DFFT(N,i+1),P_2DFFT(N,i+2),…,P_2DFFT(N,i+M-1)}
进一步的，本发明一实施方式中，所述方法还包括：
S3、根据所述能量分布矩阵中各个采样点的频移坐标值、初始发射信号的重复扫查频率、初始发射信号的中心频率，获取每个采样点在不同频移上与所述能量分布矩阵匹配的速度分布矩阵。
本发明一具体示例中，所述步骤S3具体包括：
获取由V_Matrix(x,y)组成的N*M的速度分布矩阵；与N*M的能量分布矩阵相匹配的速度分布矩阵为：
{V_MATRIX(1,i),V_MATRIX(1,i+1),V_MATRIX(1,i+2),…,V_MATRIX(1,i+M-1)}
{V_MATRIX(2,i),V_MATRIX(2,i+1),V_MATRIX(2,i+2),…,V_MATRIX(2,i+M-1)}
{V_MATRIX(3,i),V_MATRIX(3,i+1),V_MATRIX(3,i+2),…,V_MATRIX(3,i+M-1)}
……
{V_MATRIX(N,i),V_MATRIX(N,i+1),V_MATRIX(N,i+2),…,V_MATRIX(N,i+M-1)}
当fsig(y_freq)<B_iq且fsig(y_freq)>-B_iq时，
V_Matrix(x_freq,y_freq)＝fd(x_freq)*c/(2*(fsig(y_freq)+f0))，
否则V_Matrix(x_freq,y_freq)＝-Prf*c/(2*f0)，
其中，V_Matrix(x_freq,y_freq)表示频移为fd(x_freq)＝x_freq/(N-1)*Prf、对应的频率分量为fsig(y_freq)＝y_freq/(M-1)*fs–fs/2处的速度大小，B_iq表示IQ信号序列进行壁滤波处理过程中壁滤波器的带宽，Prf为初始发射信号的重复扫查频率，f0为初始发射信号的中心频率，c为超声在组织内的传播速度，fs为信号采样频率。
根据上述公式可知，所述速度分布矩阵中的各个数据点的大小与参数M、N、fs、Prf相关，其中M和fs在***中都是常数，故仅在采样门大小或者初始发射信号的重复扫查频率发生变化时，该速度分布矩阵中的各个数据点的大小才随之改变，在此不做详细赘述。
进一步的，本发明一实施方式中，所述方法还包括：
S4、调取用于显示的速度序列，以所述速度序列查询所述速度分布矩阵和与所述速度分布矩阵匹配的能量分布矩阵，获取所述速度序列对应的能量序列，以用于最终的频谱显示。
用于显示的速度序列可通过下述方式获得。在脉冲多普勒扫查时，其最大显示速度以V_max表示，所述用于显示的速度序列以V_Dis表示；
则：所述用于显示的速度序列中的每个数据
V_Dis(j)＝j*V_max/K，j＝(0,1,…K-1)
V_max＝Prf*c/(2*f0)
其中，K为***设置参数，其可以根据需要自行设定，在此不做详细赘述。
本发明一具体示例中，所述步骤S4具体包括：
M1、将所述用于显示的速度序列中每个数据点对应的初始能量值及初始记样点数均初始化为0。
M2、以用于显示的速度序列中的每个数据点依次查询速度分布矩阵，根据所述用于显示的速度序列中的每个数据点的速度大小，所述速度分布矩阵中每个矩阵速度点的大小，每个矩阵速度点的坐标值，按照预定规则获取累加系数，以及在所述能量分布矩阵中获取与各个累加系数对应的矩阵能量值；
M3、根据用于显示的速度序列中的每个数据点对应获得的所述矩阵能量值以及累加系数，获取所述用于显示的速度序列中的每个数据点对应的总能量值，以组成所述速度序列对应的能量序列，用于最终的频谱显示。
本发明一具体示例中，所述能量值以Power_V(j)表示，所述记样点数以V_Num(j)表示。
初始化后，每个数据点对应的初始能量值Power_V(0)＝0，初始记样点数V_Num(j)＝0。
以数据点V_Dis(j)依次查询速度分布矩阵的每一行，直至获得当前行中第一个大于或等于所述数据点V_Dis(j)的矩阵速度值V_Matrix(n_v,m_v)；
当然，若当前行未查询到大于或等于所述数据点V_Dis(j)的矩阵速度值V_Matrix(n_v,m_v)，则当前行对应V_Dis(j)的能量值为0，继续查询所述速度分布矩阵的下一行。
若所述数据点V_Dis(j)的值等于当前行中的矩阵速度值V_Matrix(n_v,m_v)，则直接将能量分布矩阵中、(n_v,m_v)的坐标值对应的矩阵能量值作为所述数据点V_Dis(j)在当前行的能量值；
若查询后，当前行中的速度值V_Matrix(n_v,m_v)大于所述数据点V_Dis(j)的值，则采用所述预定规则获取V_Dis(j)的累加系数；
所述预定规则，例如：采用插值法、绝对值法、三次方等方式。
为了便于理解，本发明仅以插值法为例做具体介绍。
该具体示例中，仅以查询所述速度分布矩阵的一行数据点为例做具体介绍。
对于所述速度分布矩阵的某一行：
c1＝Dis_b^2/(Dis_a^2+Dis_b^2)；
c2＝1-c1；
Dis_a＝V_Dis(i)-V_Matrix(n_v,m_v-1)；
Dis_b＝V_Matrix(n_v,m_v)-V_Dis(i)；
其中，c1、c2分别表示对应坐标点(n_v,m_v-1)、(n_v,m_v)的累加系数。
进一步的，根据每一行获得的坐标点所对应的能量分布矩阵中的矩阵能量值及其累加系数获得每个数据点在每一行对应的能量值，并将每一行获得的能量值进行累加，形成当前数据点对应的累加能量值。
∑Power_V(i)＝Power_V(0)+c1*P_2DFFT(n_v,m_v-1)+c2*P_2DFFT(n_v,m_v)+……
其中，Power_V(0)表示所述数据点的初始能量值，∑Power_V(i)表示所述数据点对应的累加能量值。
进一步的，每个数据点对应的总能量值，
Power_V(i)＝∑Power_V(i)/∑V_Num(i)，其中，Power_V(i)表示总能量值，∑V_Num(i)表示参与计算的矩阵能量值的数量。
当然，在本发明的其他实施方式中，还可以将与所述数据点V_Dis(j)对应的能量分布矩阵中的多个矩阵能量值平均值最小值、最大值、边值以及中值作为每个数据点V_Dis(j)对应的总能量，在此不做详细赘述。
进一步的，将每个数据点V_Dis(j)对应的总能量值组成能量序列，进行动态范围压缩，最后进行输出显示，在此不做详细赘述。
结合图5A-图5D所示，图5A为条件1下采用传统的脉冲多普勒超高谱分辨率成像处理方法获得的速度谱图；图5B为条件1下采用本发明的脉冲多普勒超高谱分辨率成像处理方法获得的速度谱图。图5C为条件2下采用传统的脉冲多普勒超高谱分辨率成像处理方法获得的速度谱图；图5D为条件2下采用本发明的脉冲多普勒超高谱分辨率成像处理方法获得的速度谱图。
相对于图5A、图5B，初始发射信号的中心频率为5MHz，发射周期数为20，初始发射信号的重复扫查频率为8KHz，三个散射子速度分别为0.1m/s,0.5m/s,0.8m/s，随机噪声水平为-20dB。
相对于图5C、图5D，初始发射信号的中心频率为5MHz，发射周期数为6，初始发射信号的重复扫查频率为8KHz，三个散射子速度分别为0.1m/s,0.5m/s,0.8m/s，随机噪声水平为-20dB。
分析图5A的速度谱图，可得出：散射子速度谱分辨随着速度的增大，逐渐变差，同时能量幅度也逐步下降。
分析图5A、5C的速度谱图，可得出：当发射信号带宽较宽时，传统处理方法的速度谱分辨率随之下降。
分析图5A、5B的速度谱图，可得出：本发明的速度谱图分辨率得到了极大提高，并且不受散射子速度影响。
分析图5C、5D的速度谱图，可得出：本发明的速度谱分辨率得到了极大提高的同时，不受发射信号带宽的影响；同时，频谱的能量不均一性也得到了解决，性噪比也得到了提升。
分析图5B、5D的速度谱图，可得出：本发明在发射带宽不一样的情况下，速度谱分辨率基本一致。
结合图4所示，本发明一实施方式中提供的脉冲多普勒超高谱分辨率成像处理***，所述***包括：门采样模块100、壁滤波处理模块200、FFT变换模块300、矩阵匹配模块400、处理输出模块500。
门采样模块100用于在每根扫查线上均获取N个采样点分别对应的IQ信号，以组成每根扫查线对应的IQ信号序列；
壁滤波处理模块200用于在快时间方向上，分别对每根扫查线对应的IQ信号序列进行壁滤波处理，以组成每根扫查线对应的IQ信号壁滤波序列。
本发明一具体示例中，为了方便描述，将每条回波信号上N个采样点分别对应的坐标点以(x,y)表示；
通常情况下，门采样模块100根据采样门大小在相应的回波信号上进行采样，所述采样门大小以SV表示。
相应的，对于第i根扫查线，i为正整数，其对应的IQ信号壁滤波序列可表示为：
{IQ(1,i),IQ(2,i),IQ(3,i),….IQ(N,i)}
其中，N＝2*SV/c*fs，c为超声在组织内的传播速度，fs为信号采样频率。
进一步的，本发明一优选实施方式中，所述步骤S1还包括：
预设M个存储空间，每个存储空间存储一组IQ信号序列，门采样模块100采用先进先出的方式依次存储获得M组IQ信号序列；壁滤波处理模块200对M组IQ信号序列进行壁滤波处理后，组成M*N数据矩阵以用于依次在慢时间方向和快时间方向上做FFT变换。
当M个存储空间全部存储数据后，下一个IQ信号序列需等待最前端存储空间中的IQ信号序列输出或清除后，再进行存储。
本发明具体实施方式中，壁滤波处理模块200对M组IQ信号序列进行壁滤波处理后，组成M*N数据矩阵为：
{IQ(1,i),IQ(2,i),IQ(3,i),…,IQ(N,i)}
{IQ(1,i+1),IQ(2,i+1),IQ(3,i+1),…,IQ(N,i+1)}
{IQ(1,i+2),IQ(2,i+2),IQ(3,i+2),…,IQ(N,i+2)}
……
{IQ(1,i+M-1),IQ(2,i+M-1),IQ(3,i+M-1),…,IQ(N,i+M-1)}
需要说明的是，在本发明的其他实施方式中，也可以将获取的壁滤波信号序列全部用于后续计算，在此不做详细赘述。
本发明一实施方式中，FFT变换模块300用于分别对IQ信号壁滤波序列中的每个采样点依次在慢时间方向和快时间方向上做FFT变换，以获得每个采样点在不同频移上的能量分布矩阵。
本实施方式中，FFT变换模块300用于分别对IQ信号壁滤波序列中的每个采样点依次在慢时间方向上做FFT变换，以分别获取IQ信号壁滤波序列中的每个采样点对应的不同频移上的能量分布；进一步的，FFT变换模块300用于在快时间方向上做FFT变换，可以进一步获得不同频移分量上的能量分布。
本发明具体实施方式中，FFT变换模块300具体用于将所述M*N数据矩阵依次进行加窗，2D-FFT变换，并获取M*N数据矩阵中每个数据点的复数的模平方，以形成M*N的初始能量分布矩阵；
相应的，形成的M*N的初始能量分布矩阵为：
{P_2DFFT(1,i),P_2DFFT(2,i),P_2DFFT(3,i),…,P_2DFFT(N,i)}
{P_2DFFT(1,i+1),P_2DFFT(2,i+1),P_2DFFT(3,i+1),…,P_2DFFT(N,i+1)}
{P_2DFFT(1,i+2),P_2DFFT(2,i+2),P_2DFFT(3,i+2),…,P_2DFFT(N,i+2)}
……
{P_2DFFT(1,i+M-1),P_2DFFT(2,i+M-1),P_2DFFT(3,i+M-1),…,P_2DFFT(N,i+M-1)}
进一步的，FFT变换模块300还用于：
对所述M*N的初始能量分布矩阵做转置处理，形成N*M的能量分布矩阵。
{P_2DFFT(1,i),P_2DFFT(1,i+1),P_2DFFT(1,i+2),…,P_2DFFT(1,i+M-1)}
{P_2DFFT(2,i),P_2DFFT(2,i+1),P_2DFFT(2,i+2),…,P_2DFFT(2,i+M-1)}
{P_2DFFT(3,i),P_2DFFT(3,i+1),P_2DFFT(3,i+2),…,P_2DFFT(3,i+M-1)}
……
{P_2DFFT(N,i),P_2DFFT(N,i+1),P_2DFFT(N,i+2),…,P_2DFFT(N,i+M-1)}
进一步的，本发明一实施方式中，矩阵匹配模块400用于：根据所述能量分布矩阵中各个采样点的频移坐标值、初始发射信号的重复扫查频率、初始发射信号的中心频率，获取每个采样点在不同频移上与所述能量分布矩阵匹配的速度分布矩阵。
本发明一具体示例中，矩阵匹配模块400具体用于：
获取由V_Matrix(x,y)组成的N*M的速度分布矩阵；与N*M的能量分布矩阵相匹配的速度分布矩阵为：
{V_MATRIX(1,i),V_MATRIX(1,i+1),V_MATRIX(1,i+2),…,V_MATRIX(1,i+M-1)}
{V_MATRIX(2,i),V_MATRIX(2,i+1),V_MATRIX(2,i+2),…,V_MATRIX(2,i+M-1)}
{V_MATRIX(3,i),V_MATRIX(3,i+1),V_MATRIX(3,i+2),…,V_MATRIX(3,i+M-1)}
……
{V_MATRIX(N,i),V_MATRIX(N,i+1),V_MATRIX(N,i+2),…,V_MATRIX(N,i+M-1)}
当fsig(y_freq)<B_iq且fsig(y_freq)>-B_iq时，
V_Matrix(x_freq,y_freq)＝fd(x_freq)*c/(2*(fsig(y_freq)+f0))，
否则V_Matrix(x_freq,y_freq)＝-Prf*c/(2*f0)，
V_Matrix(x_freq,y_freq)表示频移为fd(x_freq)＝x_freq/(N-1)*Prf、对应的频率分量为fsig(y_freq)＝y_freq/(M-1)*fs–fs/2处的速度大小，B_iq表示IQ信号序列进行壁滤波处理过程中壁滤波器的带宽，Prf为初始发射信号的重复扫查频率，f0为初始发射信号的中心频率，c为超声在组织内的传播速度，fs为信号采样频率。
根据上述公式可知，所述速度分布矩阵中的各个数据点的大小与参数M、N、fs、Prf相关，其中M和fs在***中都是常数，故仅在采样门大小或者初始发射信号的重复扫查频率发生变化时，该速度分布矩阵中的各个数据点的大小才随之改变，在此不做详细赘述。
进一步的，本发明一实施方式中，处理输出模块500用于：调取用于显示的速度序列，以所述速度序列查询所述速度分布矩阵和与所述速度分布矩阵匹配的能量分布矩阵，获取所述速度序列对应的能量序列，以用于最终的频谱显示。
用于显示的速度序列可通过下述方式获得。在脉冲多普勒扫查时，其最大显示速度以V_max表示，所述用于显示的速度序列以V_Dis表示；
则：所述用于显示的速度序列中的每个数据
V_Dis(j)＝j*V_max/K，j＝(0,1,…K-1)V_max＝Prf*c/(2*f0)
其中，K为***设置参数，其可以根据需要自行设定，在此不做详细赘述。
本发明一具体示例中，处理输出模块500具体用于：
将所述用于显示的速度序列中每个数据点对应的初始能量值及初始记样点数均初始化为0。
以用于显示的速度序列中的每个数据点依次查询速度分布矩阵，根据所述用于显示的速度序列中的每个数据点的速度大小，所述速度分布矩阵中每个矩阵速度点的大小，每个矩阵速度点的坐标值，按照预定规则获取累加系数，以及在所述能量分布矩阵中获取与各个累加系数对应的矩阵能量值；
根据用于显示的速度序列中的每个数据点对应获得的所述矩阵能量值以及累加系数，获取所述用于显示的速度序列中的每个数据点对应的总能量值，以组成所述速度序列对应的能量序列，用于最终的频谱显示。
本发明一具体示例中，所述能量值以Power_V(j)表示，所述记样点数以V_Num(j)表示。
初始化后，每个数据点对应的初始能量值Power_V(0)＝0，初始记样点数V_Num(j)＝0。
处理输出模块500以数据点V_Dis(j)依次查询速度分布矩阵的每一行，直至获得当前行中第一个大于或等于所述数据点V_Dis(j)的矩阵速度值V_Matrix(n_v,m_v)；
当然，若当前行未查询到大于或等于所述数据点V_Dis(j)的矩阵速度值V_Matrix(n_v,m_v)，则当前行对应V_Dis(j)的能量值为0，继续查询所述速度分布矩阵的下一行。
若所述数据点V_Dis(j)的值等于当前行中的矩阵速度值V_Matrix(n_v,m_v)，则直接将能量分布矩阵中、(n_v,m_v)的坐标值对应的矩阵能量值作为所述数据点V_Dis(j)在当前行的能量值；
若查询后，当前行中的速度值V_Matrix(n_v,m_v)大于所述数据点V_Dis(j)的值，则采用所述预定规则获取V_Dis(j)的累加系数；
所述预定规则，例如：采用插值法、绝对值法、三次方等方式。
为了便于理解，本发明仅以插值法为例做具体介绍。
该具体示例中，仅以处理输出模块500查询所述速度分布矩阵的一行数据点为例做具体介绍。
对于所述速度分布矩阵的某一行：
c1＝Dis_b^2/(Dis_a^2+Dis_b^2)；
c2＝1-c1；
Dis_a＝V_Dis(i)-V_Matrix(n_v,m_v-1)；
Dis_b＝V_Matrix(n_v,m_v)-V_Dis(i)；
其中，c1、c2分别表示对应坐标点(n_v,m_v-1)、(n_v,m_v)的累加系数。
进一步的，处理输出模块500还用于：根据每一行获得的坐标点所对应的能量分布矩阵中的矩阵能量值及其累加系数获得每个数据点在每一行对应的能量值，并将每一行获得的能量值进行累加，形成当前数据点对应的累加能量值。
∑Power_V(i)＝Power_V(0)+c1*P_2DFFT(n_v,m_v-1)+c2*P_2DFFT(n_v,m_v)+……
其中，Power_V(0)表示所述数据点的初始能量值，∑Power_V(i)表示所述数据点对应的累加能量值。
进一步的，每个数据点对应的总能量值，
Power_V(i)＝∑Power_V(i)/∑V_Num(i)，其中，Power_V(i)表示总能量值，∑V_Num(i)表示参与计算的矩阵能量值的数量。
当然，在本发明的其他实施方式中，还可以将与所述数据点V_Dis(j)对应的能量分布矩阵中的多个矩阵能量值平均值最小值、最大值、边值以及中值作为每个数据点V_Dis(j)对应的总能量，在此不做详细赘述。
进一步的，处理输出模块500将每个数据点V_Dis(j)对应的总能量值组成能量序列，进行动态范围压缩，最后进行输出显示，在此不做详细赘述。
综上所述，本发明的脉冲多普勒超高谱分辨率成像处理方法及处理***，本发明的脉冲多普勒超高谱分辨率成像处理方法及处理***，对IQ信号序列壁滤波处理后做二维快速傅里叶变换，然后在IQ信号的有效带宽内，根据频移和信号的快时间频率分量，获取相应的能量分布矩阵及速度分布矩阵；最后在根据要显示的速度标尺，查询能量分布矩阵及速度分布矩阵，获取相应的能量序列用于显示成像；本发明在发射信号带宽较宽或者散射子速度较快时，速度分辨率都得到了很大保证，极大提高了速度谱分辨率，而且性噪比也得到了提升，提高了超声成像设备的方便性和使用效率，提升了超声图像的质量。
为了描述的方便，描述以上装置时以功能分为各种模块分别描述。当然，在实施本申请时可以把各模块的功能在同一个或多个软件和/或硬件中实现。
通过以上的实施方式的描述可知，本领域的技术人员可以清楚地了解到本申请可借助软件加必需的通用硬件平台的方式来实现。基于这样的理解，本申请的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该软件产品可以保存在保存介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，信息推送服务器，或者网络设备等)执行本申请各个实施方式或者实施方式的某些部分所述的方法。
以上所描述的装置实施方式仅仅是示意性的，其中所述作为分离部件说明的模块可以是或者也可以不是物理上分开的，作为模块显示的部件可以是或者也可以不是物理模块，即可以位于一个地方，或者也可以分布到多个网络模块上。可以根据实际的需要选择其中的部分或者全部模块来实现本实施方式方案的目的。本领域普通技术人员在不付出创造性劳动的情况下，即可以理解并实施。
本申请可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构等等。也可以在分布式计算环境中实践本申请，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括保存设备在内的本地和远程计算机保存介质中。
应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施方式中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。
上文所列出的一系列的详细说明仅仅是针对本发明的可行性实施方式的具体说明，它们并非用以限制本发明的保护范围，凡未脱离本发明技艺精神所作的等效实施方式或变更均应包含在本发明的保护范围之内。

Claims (10)

1.一种脉冲多普勒超高谱分辨率成像处理方法，其特征在于，所述方法包括以下步骤：
S1、在每根扫查线上均获取N个采样点分别对应的IQ信号，以组成每根扫查线对应的IQ信号序列；
在快时间方向上，分别对每根扫查线对应的IQ信号序列进行壁滤波处理，以组成每根扫查线对应的IQ信号壁滤波序列；
S2、分别对IQ信号壁滤波序列中的每个采样点依次在慢时间方向和快时间方向上做FFT变换，以获得每个采样点在不同频移上的能量分布矩阵；
S3、根据所述能量分布矩阵中各个采样点的频移坐标值、初始发射信号的重复扫查频率、初始发射信号的中心频率，获取每个采样点在不同频移上与所述能量分布矩阵匹配的速度分布矩阵；
S4、调取用于显示的速度序列，以所述速度序列查询所述速度分布矩阵和与所述速度分布矩阵匹配的能量分布矩阵，获取所述速度序列对应的能量序列，以用于最终的频谱显示。
2.根据权利要求1所述的脉冲多普勒超高谱分辨率成像处理方法，其特征在于，所述步骤S1还包括：
预设M个存储空间，每个存储空间存储一组IQ信号序列，采用先进先出的方式依次存储获得M组IQ信号序列；
对M组IQ信号序列进行壁滤波处理后，组成M*N数据矩阵以用于依次在慢时间方向和快时间方向上做FFT变换。
3.根据权利要求2所述的脉冲多普勒超高谱分辨率成像处理方法，其特征在于，所述步骤S2具体包括：
P1、将所述M*N数据矩阵依次进行加窗，2D-FFT变换，并获取M*N数据矩阵中每个数据点的复数的模平方，以形成M*N的初始能量分布矩阵；
P2、对所述M*N的初始能量分布矩阵做转置处理，形成N*M的能量分布矩阵。
4.根据权利要求3所述的脉冲多普勒超高谱分辨率成像处理方法，其特征在于，所述步骤S3具体包括：
获取由V_Matrix(x_freq,y_freq)组成的N*M的速度分布矩阵；
当fsig(y_freq)<B_iq且fsig(y_freq)>-B_iq时，
V_Matrix(x_freq,y_freq)＝fd(x_freq)*c/(2*(fsig(y_freq)+f0))，
否则V_Matrix(x_freq,y_freq)＝-Prf*c/(2*f0)，
其中，V_Matrix(x_freq,y_freq)表示频移为fd(x_freq)＝x_freq/(N-1)*Prf、对应的频率分量为fsig(y_freq)＝y_freq/(M-1)*fs–fs/2处的速度大小，B_iq表示IQ信号序列进行壁滤波处理过程中壁滤波器的带宽，Prf为初始发射信号的重复扫查频率，f0为初始发射信号的中心频率，c为超声在组织内的传播速度，fs为信号采样频率。
5.根据权利要求1所述的脉冲多普勒超高谱分辨率成像处理方法，其特征在于，所述步骤S4具体包括：
M1、将所述用于显示的速度序列中每个数据点对应的初始能量值及初始记样点数均初始化为0；
M2、以用于显示的速度序列中的每个数据点依次查询速度分布矩阵，根据所述用于显示的速度序列中的每个数据点的速度大小，所述速度分布矩阵中每个矩阵速度点的大小，每个矩阵速度点的坐标值，按照预定规则获取累加系数，以及在所述能量分布矩阵中获取与各个累加系数对应的矩阵能量值；
M3、根据用于显示的速度序列中的每个数据点对应获得的所述矩阵能量值以及累加系数，获取所述用于显示的速度序列中的每个数据点对应的总能量值，以组成所述速度序列对应的能量序列，用于最终的频谱显示。
6.一种脉冲多普勒超高谱分辨率成像处理***，其特征在于，所述***包括：
门采样模块，用于每根扫查线上均获取N个采样点分别对应的IQ信号，以组成每根扫查线对应的IQ信号序列；
壁滤波处理模块，用于在快时间方向上，分别对每根扫查线对应的IQ信号序列进行壁滤波处理，以组成每根扫查线对应的IQ信号壁滤波序列；
FFT变换模块，用于分别对IQ信号壁滤波序列中的每个采样点依次在慢时间方向和快时间方向上做FFT变换，以获得每个采样点在不同频移上的能量分布矩阵；
矩阵匹配模块，用于根据所述能量分布矩阵中各个采样点的频移坐标值、初始发射信号的重复扫查频率、初始发射信号的中心频率，获取每个采样点在不同频移上与所述能量分布矩阵匹配的速度分布矩阵；
处理输出模块，用于调取用于显示的速度序列，以所述速度序列查询所述速度分布矩阵和与所述速度分布矩阵匹配的能量分布矩阵，获取所述速度序列对应的能量序列，以用于最终的频谱显示。
7.根据权利要求6所述的脉冲多普勒超高谱分辨率成像处理***，其特征在于，
预设M个存储空间，每个存储空间存储一组IQ信号序列，所述门采样模块采用先进先出的方式依次存储获得M组IQ信号序列；
所述壁滤波处理模块对M组IQ信号序列进行壁滤波处理后，组成M*N数据矩阵以用于依次在慢时间方向和快时间方向上做FFT变换。
8.根据权利要求7所述的脉冲多普勒超高谱分辨率成像处理***，其特征在于，所述FFT变换模块具体用于：
将所述M*N数据矩阵依次进行加窗，2D-FFT变换，并获取M*N数据矩阵中每个数据点的复数的模平方，以形成M*N的初始能量分布矩阵；
对所述M*N的初始能量分布矩阵做转置处理，形成N*M的能量分布矩阵。
9.根据权利要求8所述的脉冲多普勒超高谱分辨率成像处理***，其特征在于，所述矩阵匹配模块具体用于：
获取由V_Matrix(x_freq,y_freq)组成的N*M的速度分布矩阵；
当fsig(y_freq)<B_iq且fsig(y_freq)>-B_iq时，
V_Matrix(x_freq,y_freq)＝fd(x_freq)*c/(2*(fsig(y_freq)+f0))，
否则V_Matrix(x_freq,y_freq)＝-Prf*c/(2*f0)，
其中，V_Matrix(x_freq,y_freq)表示频移为fd(x_freq)＝x_freq/(N-1)*Prf、对应的频率分量为fsig(y_freq)＝y_freq/(M-1)*fs–fs/2处的速度大小，B_iq表示IQ信号序列进行壁滤波处理过程中壁滤波器的带宽，Prf为初始发射信号的重复扫查频率，f0为初始发射信号的中心频率，c为超声在组织内的传播速度，fs为信号采样频率。
10.根据权利要求6所述的脉冲多普勒超高谱分辨率成像处理***，其特征在于，所述处理输出模块具体用于：
M1、将所述用于显示的速度序列中每个数据点对应的初始能量值及初始记样点数均初始化为0；
M2、以用于显示的速度序列中的每个数据点依次查询速度分布矩阵，根据所述用于显示的速度序列中的每个数据点的速度大小，所述速度分布矩阵中每个矩阵速度点的大小，每个矩阵速度点的坐标值，按照预定规则获取累加系数，以及在所述能量分布矩阵中获取与各个累加系数对应的矩阵能量值；
M3、根据用于显示的速度序列中的每个数据点对应获得的所述矩阵能量值以及累加系数，获取所述用于显示的速度序列中的每个数据点对应的总能量值，以组成所述速度序列对应的能量序列，用于最终的频谱显示。
CN201610356317.8A 2016-05-26 2016-05-26 脉冲多普勒超高谱分辨率成像处理方法及处理*** Active CN105997148B (zh)

Priority Applications (4)

Application Number Priority Date Filing Date Title
CN201610356317.8A CN105997148B (zh) 2016-05-26 2016-05-26 脉冲多普勒超高谱分辨率成像处理方法及处理***
PCT/CN2017/073774 WO2017202068A1 (zh) 2016-05-26 2017-02-16 脉冲多普勒超高谱分辨率成像处理方法及处理***
EP17801924.6A EP3466343B1 (en) 2016-05-26 2017-02-16 Pulse doppler ultrahigh spectrum resolution imaging processing method and processing system
US16/304,268 US11197656B2 (en) 2016-05-26 2017-02-16 Pulsed doppler ultra-high spectral resolution imaging processing method and processing system

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610356317.8A CN105997148B (zh) 2016-05-26 2016-05-26 脉冲多普勒超高谱分辨率成像处理方法及处理***

Publications (2)

Publication Number Publication Date
CN105997148A CN105997148A (zh) 2016-10-12
CN105997148B true CN105997148B (zh) 2019-01-29

Family

ID=57094165

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610356317.8A Active CN105997148B (zh) 2016-05-26 2016-05-26 脉冲多普勒超高谱分辨率成像处理方法及处理***

Country Status (4)

Country Link
US (1) US11197656B2 (zh)
EP (1) EP3466343B1 (zh)
CN (1) CN105997148B (zh)
WO (1) WO2017202068A1 (zh)

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105997147B (zh) * 2016-07-28 2019-04-16 飞依诺科技(苏州)有限公司 一种超声波脉冲多普勒成像方法及装置
CN106991708B (zh) * 2017-04-27 2020-04-14 飞依诺科技(苏州)有限公司 超声多普勒血流成像的处理方法及处理***
CN106955098B (zh) * 2017-05-05 2020-10-16 飞依诺科技(苏州)有限公司 一种血管流速计算方法及装置
CN110349256B (zh) * 2019-07-16 2023-05-23 深圳大学 血管重建方法、装置及计算机终端
CN110750347A (zh) * 2019-10-22 2020-02-04 上海创远仪器技术股份有限公司 长曝光对比频谱数据处理***及其方法
CN111898476B (zh) * 2020-07-12 2022-04-26 西北工业大学 一种耦合随机共振的自适应线谱增强方法
CN112087272B (zh) * 2020-08-04 2022-07-19 中电科思仪科技股份有限公司 一种电磁频谱监测接收机信号自动检测方法
WO2022051992A1 (zh) * 2020-09-10 2022-03-17 华为技术有限公司 基于回波信号的速度探测方法和装置
CN113822329B (zh) * 2021-08-10 2023-11-03 国网新源控股有限公司 一种水电机组主轴摆度信号处理方法及装置
CN114376606B (zh) * 2022-01-18 2023-05-09 武汉联影医疗科技有限公司 一种超声成像的滤波方法和***
CN114578093B (zh) * 2022-03-10 2023-08-18 中国计量科学研究院 一种基于混合基fft的激光多普勒测速仪测速方法
CN116982942B (zh) * 2023-08-07 2024-04-02 武汉迅检科技有限公司 一种基于红外热成像的口腔测温方法及***

Family Cites Families (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4930513A (en) * 1988-07-26 1990-06-05 U.S. Philips Corporation Two dimensional processing of pulsed Doppler signals
US5122803A (en) * 1991-11-06 1992-06-16 The United States Of America As Represented By The Secretary Of The Army Moving target imaging synthetic aperture radar
US6390984B1 (en) * 2000-09-14 2002-05-21 Ge Medical Systems Global Technology Company, Llc Method and apparatus for locking sample volume onto moving vessel in pulsed doppler ultrasound imaging
FR2833445B1 (fr) * 2001-12-12 2004-02-27 Pascal Hannequin Procede et dispositif de traitement d'images a bas niveau
CN100367915C (zh) * 2004-03-15 2008-02-13 深圳迈瑞生物医疗电子股份有限公司 基于反向初始化iir的彩色血流成像壁滤波方法和装置
RU2536418C2 (ru) * 2009-05-13 2014-12-20 Конинклейке Филипс Электроникс Н.В. Ультразвуковое допплеровское аудиоустройство контроля кровотока со смещением основного тона
WO2012131340A2 (en) * 2011-03-25 2012-10-04 Norwegian University Of Science And Technology (Ntnu) Methods and apparatus for multibeam doppler ultrasound display
US20140088429A1 (en) * 2011-05-25 2014-03-27 Orcasonix Ltd. Ultrasound imaging system and method
CN103142253B (zh) * 2013-03-19 2014-11-12 飞依诺科技(苏州)有限公司 超声成像***及该***的波束叠加方法
CN103399301B (zh) * 2013-07-01 2015-08-05 北京航空航天大学 一种宽带sar信号的接收装置及接收方法
CN103675759B (zh) * 2013-11-27 2016-03-09 杭州电子科技大学 一种改进的分数阶傅里叶变换机动弱目标检测方法
JP6342212B2 (ja) * 2014-05-12 2018-06-13 キヤノンメディカルシステムズ株式会社 超音波診断装置
CN105476665B (zh) * 2016-01-27 2019-01-25 成都思多科医疗科技有限公司 一种基于超声的血流速度测量及血流流量测量方法

Also Published As

Publication number Publication date
US11197656B2 (en) 2021-12-14
WO2017202068A1 (zh) 2017-11-30
US20190083066A1 (en) 2019-03-21
EP3466343A4 (en) 2019-06-19
EP3466343A1 (en) 2019-04-10
EP3466343B1 (en) 2021-03-03
CN105997148A (zh) 2016-10-12

Similar Documents

Publication Publication Date Title
CN105997148B (zh) 脉冲多普勒超高谱分辨率成像处理方法及处理***
EP3132281B1 (en) Ultrasonic imaging compression methods and apparatus
US8684934B2 (en) Adaptively performing clutter filtering in an ultrasound system
US8858446B2 (en) Color doppler ultrasonic diagnosis apparatus which can calculate blood flow component information
US11346929B2 (en) Systems and methods for ultrafast ultrasound imaging
US11432793B2 (en) High resolution compound ultrasound flow imaging
CN107303186B (zh) 弹性成像中的频率复合
EP2610641B1 (en) Ultrasound and system for forming a Doppler ultrasound image
Ramalli et al. A real-time chirp-coded imaging system with tissue attenuation compensation
CN105919625B (zh) 脉冲多普勒壁滤波处理方法及处理***
JP3093823B2 (ja) 超音波ドプラ診断装置
CN105212964B (zh) 基于rf数据超声成像处理方法及***
JP6305470B2 (ja) 改善されたドップライメージング
JP5247322B2 (ja) 超音波撮像装置
CN108652666B (zh) 一种多普勒血流成像的生成方法和装置
JP2007507271A (ja) 超音波イメージングにおけるスモールアンサンブル長によるクラッタフィルタリング
CN105919624B (zh) 一种高脉冲重复扫查频率信号的编码和装置
EP2341364B1 (en) Adaptive clutter filtering method and ultrasound system for the same
Di Ianni et al. Real-time implementation of synthetic aperture vector flow imaging on a consumer-level tablet
KR101809358B1 (ko) 새로운 평면파 합성을 이용한 초음파 도플러 영상 장치 및 그 제어 방법
WO2018172882A1 (en) Doppler ultrasound using sparse emission pattern
Schiffner et al. A low-rate parallel Fourier domain beamforming method for ultrafast pulse-echo imaging
US20240027406A1 (en) Method and system for imaging a target from coherent waves
Bar-Zion et al. Towards sub-Nyquist tissue Doppler imaging using non-uniformly spaced stream of pulses
Bayat et al. Pixel-based ultrasound image reconstruction: Impact of grid size on signal frequency content

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CP01 Change in the name or title of a patent holder

Address after: 215123 5th floor, building a, 4th floor, building C, No. 27, Xinfa Road, Suzhou Industrial Park, Jiangsu Province

Patentee after: Feiyinuo Technology Co.,Ltd.

Address before: 215123 5th floor, building a, 4th floor, building C, No. 27, Xinfa Road, Suzhou Industrial Park, Jiangsu Province

Patentee before: Feiyinuo Technology (Suzhou) Co.,Ltd.

CP01 Change in the name or title of a patent holder
CP03 Change of name, title or address

Address after: 215123 5th floor, building a, 4th floor, building C, No. 27, Xinfa Road, Suzhou Industrial Park, Jiangsu Province

Patentee after: Feiyinuo Technology (Suzhou) Co.,Ltd.

Address before: 215123 5F, building a, No. 27, Xinfa Road, Suzhou Industrial Park, Jiangsu Province

Patentee before: VINNO TECHNOLOGY (SUZHOU) Co.,Ltd.

CP03 Change of name, title or address