CN110559014A - 基于探头补偿的分数阶傅里叶变换回波成像方法与系统 - Google Patents
基于探头补偿的分数阶傅里叶变换回波成像方法与系统 Download PDFInfo
- Publication number
- CN110559014A CN110559014A CN201910803651.7A CN201910803651A CN110559014A CN 110559014 A CN110559014 A CN 110559014A CN 201910803651 A CN201910803651 A CN 201910803651A CN 110559014 A CN110559014 A CN 110559014A
- Authority
- CN
- China
- Prior art keywords
- signal
- probe
- fourier transform
- echo
- imaging
- 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.)
- Granted
Links
- 239000000523 sample Substances 0.000 title claims abstract description 147
- 238000003384 imaging method Methods 0.000 title claims abstract description 81
- 230000004044 response Effects 0.000 claims abstract description 83
- 238000000034 method Methods 0.000 claims abstract description 43
- 238000007906 compression Methods 0.000 claims abstract description 39
- 230000006835 compression Effects 0.000 claims abstract description 39
- 238000001514 detection method Methods 0.000 claims abstract description 37
- 238000005259 measurement Methods 0.000 claims abstract description 6
- 239000011159 matrix material Substances 0.000 claims description 24
- 238000006243 chemical reaction Methods 0.000 claims description 20
- 230000014509 gene expression Effects 0.000 claims description 17
- 238000001228 spectrum Methods 0.000 claims description 12
- 238000005070 sampling Methods 0.000 claims description 8
- 238000002604 ultrasonography Methods 0.000 claims description 8
- 230000002776 aggregation Effects 0.000 claims description 6
- 238000004220 aggregation Methods 0.000 claims description 6
- 239000007822 coupling agent Substances 0.000 claims description 6
- 239000000126 substance Substances 0.000 claims description 5
- 230000008878 coupling Effects 0.000 claims description 3
- 238000010168 coupling process Methods 0.000 claims description 3
- 238000005859 coupling reaction Methods 0.000 claims description 3
- 238000010276 construction Methods 0.000 claims description 2
- 238000005457 optimization Methods 0.000 claims description 2
- 238000005516 engineering process Methods 0.000 description 9
- 238000010606 normalization Methods 0.000 description 7
- 210000001519 tissue Anatomy 0.000 description 4
- 238000012285 ultrasound imaging Methods 0.000 description 4
- 238000009825 accumulation Methods 0.000 description 3
- 238000001914 filtration Methods 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000002059 diagnostic imaging Methods 0.000 description 1
- 230000002526 effect on cardiovascular system Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000000004 hemodynamic effect Effects 0.000 description 1
- 230000010365 information processing Effects 0.000 description 1
- 210000000056 organ Anatomy 0.000 description 1
- 210000004872 soft tissue Anatomy 0.000 description 1
- 238000004148 unit process Methods 0.000 description 1
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/44—Constructional features of the ultrasonic, sonic or infrasonic diagnostic device
- A61B8/4444—Constructional features of the ultrasonic, sonic or infrasonic diagnostic device related to the probe
-
- 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
-
- 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
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Heart & Thoracic Surgery (AREA)
- Molecular Biology (AREA)
- Biophysics (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Pathology (AREA)
- Radiology & Medical Imaging (AREA)
- Biomedical Technology (AREA)
- Veterinary Medicine (AREA)
- Medical Informatics (AREA)
- Physics & Mathematics (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Ultra Sonic Daignosis Equipment (AREA)
- Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)
Abstract
本发明提供一种基于探头补偿的分数阶傅里叶变换回波成像方法与系统，其特征在于先进行超声探头的频响补偿，再利用LFM信号在分数阶傅里叶域中的能量聚集特性进行脉冲压缩并成像。所述超声探头的频响补偿包括超声探头频响测定与回波信号探头频响补偿两个部分。超声探头频响测定可以采用基于时域、频域或时频域的方法来求解超声探头的冲激响应；对完成补偿后的回波信号做最佳阶数的分数阶傅里叶变换变换至能量聚集域，实现脉冲压缩，对脉冲压缩后的结果进行包络检测，再将分数阶傅里叶变换域横坐标u转换为时延τ，并将幅值转换为灰度值，最后根据时延和灰度进行成像。本发明具有回波检测精度高、成像质量高、抗噪性能好等优点。
Description
技术领域
本发明涉及超声成像技术领域，更具体地，涉及一种基于探头补偿的分数阶傅里叶变换回波成像方法与系统。
背景技术
在医学超声诊断系统中，超声成像技术是最重要的技术之一。医学超声成像因其具有强度低、对人体无损伤、成本低和易于使用的特点而被广泛使用。尤其是在人体软组织的探测和心血管脏器的血流动力学观察方面具有不可替代的优势。随着电子学、生物工程学、图像与信息处理技术和计算机技术等学科的迅猛发展，超声成像已经发展成为一门成熟的医学影像学科，具有非常广阔的发展前景。
现今的超声成像技术主要利用了信号传播的时延特性：信号由发射探头发出，由不同位置的反射体反射产生回波。接收探头接收到回波后，估计时延并利用时延信息来表征反射体的实际位置进行成像。现有的常用的探测技术是基于短脉冲信号的探测，信号持续时间短、占空比低，此种方法会导致以下不足：
(1)空间分辨率低。当前的超声成像系统采用短脉冲基波成像原理，空间分辨率为毫米量级，难以区分距离更近的成像目标；
(2)探测深度短。脉冲持续时间越短，所携带的能量就越低，此种方法使用的短脉冲信号将导致较短的探测深度；
(3)图像质量差。超声成像基本原理是从大量背景噪声中提取有用的信息。此种方法受噪声干扰严重，导致图像清晰度、对比度差。
根据香农信息论和参数估计理论，信号的时宽-带宽积越大，所携带的信息量越多，回波检测的精度以及成像的精度就越高。在医学超声探测成像技术中的脉冲编码技术就是为了拓展信号的时宽-带宽积而提出的。线性调频(LFM)信号作为一种时宽-带宽积大的信号，应用于超声成像系统中被证明可以提高成像精度，然而现有的基于LFM信号的探测在接收端进行回波检测时依然存在下述问题，使得LFM信号在医学探测中尚未广泛应用：
(1)探头幅频响应对探测信号的影响。探头自身的脉冲响应相当于一个通带较小的带通滤波器，可近似为高斯包络的余弦波，会对发射信号的低频和高频成分有比较大地衰减，使得回波信号的带宽受限于探头的带宽。假设发射超声探头与接收超声探头的的冲激响应均为h1(t)，若系统发送的宽频模拟信号为x(t)，探测时，若不对探头进行补偿，接收到的回波信号y(t)可表示为：
y(t)＝x(t-τ)*h1(t)*h1(t)＝x(t-τ)*h(t)；
即相当于发送信号x(t)经过了一个冲激响应为h(t)的信道得到回波信号y(t)，若此时直接对接收信号进行匹配回波检测，将影响最后的成像效果。
(2)回波检测技术精度不高。超声系统成像时往往需要对回波信号先进行脉冲压缩。传统的脉冲压缩方法是匹配滤波，匹配滤波器的单位冲激响应是发送信号的时反序列，因此匹配滤波过程本质是信号间的互相关，没有利用到信号变换域的特点，在信噪比高的时候性能尚可，然而医学超声经过多层人体组织的反射和传播，探测目标信号微弱，导致成像精度不够，与传统脉冲法相比成像质量并无显著提升，这也是目前LFM信号尚未在医学超声探测领域得到大量使用的主要原因。
LFM信号在分数阶傅里叶域中的特殊性质十分利于回波信号的时延估计与脉冲压缩。分数阶傅里叶变换可以理解是傅里叶变换的一种推广形式，其中一点不同是分数阶傅里叶变换中具有阶数p。分数阶傅里叶域是通过分数阶傅里叶变换操作对时频域同时进行一定角度的旋转，其中旋转的角度与阶数满足α＝pπ/2，旋转后的域仍然具有正交性。LFM信号在分数阶傅里叶变换中可旋转一个特殊的角度至能量聚集域中，在该域里LFM信号具有良好的能量聚集特性，即该信号可在分数阶傅里叶域里形成一个冲激的形式，使其能很好地与噪声进行分离。经过时延后的LFM信号在分数阶傅里叶域里形成冲激的位置坐标与时延和信号的参数有关，在发送信号选定的情况下，冲激的位置坐标就只与时延有关。
为将大时宽-带宽积的LFM信号应用于超声成像系统中，提高成像精度，本发明提出的思路是，首先对探头的频响特性进行均衡，然后利用分数阶傅里叶变换实现脉冲压缩，再进行成像。
发明内容
针对目前超声探测系统对LFM信号的限制问题，本发明提供一种基于探头补偿的分数阶傅里叶变换回波成像方法与系统，以达到更高的成像精度。
本发明的目的至少通过如下技术方案之一实现。
基于探头补偿的分数阶傅里叶变换回波成像方法，包括以下步骤：
S1、先进行超声探头的频响补偿，包括超声探头频响测定与回波信号探头频响补偿两个部分；
S2、再进行基于分数阶傅里叶变换的脉冲压缩与成像。
进一步地，所述超声探头频响测定的步骤如下：
S1.1.1、将发射探头与接收探头通过耦合剂耦合，两个超声探头总的冲激响应用h(n)表示；
S1.1.2、利用数据处理单元中的信号发送单元发送一个已知宽带线性调频信号即发送信号x(n)，由D/A转换器完成数模转换，再通过发射探头转换为声信号，该声信号经过耦合剂后由接收探头接收并转换为电信号，再由A/D转换器完成模数转换得到接收信号y(n)，将其发送给数据处理单元中的信号接收单元；
S1.1.3、根据接收信号y(n)和发送信号x(n)求解超声探头的冲激响应h(n)。
进一步地，求取探头冲激响应的方法可采用频域最优逼近法、时域最优逼近法、时频域结合法；
所述三种方法即是求解三种不同代价函数的最优化问题；采用频域最优逼近法时，代价函数J设置为采用时域最优逼近法时，代价函数J设置为采用时频域结合法时，代价函数I设置为
其中，在上述三个式子中，k、l、n为参变量；Y(k)、X(k)、H(k)分别为接收信号、发送信号和超声探头冲激响应的频域表达式；Y(k)，k＝1，2，...，2N，是通过对y(n)，n＝1，2，...，N进行2N点的快速傅里叶变换得到的，N为所取信号长度；X(k)，k＝1，2，...，2N，是通过对x(n)，n＝1，2，...，N进行2N点的快速傅里叶变换得到的；参数k1、k2为在发送信号和接收信号的傅里叶变换谱上截取发送信号带宽所在的一段频谱，频谱的起始点和截止点，起始点k1＝int(f1N/fs)，截止点k2＝int(f2N/fs)，函数int(x)表示返回x的整数部分，f1、f2分别为发送信号的起始频率和截止频率，fs为采样频率；e为自然指数；为频域的频率间隔；L为h(n)的阶数。
进一步地，所述频域最优逼近法的实现方法是：令
使上式取得最小值的H(k)，即为发射探头和接收探头总的冲激响应h(n)的傅里叶变换谱，其中Jr、Ji分别为J的实部和虚部，具体写为：
其中，Yr(k)、Xr(k)、Hr(k)分别为接收信号、发送信号和超声探头冲激响应的频域表达式在第k点的实部，Yi(k)、Xi(k)、Hi(k)分别为接收信号、发送信号和超声探头冲激响应的频域表达式在第k点的虚部，下标r和i分别表示实部和虚部。
构造矩阵
其中，
此时有XXTH＝XY，解出矩阵H，并对其作快速反傅里叶变换，得到两个超声探头的总冲激响应h(n)。
进一步地，所述的时域最优逼近法的实现方法为：信号接收单元接收到的数字信号写为令
其中L为h(n)的阶数，该值越大，最后所得结果越精确，但会增加相应的计算量。求使上式取得最小值的h(n)，即为两个超声探头的总冲激响应；令 将上式写为矩阵形式为T＝Y-XTH；
h1～hL为信道冲激响应序列，x1～xN为发送信号的样本，y1～yN+L-1为接收信号的样本；
上式稍做调整并两边左乘X，得
XXTH＝X(Y-T)＝XY-XT；
由得XT＝0，于是上式写为XXTH＝XY，求解该矩阵方程得到矩阵H，然后得到两个超声探头总的冲激响应h(n)。
进一步地，所述的时频域结合法为：给定
其中则令上式取得最小值的h(n)即为两个超声探头总的冲激响应，于是
其中，
Jr、Ji分别为J的实部和虚部；Yr(k)、Xr(k)分别为接收信号、发送信号的频域表达式在第k点的实部，Yi(k)、Xi(k)分别为接收信号、发送信号的频域表达式在第k点的虚部，下标r和i分别表示实部和虚部；
为将目标代价函数I写为矩阵形式，构造矩阵其中
同理，构造矩阵其中
写出h(n)的矩阵形式，h0～hL-1为信道冲激响应序列，
根据最小二乘法，得到关系式H＝(XXT)-1XY，求解H，然后求得两个超声探头总的冲激响应h(n)。
进一步地，所述回波信号探头频响补偿的步骤如下：
S1.2.1、求发射探头和接收探头的总冲激响应h(n)的逆系统h-1(n)，两者的关系是：h-1(n)与h(n)的卷积为冲激函数，即满足h-1(n)*h(n)＝δ(n)；
S1.2.2、在数据处理单元中，使回波信号经过冲激响应为h-1(n)的均衡器，完成对超声探头的频响补偿。
进一步地，所述基于分数阶傅里叶变换的脉冲压缩与成像是利用线性调频信号在分数阶傅里叶域中的能量聚集特性进行脉冲压缩并成像，包括以下步骤：
S2.1、根据发送信号的参数确定分数阶傅里叶变换的最佳阶数，对经过补偿后的回波信号作最佳阶数的分数阶傅里叶变换，将回波信号转换到能量聚集域中，使其呈现冲激形式，实现脉冲压缩；
S2.2、对脉冲压缩后的结果进行包络检测；
S2.3、将分数阶傅里叶变换域横坐标u转换为时延τ，并将幅值转换为灰度值；
横坐标u根据下式转换为时延τ：
符号Fα[x(t)]表示对信号x(t)作基于旋转角度α的分数阶傅里叶变换，f0为所选信号的初始频率，k为调频率，A为信号幅度，最佳旋转角度α0＝-arccotk，T为有限观测时间，均为定值；与u一一对应，均可从坐标轴中得到，唯一未知数为τ，将每一个u转换为对应的时延τ，此时每一个时延值对应一个幅值
将幅值的最大值设置为255，则幅值与灰度值的转换关系式为：
其中，G是幅值对应为时的灰度值；
S2.4、时延与灰度值具有一一对应的关系，此时以时延为横坐标，灰度值为纵坐标，即可完成成像。
基于探头补偿的分数阶傅里叶变换回波成像系统，包括数据处理单元、D/A转换器、功率放大器、发射超声探头、接收超声探头、低噪声放大器、A/D转换器；
其中，数据处理单元、D/A转换器、功率放大器、发射超声探头顺次相连，实现发送信号的功能；接收超声探头、低噪声放大器、A/D转换器和数据处理单元顺次相连，完成对回波信号的接收、检测和成像工作；
所述数据处理单元完成系统中的多项工作，包括：产生参数可任意调节的数字宽频LFM信号；根据探头频响进行超声探头的频响补偿；对已补偿的回波信号作最佳阶数的分数阶傅里叶变换；包络检测；快速准确的求出时延，实现基于分数阶傅里叶变换的回波检测与成像；
所述D/A转换器用于将数据处理单元产生的数字宽频LFM信号转换为模拟宽频LFM信号；所述功率放大器用于将模拟宽频LFM信号的功率放大；所述发射超声探头用于发送经过放大的模拟LFM信号；所述接收超声探头用于接收回波信号；所述低噪声放大器用于对回波信号进行放大；所述A/D转换器用于将回波信号转换为数字信号再发送到数据处理单元。
进一步地，基于探头补偿的分数阶傅里叶变换回波成像系统的工作过程包括以下步骤：
步骤1、数据处理单元生成数字宽频LFM信号输入至D/A转换器，由D/A转换器转换为模拟信号，经过功率放大器放大后送至发射超声探头，然后由发射超声探头发出；
步骤2、人体组织反射回模拟回波信号，由接收装置中的接收超声探头接收，回波信号经过低噪声放大器放大后再经过A/D转换器完成模数转换，最后被送至数据处理单元；
步骤3、数据处理单元接收回波信号并进行探头频响补偿；
步骤4、对经过补偿后的回波信号作基于分数阶傅里叶变换的脉冲压缩与成像；
其中，基于分数阶傅里叶变换的脉冲压缩与成像的步骤如下：
步骤4.1、脉冲压缩；根据发送信号的参数确定分数阶傅里叶变换的最佳阶数，对经过补偿后的回波信号作最佳阶数的分数阶傅里叶变换，将回波信号转换到能量聚集域中，使其呈现冲激形式，实现脉冲压缩；
步骤4.2、包络检测；完成脉冲压缩后的回波信号，在数据处理单元中进行基于希尔伯特变换的包络检测，以消除旁瓣对成像的影响；
步骤4.3、坐标转换；数据处理单元根据下式完成横坐标u与时延τ的转换：
符号Fα[x(t)]表示对信号x(t)作基于旋转角度α的分数阶傅里叶变换，f0为所选信号的初始频率，k为调频率，A为信号幅度，最佳旋转角度α0＝-arccotk，T为有限观测时间，均为定值，与u一一对应，均可从坐标轴中得到，唯一未知数为τ，将每一个u转换为对应的时延τ；
将幅值的最大值设置为255，则幅值与灰度值的转换关系式为：
其中，G是幅值对应为时的灰度值；
步骤4.4、成像；数据处理单元根据时延与灰度值一一对应的关系，以时延为横坐标，灰度值为纵坐标，完成成像。
本发明与现有成像技术相比，具有如下优点：
(1)应用宽频的LFM信号，时宽-带宽积较大，有利于提高探测与成像精度。
(2)实现了对探头频响的自适应补偿，消除了探头频响对回波检测估计精度的影响。
(3)基于分数阶傅里叶变换的脉冲压缩算法精度高于匹配滤波法，回波成像精度高于传统方法。
(4)抗噪性能好，LFM信号的能量聚集特性可有效地分离信号中的噪声。
附图说明
图1是本发明实施例的进行超声探头频响测定的装置结构图；
图2是本发明实施例的基于探头补偿的分数阶傅里叶变换回波检测与成像的系统框图；
图3是本发明实施例的基于探头补偿的分数阶傅里叶变换回波检测与成像系统的工作流程图
图4是本发明实施例的基于分数阶傅里叶变换的脉冲压缩与成像流程图。
具体实施方式
下面结合实施例及附图对本发明的具体实施作进一步详细的描述，但本发明的实施方式不限于此。
基于探头补偿的分数阶傅里叶变换回波成像方法，包括以下步骤：
S1、先进行超声探头的频响补偿，包括超声探头频响测定与回波信号探头频响补偿两个部分；
S2、再进行基于分数阶傅里叶变换的脉冲压缩与成像。
如图1所示，进行超声探头频响测定的装置包括数据处理单元中的信号发送单元、D/A转换器、A/D转换器、发射超声探头、接收超声探头、数据处理单元中的信号接收单元。
所述的超声探头频响测定步骤如下：
S1.1.1、将发射探头与接收探头通过耦合剂耦合，两个超声探头总的冲激响应用h(n)表示；
S1.1.2、利用信号发送单元发送一个已知宽带线性调频信号即发送信号x(n)，由D/A转换器完成数模转换，再通过发射探头转换为声信号，该声信号经过耦合剂后由接收探头接收并转换为电信号，再由A/D转换器完成模数转换得到接收信号y(n)，将其发送给信号接收单元；
S1.1.3根据接收信号y(n)和发送信号x(n)求解超声探头的冲激响应h(n)。
步骤S1.1.3中，本实施例采用时频域结合法来求解两个超声探头总的冲激响应h(n)，具体如下：
对信号接收单元接收到的数字信号为y(n)，n＝1，2，...，N进行2N点的快速傅里叶变换，得到接收信号的傅里叶变换谱Y(k)，k＝1，2，...，2N，N为所取信号长度。对发射信号x(n)，n＝1，2，...，N进行2N点的快速傅里叶变换，得到发送信号的傅里叶变换谱X(k)，k＝1，2，...，2N。在发送信号和接收信号的傅里叶变换谱上截取发送信号带宽所在的一段频谱，频谱的起始点为k1＝int(f1N/fs)，截止点为k2＝int(f2N/fs)，其中函数int(x)表示返回x的整数部分，f1、f2分别为发送信号的起始频率和截止频率，fs为采样频率。给定
其中则令上式取得最小值的h(n)即为两个超声探头总的冲激响应，于是
其中，
Jr、Ji分别为J的实部和虚部；Yr(k)、Xr(k)分别为接收信号、发送信号的频域表达式在第k点的实部，Yi(k)、Xi(k)分别为接收信号、发送信号的频域表达式在第k点的虚部，下标r和i分别表示实部和虚部；
为将目标代价函数I写为矩阵形式，构造矩阵其中
同理，构造矩阵其中
写出h(n)的矩阵形式，h0～hL-1为信道冲激响应序列，
根据最小二乘法，得到关系式H＝(XXT)-1XY，求解H，然后求得两个超声探头总的冲激响应h(n)。
进一步地，所述回波信号探头频响补偿的步骤如下：
S1.2.1、求发射探头和接收探头的总冲激响应h(n)的逆系统h-1(n)，两者的关系是：h-1(n)与h(n)的卷积为冲激函数，即满足h-1(n)*h(n)＝δ(n)；
S1.2.2、在数据处理单元中，使回波信号经过冲激响应为h-1(n)的均衡器，完成对超声探头的频响补偿。
进一步地，所述基于分数阶傅里叶变换的脉冲压缩与成像是利用线性调频信号在分数阶傅里叶域中的能量聚集特性进行脉冲压缩并成像，包括以下步骤：
S2.1、根据发送信号的参数确定分数阶傅里叶变换的最佳阶数，对经过补偿后的回波信号作最佳阶数的分数阶傅里叶变换，将回波信号转换到能量聚集域中，使其呈现冲激形式，实现脉冲压缩；
S2.2、对脉冲压缩后的结果进行包络检测；
S2.3、将分数阶傅里叶变换域横坐标u转换为时延τ，并将幅值转换为灰度值；
横坐标u根据下式转换为时延τ：
符号Fα[x(t)]表示对信号x(t)作基于旋转角度α的分数阶傅里叶变换，f0为所选信号的初始频率，k为调频率，A为信号幅度，最佳旋转角度α0＝-arccotk，T为有限观测时间，均为定值；与u一一对应，均可从坐标轴中得到，唯一未知数为τ，将每一个u转换为对应的时延τ，此时每一个时延值对应一个幅值
将幅值的最大值设置为255，则幅值与灰度值的转换关系式为：
其中，G是幅值对应为时的灰度值；
S2.4、时延与灰度值具有一一对应的关系，此时以时延为横坐标，灰度值为纵坐标，即可完成成像。
如图2所示，基于探头补偿的分数阶傅里叶变换回波成像系统，包括数据处理单元、D/A转换器、功率放大器、发射超声探头、接收超声探头、低噪声放大器、A/D转换器。
其中，数据处理单元、D/A转换器、功率放大器、发射超声探头顺次相连，实现发送信号的功能；接收超声探头、低噪声放大器、A/D转换器和数据处理单元顺次相连，完成对回波信号的接收、检测和成像工作。
所述数据处理单元完成系统中的多项工作，包括：产生参数可任意调节的数字宽频LFM信号；根据探头频响进行超声探头的频响补偿；对已补偿的回波信号作最佳阶数的分数阶傅里叶变换；包络检测；快速准确的求出时延，实现基于分数阶傅里叶变换的回波检测与成像。
所述D/A转换器用于将数据处理单元产生的数字宽频LFM信号转换为模拟宽频LFM信号；所述功率放大器用于将模拟宽频LFM信号的功率放大；所述发射超声探头用于发送经过放大的模拟LFM信号；所述接收超声探头用于接收回波信号；所述低噪声放大器用于对回波信号进行放大；所述A/D转换器用于将回波信号转换为数字信号再发送到数据处理单元。
该系统发送信号为宽频LFM信号，本实施例中使用的是模拟频率从f0＝4MHz至f1＝6MHz、带宽B＝2MHz、持续时间为T＝20ms的LFM信号。数据处理单元处理的是离散的数字信号，因此需要对回波信号进行离散分数阶傅里叶变换，本实施例采用Ozaktas提出的采样型离散分数阶傅里叶变换。该算法以量纲归一化理论为基础建立算法模型，参量均为量纲归一化量(无量纲量)，以符号表示量纲归一化，为量纲归一化因子，分别为量纲归一化的LFM信号的调频率、持续时间和起始频率，能量聚集域的旋转角度在量纲归一化后，时域上的采样间隔Ts和u域的采样间隔一致，为Ts＝Us＝1/X，令n＝0，1，...，N-1表示采样点的序号，进行量纲归一化时将信号左移(N-1)/2，不考虑噪声，则经过时延τ后并经过补偿的回波信号表达式为：
其中，A为信号幅度，n0表示归一化后接收的回波信号的起始采样点序号，满足
为归一化时延，表示LFM信号的持续时间对应的样点数。对x(n)做最佳阶数的离散分数阶傅里叶变换并取绝对值，有
经过包络检测后，上式中的与u一一对应，均可从坐标轴中得到，由公式(2)可计算出每一个u对应的n0，再由公式(1)，即可计算出对应的时延τ，完成u与τ的转换。并根据以下规则将幅值转换为灰度值：将幅值的最大值设置为255，则转换关系式为：
其中，G是幅值对应为时的灰度值。根据时延与灰度值一一对应的关系，以时延为横坐标，灰度值为纵坐标，即可完成成像。
如图3所示，基于探头补偿的分数阶傅里叶变换回波成像系统的工作过程包括以下步骤：
步骤1、数据处理单元生成数字宽频LFM信号输入至D/A转换器，由D/A转换器转换为模拟信号，经过功率放大器放大后送至发射超声探头，然后由发射超声探头发出；
步骤2、人体组织反射回模拟回波信号，由接收装置中的接收超声探头接收，回波信号经过低噪声放大器放大后再经过A/D转换器完成模数转换，最后被送至数据处理单元；
步骤3、数据处理单元接收回波信号并进行探头频响补偿；
步骤4、对经过补偿后的回波信号作基于分数阶傅里叶变换的脉冲压缩与成像。
如图4所示，基于分数阶傅里叶变换的脉冲压缩与成像的步骤如下：
步骤1：脉冲压缩。根据发送信号的参数确定分数阶傅里叶变换的最佳阶数：对经过补偿后的回波信号作p0阶的分数阶傅里叶变换，将回波信号转换到能量聚集域中，使其呈现冲激形式，实现脉冲压缩。
步骤2：包络检测。完成脉冲压缩后的回波信号，将在数据处理单元中进行基于希尔伯特变换的包络检测，以消除旁瓣对成像的影响。
步骤3：坐标转换。数据处理单元根据公式(1)、公式(2)完成横坐标u与时延τ的转换，并根据公式(3)将幅值转换为灰度值。
步骤4：成像。数据处理单元根据时延与灰度值一一对应的关系，以时延为横坐标，灰度值为纵坐标，完成成像。人体组织对信号的反射越强烈，灰度值越大，所成像就越亮，反之越暗。
Claims (10)
1.基于探头补偿的分数阶傅里叶变换回波成像方法，其特征在于，包括以下步骤：
S1、先进行超声探头的频响补偿，包括超声探头频响测定与回波信号探头频响补偿两个部分；
S2、再进行基于分数阶傅里叶变换的脉冲压缩与成像。
2.根据权利要求1所述的基于探头补偿的分数阶傅里叶变换回波成像方法，其特征在于，所述超声探头频响测定的步骤如下：
S1.1.1、将发射探头与接收探头通过耦合剂耦合，两个超声探头总的冲激响应用h(n)表示；
S1.1.2、利用数据处理单元中的信号发送单元发送一个已知宽带线性调频信号即发送信号x(n)，由D/A转换器完成数模转换，再通过发射探头转换为声信号，该声信号经过耦合剂后由接收探头接收并转换为电信号，再由A/D转换器完成模数转换得到接收信号y(n)，将其发送给数据处理单元中的信号接收单元；
S1.1.3、根据接收信号y(n)和发送信号x(n)求解超声探头的冲激响应h(n)。
3.根据权利要求2所述的基于探头补偿的分数阶傅里叶变换回波成像方法，其特征在于，求取探头冲激响应的方法可采用频域最优逼近法、时域最优逼近法、时频域结合法；
所述三种方法即是求解三种不同代价函数的最优化问题；采用频域最优逼近法时，代价函数J设置为采用时域最优逼近法时，代价函数J设置为采用时频域结合法时，代价函数J设置为
其中，在上述三个式子中，k、l、n为参变量；Y(k)、X(k)、H(k)分别为接收信号、发送信号和超声探头冲激响应的频域表达式；Y(k)，k＝1，2，...，2N，是通过对y(n)，n＝1，2，...，N进行2N点的快速傅里叶变换得到的，N为所取信号长度；X(k)，k＝1，2，...，2N，是通过对x(n)，n＝1，2，...，N进行2N点的快速傅里叶变换得到的；参数k1、k2为在发送信号和接收信号的傅里叶变换谱上截取发送信号带宽所在的一段频谱，频谱的起始点和截止点，起始点k1＝int(f1N/fs)，截止点k2＝int(f2N/fs)，函数int(x)表示返回x的整数部分，f1、f2分别为发送信号的起始频率和截止频率，fs为采样频率；e为自然指数；为频域的频率间隔；L为h(n)的阶数。
4.根据权利要求3所述的基于探头补偿的分数阶傅里叶变换回波成像方法，其特征在于，所述频域最优逼近法的实现方法是：令
使上式取得最小值的H(k)，即为发射探头和接收探头总的冲激响应h(n)的傅里叶变换谱，其中Jr、Ji分别为J的实部和虚部，具体写为：
其中，Yr(k)、Xr(k)、Hr(k)分别为接收信号、发送信号和超声探头冲激响应的频域表达式在第k点的实部，Yi(k)、Xi(k)、Hi(k)分别为接收信号、发送信号和超声探头冲激响应的频域表达式在第k点的虚部，下标r和i分别表示实部和虚部；
构造矩阵
其中，
此时有XXTH＝XY，解出矩阵H，并对其作快速反傅里叶变换，得到两个超声探头的总冲激响应h(n)。
5.根据权利要求3所述的基于探头补偿的分数阶傅里叶变换回波成像方法，其特征在于，所述的时域最优逼近法的实现方法为：信号接收单元接收到的数字信号写为 令
其中L为h(n)的阶数，求使上式取得最小值的h(n)，即为两个超声探头的总冲激响应；令将上式写为矩阵形式为
T＝Y-XTH；
h1～hL为信道冲激响应序列，x1～xN为发送信号的样本，y1～yN+L-1为接收信号的样本；
上式稍做调整并两边左乘X，得
XXTH＝X(Y-T)＝XY-XT；
由得XT＝0，于是上式写为XXTH＝XY，求解该矩阵方程得到矩阵H，然后得到两个超声探头总的冲激响应h(n)。
6.根据权利要求3所述的基于探头补偿的分数阶傅里叶变换回波成像方法，其特征在于，所述的时频域结合法为：给定
其中则令上式取得最小值的h(n)即为两个超声探头总的冲激响应，于是
其中，
Jr、Ji分别为J的实部和虚部；Yr(k)、Xr(k)分别为接收信号、发送信号的频域表达式在第k点的实部，Yi(k)、Xi(k)分别为接收信号、发送信号的频域表达式在第k点的虚部，下标r和i分别表示实部和虚部；
为将目标代价函数J写为矩阵形式，构造矩阵其中
同理，构造矩阵其中
写出h(n)的矩阵形式，h0～hL-1为信道冲激响应序列，
根据最小二乘法，得到关系式H＝(XXT)-1XY，求解H，然后求得两个超声探头总的冲激响应h(n)。
7.根据权利要求1所述的基于探头补偿的分数阶傅里叶变换回波成像方法，其特征在于，所述回波信号探头频响补偿的步骤如下：
S1.2.1、求发射探头和接收探头的总冲激响应h(n)的逆系统h-1(n)，两者的关系是：h-1(n)与h(n)的卷积为冲激函数，即满足h-1(n)*h(n)＝δ(n)；
S1.2.2、在数据处理单元中，使回波信号经过冲激响应为h-1(n)的均衡器，完成对超声探头的频响补偿。
8.根据权利要求1所述的基于探头补偿的分数阶傅里叶变换回波成像方法，其特征在于，所述基于分数阶傅里叶变换的脉冲压缩与成像是利用线性调频信号在分数阶傅里叶域中的能量聚集特性进行脉冲压缩并成像，包括以下步骤：
S2.1、根据发送信号的参数确定分数阶傅里叶变换的最佳阶数，对经过补偿后的回波信号作最佳阶数的分数阶傅里叶变换，将回波信号转换到能量聚集域中，使其呈现冲激形式，实现脉冲压缩；
S2.2、对脉冲压缩后的结果进行包络检测；
S2.3、将分数阶傅里叶变换域横坐标u转换为时延τ，并将幅值转换为灰度值；
横坐标u根据下式转换为时延τ：
符号Fα[x(t)]表示对信号x(t)作基于旋转角度α的分数阶傅里叶变换，f0为所选信号的初始频率，k为调频率，A为信号幅度，最佳旋转角度α0＝-arccotk，T为有限观测时间；
将幅值的最大值设置为255，则幅值与灰度值的转换关系式为：
其中，G是幅值对应为时的灰度值；
S2.4、时延与灰度值具有一一对应的关系，此时以时延为横坐标，灰度值为纵坐标，即可完成成像。
9.基于探头补偿的分数阶傅里叶变换回波成像系统，其特征在于，包括数据处理单元、D/A转换器、功率放大器、发射超声探头、接收超声探头、低噪声放大器、A/D转换器；
其中，数据处理单元、D/A转换器、功率放大器、发射超声探头顺次相连，实现发送信号的功能；接收超声探头、低噪声放大器、A/D转换器和数据处理单元顺次相连，完成对回波信号的接收、检测和成像工作；
所述数据处理单元完成系统中的多项工作，包括：产生参数可任意调节的数字宽频LFM信号；根据探头频响进行超声探头的频响补偿；对已补偿的回波信号作最佳阶数的分数阶傅里叶变换；包络检测；快速准确的求出时延，实现基于分数阶傅里叶变换的回波检测与成像；
所述D/A转换器用于将数据处理单元产生的数字宽频LFM信号转换为模拟宽频LFM信号；所述功率放大器用于将模拟宽频LFM信号的功率放大；所述发射超声探头用于发送经过放大的模拟LFM信号；所述接收超声探头用于接收回波信号；所述低噪声放大器用于对回波信号进行放大；所述A/D转换器用于将回波信号转换为数字信号再发送到数据处理单元。
10.根据权利要求9所述的基于探头补偿的分数阶傅里叶变换回波成像系统，其特征在于，基于探头补偿的分数阶傅里叶变换回波成像系统的工作过程包括以下步骤：
步骤1、数据处理单元生成数字宽频LFM信号输入至D/A转换器，由D/A转换器转换为模拟信号，经过功率放大器放大后送至发射超声探头，然后由发射超声探头发出；
步骤2、人体组织反射回模拟回波信号，由接收装置中的接收超声探头接收，回波信号经过低噪声放大器放大后再经过A/D转换器完成模数转换，最后被送至数据处理单元；
步骤3、数据处理单元接收回波信号并进行探头频响补偿；
步骤4、对经过补偿后的回波信号作基于分数阶傅里叶变换的脉冲压缩与成像；
其中，基于分数阶傅里叶变换的脉冲压缩与成像的步骤如下：
步骤4.1、脉冲压缩；根据发送信号的参数确定分数阶傅里叶变换的最佳阶数，对经过补偿后的回波信号作最佳阶数的分数阶傅里叶变换，将回波信号转换到能量聚集域中，使其呈现冲激形式，实现脉冲压缩；
步骤4.2、包络检测；完成脉冲压缩后的回波信号，在数据处理单元中进行基于希尔伯特变换的包络检测，以消除旁瓣对成像的影响；
步骤4.3、坐标转换；数据处理单元根据下式完成横坐标u与时延τ的转换：
符号Fα[x(t)]表示对信号x(t)作基于旋转角度α的分数阶傅里叶变换，f0为所选信号的初始频率，k为调频率，A为信号幅度，最佳旋转角度α0＝-arccotk，T为有限观测时间；
将幅值的最大值设置为255，则幅值与灰度值的转换关系式为：
其中，G是幅值对应为时的灰度值；
步骤4.4、成像；数据处理单元根据时延与灰度值一一对应的关系，以时延为横坐标，灰度值为纵坐标，完成成像。
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910803651.7A CN110559014B (zh) | 2019-08-28 | 2019-08-28 | 基于探头补偿的分数阶傅里叶变换回波成像方法与系统 |
| PCT/CN2020/112295 WO2021037252A1 (zh) | 2019-08-28 | 2020-08-28 | 基于探头补偿的分数阶傅里叶变换回波成像方法与系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910803651.7A CN110559014B (zh) | 2019-08-28 | 2019-08-28 | 基于探头补偿的分数阶傅里叶变换回波成像方法与系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110559014A true CN110559014A (zh) | 2019-12-13 |
| CN110559014B CN110559014B (zh) | 2021-10-26 |
Family
ID=68776579
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910803651.7A Expired - Fee Related CN110559014B (zh) | 2019-08-28 | 2019-08-28 | 基于探头补偿的分数阶傅里叶变换回波成像方法与系统 |
Country Status (2)
| Country | Link |
|---|---|
| CN (1) | CN110559014B (zh) |
| WO (1) | WO2021037252A1 (zh) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111012317A (zh) * | 2020-01-18 | 2020-04-17 | 四川知周光声医疗科技有限公司 | 一种光声乳腺的图像重建方法及系统 |
| CN111796259A (zh) * | 2020-06-08 | 2020-10-20 | 华南理工大学 | 一种基于lfm互相关模特性逼近的子样本时延估计方法 |
| CN111948621A (zh) * | 2020-07-14 | 2020-11-17 | 浙江大学 | 一种极低信噪比下最佳稀疏域变换的lfm雷达信号处理方法 |
| WO2021037252A1 (zh) * | 2019-08-28 | 2021-03-04 | 华南理工大学 | 基于探头补偿的分数阶傅里叶变换回波成像方法与系统 |
| CN113687364A (zh) * | 2021-07-13 | 2021-11-23 | 中国船舶重工集团公司第七一五研究所 | 一种三维成像声纳图像假底的抑制方法 |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113030861B (zh) * | 2021-03-07 | 2023-06-23 | 西北工业大学 | 一种时间反转和分数阶小波变换的水声信号检测方法 |
| CN113189543B (zh) * | 2021-04-27 | 2023-07-14 | 哈尔滨工程大学 | 一种基于运动补偿鲁棒主成分分析的干扰抑制方法 |
| CN116524051B (zh) * | 2023-04-10 | 2024-01-09 | 哈尔滨工业大学 | 一种基于分数阶傅里叶变换域模态分解的高分辨isar成像方法 |
| CN116633323B (zh) * | 2023-04-25 | 2024-04-16 | 中国计量科学研究院 | 基于光电导技术高速数字采集系统响应特性校准方法及系统 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH1033529A (ja) * | 1996-07-22 | 1998-02-10 | Ge Yokogawa Medical Syst Ltd | 超音波撮像方法および装置 |
| CN1815552A (zh) * | 2006-02-28 | 2006-08-09 | 安徽中科大讯飞信息科技有限公司 | 基于线谱频率及其阶间差分参数的频谱建模与语音增强方法 |
| US20150348531A1 (en) * | 2012-12-19 | 2015-12-03 | University Of Leeds | Ultrasound generation |
| CN105919624A (zh) * | 2016-05-25 | 2016-09-07 | 飞依诺科技(苏州)有限公司 | 一种高脉冲重复扫查频率信号的编码、成像方法和装置 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105997142A (zh) * | 2016-06-12 | 2016-10-12 | 飞依诺科技(苏州)有限公司 | 一种超声系统发射信号的复合、成像方法及装置 |
| CN110559014B (zh) * | 2019-08-28 | 2021-10-26 | 华南理工大学 | 基于探头补偿的分数阶傅里叶变换回波成像方法与系统 |
-
2019
- 2019-08-28 CN CN201910803651.7A patent/CN110559014B/zh not_active Expired - Fee Related
-
2020
- 2020-08-28 WO PCT/CN2020/112295 patent/WO2021037252A1/zh active Application Filing
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH1033529A (ja) * | 1996-07-22 | 1998-02-10 | Ge Yokogawa Medical Syst Ltd | 超音波撮像方法および装置 |
| CN1815552A (zh) * | 2006-02-28 | 2006-08-09 | 安徽中科大讯飞信息科技有限公司 | 基于线谱频率及其阶间差分参数的频谱建模与语音增强方法 |
| US20150348531A1 (en) * | 2012-12-19 | 2015-12-03 | University Of Leeds | Ultrasound generation |
| CN105919624A (zh) * | 2016-05-25 | 2016-09-07 | 飞依诺科技(苏州)有限公司 | 一种高脉冲重复扫查频率信号的编码、成像方法和装置 |
Non-Patent Citations (3)
| Title |
|---|
| 李进等: "基于分数阶Fourier变换的线性调频超声回波信号的滤波", 《生物医学工程研究》 * |
| 杨娜等: "《传感器与测试技术》", 31 July 2012, 航空工业出版社 * |
| 陈敏: "声频定向系统理论与关键技术研究", 《中国博士学位论文全文数据库 工程科技Ⅱ辑》 * |
Cited By (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2021037252A1 (zh) * | 2019-08-28 | 2021-03-04 | 华南理工大学 | 基于探头补偿的分数阶傅里叶变换回波成像方法与系统 |
| CN111012317A (zh) * | 2020-01-18 | 2020-04-17 | 四川知周光声医疗科技有限公司 | 一种光声乳腺的图像重建方法及系统 |
| CN111012317B (zh) * | 2020-01-18 | 2022-10-25 | 中川新迈科技有限公司 | 一种光声乳腺的图像重建方法及系统 |
| CN111796259A (zh) * | 2020-06-08 | 2020-10-20 | 华南理工大学 | 一种基于lfm互相关模特性逼近的子样本时延估计方法 |
| CN111796259B (zh) * | 2020-06-08 | 2023-04-07 | 华南理工大学 | 一种基于lfm互相关模特性逼近的子样本时延估计方法 |
| CN111948621A (zh) * | 2020-07-14 | 2020-11-17 | 浙江大学 | 一种极低信噪比下最佳稀疏域变换的lfm雷达信号处理方法 |
| CN111948621B (zh) * | 2020-07-14 | 2023-06-27 | 浙江大学 | 一种极低信噪比下最佳稀疏域变换的lfm雷达信号处理方法 |
| CN113687364A (zh) * | 2021-07-13 | 2021-11-23 | 中国船舶重工集团公司第七一五研究所 | 一种三维成像声纳图像假底的抑制方法 |
| CN113687364B (zh) * | 2021-07-13 | 2024-03-19 | 中国船舶重工集团公司第七一五研究所 | 一种三维成像声纳图像假底的抑制方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| WO2021037252A1 (zh) | 2021-03-04 |
| CN110559014B (zh) | 2021-10-26 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110559014B (zh) | 基于探头补偿的分数阶傅里叶变换回波成像方法与系统 | |
| US20070291591A1 (en) | High-resolution ultrasound displacement measurement apparatus and method | |
| CN107714091B (zh) | 经颅低频超声线性调频脉冲逆转微泡成像方法 | |
| CN103490754B (zh) | 一种大时间带宽积的超声信号及其脉冲压缩方法及系统 | |
| Yücel et al. | Coded waveform excitation for high-resolution ultrasonic guided wave response | |
| Ramalli et al. | A real-time chirp-coded imaging system with tissue attenuation compensation | |
| Chernyakova et al. | Fourier-domain beamforming and structure-based reconstruction for plane-wave imaging | |
| CN107970042B (zh) | 一种超声非线性定量实时成像方法与系统 | |
| Liu et al. | Compressed sensing based synthetic transmit aperture imaging: Validation in a convex array configuration | |
| Tang et al. | A method based on SVD for detecting the defect using the magnetostrictive guided wave technique | |
| JPS62249058A (ja) | 超音波検出方法及び装置 | |
| Tang et al. | Non-contact phase coded excitation of ultrasonic Lamb wave for blind hole inspection | |
| CN1571649A (zh) | 用于指示身体组织机械硬度特性的设备和方法 | |
| Sheng et al. | Application of pulse compression technology in electromagnetic ultrasonic thickness measurement | |
| Fan et al. | Coded excitation with nonlinear frequency modulation carrier in ultrasound imaging system | |
| GB2436655A (en) | Deformation determination in pulse echo or ultrasonic imaging | |
| Guo et al. | Reconstructing outside pass-band data to improve time resolution in ultrasonic detection | |
| CN113822363B (zh) | 一种基于二次有理核函数的卷积神经网络时延估计方法 | |
| Yücel et al. | Pulse-compression based iterative time-of-flight extraction of dispersed ultrasonic guided waves | |
| JPH11309146A (ja) | 超音波散乱体の流れをイメ―ジングするシステムおよび方法 | |
| CN111812622A (zh) | 一种基于透镜回波的超声探头带宽检测方法及系统 | |
| CN103565469B (zh) | 一种提高二维阵列实时超声体成像质量的方法 | |
| CN113050043A (zh) | 基于非理想lpf的探地雷达超宽带高斯脉冲fri采样方法 | |
| Peng et al. | Chirp-coded pulse excitation for ultrasound elasticity imaging | |
| LIU et al. | Simulation of the Application of Compressive Sensing to Ultrasound Total Focusing Method Imaging |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20211026 |