CN111856401A - 一种基于互谱相位拟合的时延估计方法 - Google Patents

一种基于互谱相位拟合的时延估计方法 Download PDF

Info

Publication number
CN111856401A
CN111856401A CN202010633065.5A CN202010633065A CN111856401A CN 111856401 A CN111856401 A CN 111856401A CN 202010633065 A CN202010633065 A CN 202010633065A CN 111856401 A CN111856401 A CN 111856401A
Authority
CN
China
Prior art keywords
time delay
cross
frequency
spectrum
delay estimation
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.)
Withdrawn
Application number
CN202010633065.5A
Other languages
English (en)
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.)
Nanjing University
Original Assignee
Nanjing University
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 Nanjing University filed Critical Nanjing University
Priority to CN202010633065.5A priority Critical patent/CN111856401A/zh
Publication of CN111856401A publication Critical patent/CN111856401A/zh
Withdrawn legal-status Critical Current

Links

Images

Classifications

    • 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
    • G01S5/00—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
    • G01S5/18—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using ultrasonic, sonic, or infrasonic waves
    • G01S5/22—Position of source determined by co-ordinating a plurality of position lines defined by path-difference measurements

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Circuit For Audible Band Transducer (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

本发明公开了一种基于互谱相位拟合的时延估计方法，包括以下步骤：利用麦克风阵列采集空间声源信号；通过传统GCC时延估计算法计算每两组阵元接收信号之间的时延整数部分；重新抽样信号，计算互功率谱；对互谱进行频域滤波，利用互谱相位跟频率的线性关系计算每个频点的时延值、拟合得到分数时延估计值，综合确定整体时延。本发明支持整数时延估计误差不大的情况下，实现对两路麦克风接收信号时延值较为准确的估计。

Description

一种基于互谱相位拟合的时延估计方法
技术领域
本发明属于阵列信号处理领域，具体涉及一种基于互谱相位拟合的时延估计方法。
背景技术
基于麦克风阵列时延估计算法的声源定位技术具有原理简单、实时性强等特点，因此应用广泛。该技术主要分为两步：第一步是根据麦克风阵列各阵元的接收信号，估计每两路接收信号的相对时延，目前主要有GCC、基于LMS自适应滤波器等时延估计算法；第二步结合声源与麦克风阵列的空间几何关系进行声源定位。实际声源定位的精度取决于时延估计精度，而由于环境噪声、系统采样等各种因素影响，时延估计误差会非常大。利用传统的时域滤波器会产生一定相位时延，而线性相位滤波器对系统有较高的要求，不便于实际应用；基于插值、多项式拟合等算法的分数阶时延估计存在精度不够高、计算量过于复杂等一系列问题。
发明内容
发明目的：针对上述现有技术存在的问题和不足，本发明的目的是提供一种基于互谱相位拟合的时延估计的方法，支持整数延迟误差不大的情况下，实现对两路信号相对时延较为准确的估计。
技术方案：为实现上述发明目的，本发明采用的技术方案为一种基于互谱相位拟合的时延估计方法，包括以下步骤：
1)利用麦克风阵列采集声源信号，计算每两路接收信号整数部分时延；
2)根据整数部分时延，重新抽样，构建延迟点数小于1的两路接收信号，求出互功率谱；
3)对所述互功率谱进行频域滤波，去除噪声、杂波干扰；
4)利用互谱相位跟频率的线性关系计算每个频点的时延，拟合得到分数部分时延估计值，综合得到整体时延。
进一步地，所述步骤1)包括以下步骤：
11)利用麦克风阵列采集声源信息，在同一时刻抽样长度为N的两路接收信号
Figure BDA0002566487730000021
12)根据传统GCC(Generalized Cross Correlation,广义互相关)时延估计算法计算相对时延整数部分m。
进一步地，所述步骤2)包括以下步骤：
21)信号
Figure BDA0002566487730000022
不变，将信号
Figure BDA0002566487730000023
抽样位置顺延m个采样单元，重新抽样长度N得到
Figure BDA0002566487730000024
22)对信号
Figure BDA0002566487730000025
分别进行FFT(Fast Fourier Transform,快速傅里叶变换)转换到频域，得到
Figure BDA0002566487730000026
再将
Figure BDA0002566487730000027
取共轭与
Figure BDA0002566487730000028
相乘得到互功率谱
Figure BDA0002566487730000029
进一步地，所述步骤3)包括以下步骤：
31)将声源信号频率范围[f1,f2]除以频率分辨率f0、取整，得到相应的正频率部分频点k值范围[k1,k2]；
32)设计一个长度为N的滤波器，系数
Figure BDA00025664877300000210
在
Figure BDA00025664877300000211
Figure BDA00025664877300000212
范围内的值取1，其余为0，而后与互谱进行点乘运算得到新互谱，实现频域滤波；也可以区分正负频率，直接将互谱
Figure BDA00025664877300000213
在频点范围[k1,k2]、[-k2,-k1]以外的值置0的方式达到滤波效果；
33)利用angle函数获得滤波后的互谱相位
Figure BDA00025664877300000214
N'为滤波后保留的互谱长度。
进一步地，所述步骤4)包括以下步骤：
41)根据互谱相位与频率的线性关系，将
Figure BDA00025664877300000215
除以
Figure BDA00025664877300000216
得到每个频点的时延值
Figure BDA00025664877300000217
42)对
Figure BDA00025664877300000218
进行最小二乘线性拟合，确定两路接收信号
Figure BDA00025664877300000219
Figure BDA00025664877300000220
分数部分时延估计值
Figure BDA00025664877300000221
43)与步骤1)得到的整数部分m求和得到整体时延值△n。
有益效果：本发明提出的基于互谱相位拟合的时延估计算法不仅具有原理简单、实时性强等优点，还解决了传统时域滤波器产生相位延时以及滤波不彻底的问题。同时，针对传统分数阶时延估计算法精度差、计算量大等问题，提出一种利用两路信号互谱相位拟合的新型时延估计方法；实验表明，整数时延估计精度误差不大的情况下，本发明都可以保证较高的时延估计精度。
附图说明
图1为本发明提供的基于互谱相位拟合的时延估计方法的场景结构图；
图2为本发明提供的基于互谱相位拟合的时延估计方法的处理流程图。
图3为本发明提供的基于互谱相位拟合的时延估计方法的频域滤波流程图。
具体实施方式
下面结合附图和具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。
本发明公开了一种基于互谱相位拟合的时延估计方法，适用于对麦克风阵列接收的声源信号进行相对时延估计，实现对声源测向或定位。通过频域滤波，可以理想地剔除噪声、杂波影响，获得纯净的有用信号，而且不会产生相位延时；传统的时延估计算法在时域进行细化处理，不可避免受到噪声等因素干扰，同时存在计算量大、精度不够高等问题。基于互谱相位拟合的时延估计算法在频域进行处理，可以完美地解决该问题，能够快速、准确的估计出相对时延值。
下面结合附图，详细描述本发明的技术方案：
本发明的应用场景如图1所示，以空间正四面体麦克风阵列模型为例，采集声源信号，从中提取数据信息。该系统的特征是麦克风阵列处于被动式接收信号，具有隐蔽性和安全性。我们在MATLAB环境中运行脚本。
在本发明的方案中，设采样频率Fs，抽样信号长度N，声源s(t)，第i个麦克风接收信号xi(t)为：
xi(t)＝αis(t-ti)+ni(t) (1)
其中，αi表示声源到达第i个阵元的幅度衰减，s(t-ti)表示声源s(t)经过ti时间到达麦克风i，ni(t)为第i个麦克风接收到的相互独立的高斯白噪声。
整个发明流程如图2所示，首先在同一时刻，抽取长度为N的两路麦克风接收信号
Figure BDA0002566487730000041
n表示信号x[n]的索引，利用传统GCC时延估计算法获得整数部分时延值。GCC时延估计算法主要流程是先将两路接收信号分别进行FFT转换到频域，得到
Figure BDA0002566487730000042
k表示X[k]的索引，再将
Figure BDA0002566487730000043
取共轭与
Figure BDA0002566487730000044
相乘得到互功率谱，而后在互谱上加权通过IFFT(Inverse Fast Fourier Transform，快速逆傅里叶变换)到时域进行峰值检测得到整数时延值m。即：
Figure BDA0002566487730000045
Figure BDA0002566487730000046
其中，“*”表示共轭运算，
Figure BDA00025664877300000417
表示两路信号的互相关函数，其峰值对应的位置即两路信号的时延值，
Figure BDA0002566487730000047
表示频域加权函数，X1[k]表示x1[n]的快速傅里叶变换，X2 *(k-m)表示将信号x2[n]的快速傅里叶变换X2[k]取共轭，平移m个频域间隔。
根据整数部分时延，将信号
Figure BDA0002566487730000048
抽样位置顺延m个采样单元，得到
Figure BDA0002566487730000049
在整数时延准确的前提下，
Figure BDA00025664877300000410
与
Figure BDA00025664877300000411
的相对延迟点数在-1与1之间。将重新抽样后的两路信号加窗、FFT并写成指数形式：
Figure BDA00025664877300000412
其中，X(k)表示x[n]的快速傅里叶变换，|X(k)|表示幅度部分，
Figure BDA00025664877300000413
表示相位。
根据FFT的叠加性和时移特性，二者的互功率谱展开式为：
Figure BDA00025664877300000414
其中，S(k)、Ni(k)分别表示声源s(t)、噪声ni(t)的快速傅里叶变换，S*(k)、Ni *(k)为对应的共轭运算，exp表示以自然常数e为底的指数，p为待估计的分数部分时延值。
带入频率f表达式：
Figure BDA00025664877300000415
其中，f0为抽样信号频率分辨率。用W(k)代替
Figure BDA00025664877300000416
中含有噪声部分，即：
Figure BDA0002566487730000051
进一步地，式(6)简化为：
Figure BDA0002566487730000052
抽样信号确定，p为常数，如果不考虑噪声W(k)，则互谱相位
Figure BDA0002566487730000053
与频率呈线性关系，即：
Figure BDA0002566487730000054
且，每个频点的时延值应该相等。实际中，由于噪声、杂波干扰，会出现抖动误差。
如图3所示，考虑到声源信号与噪声频率差异，在频谱图上出现的位置不同，因此通过频域滤波的方式可以剔除噪声、其他波段杂波干扰。设声源信号频率范围为[f1,f2]，分别除以频率分辨率f0、取整，得到相应的k值范围为[k1,k2]。
设计一个长度为N的滤波器，其系数为序列
Figure BDA0002566487730000055
将
Figure BDA0002566487730000056
Figure BDA0002566487730000057
范围内的系数值h取1，其余为0，而后与互谱进行点乘运算得到新互谱，实现频域滤波；也可以区分正负频率，直接将互谱
Figure BDA0002566487730000058
在频点范围[k1,k2]、[-k2,-k1]以外的值置0，与频域滤波效果是一致的。得到的即只含有声源信号的新互谱，设滤波后的互谱长度为N'。即：
Figure BDA0002566487730000059
其中，fbandpass为信号频带范围。同时，对互谱直接进行频域滤波与对两路信号分别滤波求互谱，结果是等效的，且前者更便于计算。
进一步地，利用相位函数angle得到滤波后互谱
Figure BDA00025664877300000510
的相位，为了防止整数延迟点数误差，可以加相位解卷绕函数unwrap，而后除以
Figure BDA00025664877300000511
得到每个频点的时延值，即
Figure BDA00025664877300000512
Figure BDA00025664877300000513
如果每个频点的时延值仍不同，说明还有残留噪声或者杂波；考虑到互谱的对称性，可以对正频率部分
Figure BDA00025664877300000514
进行最小二乘线性拟合、取中值的方式确定信号
Figure BDA00025664877300000515
的时延估计值
Figure BDA00025664877300000516
即
Figure BDA00025664877300000517
的分数部分时延估计值。
综合整数时延m，求和得到两路接收信号
Figure BDA0002566487730000061
整体时延值△n：
Figure BDA0002566487730000062
利用图1所示的正四面体麦克风阵列模型对空间声源进行仿真实验，阵元间距d＝1m，采样频率Fs＝32kHz，抽样信号长度N＝1024，取一段狙击枪声信号作为声源信号，频率下限f1＝100Hz，频率上限f2＝1000Hz，模拟不同的时延值，结果见表1所示。由于每个频点只能取抽样信号频率分辨率f0的整数倍，实际滤波后保留的频率范围是：125:1000Hz。
表1 MATLAB仿真结果
Figure BDA0002566487730000063
仿真证明，本发明提出的基于互谱相位拟合的时延估计算法不仅原理简单、时延估计精度高，可以精确到0.01采样单元；还解决了时域滤波不彻底、附加相位延时的问题，得到了纯净的声源信号；对于不同频带的多声源，也可以通过频域滤波的方式进行区分，而后分别进行定位。同时，利用GCC获得的整数时延估计误差不大的情况下，本发明都可以准确估计出整体时延值；相应的，对于整数时延比较小的两路接收信号，可以省略GCC求整数时延部分、直接利用互谱相位拟合计算整体时延值，但是精度会有所下降。
本发明具体应用途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进，这些改进也应视为本发明的保护范围。

Claims (6)

1.一种基于互谱相位拟合的时延估计方法，其特征在于，包括以下步骤：
1)利用麦克风阵列采集声源信号，计算每两路接收信号相对时延的整数部分；
2)根据整数部分时延，重新抽样，构建延迟点数小于1的两路接收信号，求出互功率谱；
3)对所述互功率谱进行频域滤波，去除噪声、杂波干扰，计算滤波后的互谱相位；
4)利用互谱相位跟频率的线性关系计算每个频点的时延，拟合得到分数部分时延估计值，综合得到整体时延。
2.根据权利要求1所述一种基于互谱相位拟合的时延估计方法，其特征在于：所述步骤1)包括以下步骤：
11)利用麦克风阵列采集声源信息，在同一时刻抽样长度为N的两路接收信号
Figure FDA0002566487720000011
其中x1[n]表示第1个麦克风接收信号，n表示信号x[n]的索引，x2[n]表示第2个麦克风接收信号；
12)根据传统GCC时延估计算法计算相对时延整数部分m。
3.根据权利要求2所述一种基于互谱相位拟合的时延估计方法，其特征在于：所述步骤2)包括以下步骤：
21)信号
Figure FDA0002566487720000012
不变，将信号
Figure FDA0002566487720000013
抽样位置顺延m个采样单元，重新抽样长度N得到
Figure FDA0002566487720000014
22)对信号
Figure FDA0002566487720000015
分别进行FFT转换到频域，得到
Figure FDA0002566487720000016
Figure FDA0002566487720000017
k表示X[k]的索引，再将
Figure FDA0002566487720000018
取共轭与
Figure FDA0002566487720000019
相乘得到互功率谱
Figure FDA00025664877200000110
4.根据权利要求2所述一种基于互谱相位拟合的时延估计方法，其特征在于：所述步骤3)包括以下步骤：
31)将声源信号频率范围[f1,f2]除以频率分辨率f0、取整，得到相应的正频率部分频点k值范围[k1,k2]，其中f1表示频率下限，f2表示频率上限；
32)设计一个长度为N的滤波器，系数
Figure FDA0002566487720000021
在
Figure FDA0002566487720000022
Figure FDA0002566487720000023
范围内的值取1，其余为0，而后与互谱进行点乘运算得到新互谱，实现频域滤波；
33)利用angle函数获得滤波后的互谱相位
Figure FDA0002566487720000024
N'为滤波后保留的互谱长度，
Figure FDA0002566487720000025
表示第k个频点的相位。
5.根据权利要求2所述一种基于互谱相位拟合的时延估计方法，其特征在于：所述步骤3)包括以下步骤：
31)将声源信号频率范围[f1,f2]除以频率分辨率f0、取整，得到相应的正频率部分频点k值范围[k1,k2]，其中f1表示频率下限，f2表示频率上限；
32)区分正负频率，直接将互谱
Figure FDA0002566487720000026
在频点范围[k1,k2]、[-k2,-k1]以外的值置0的方式达到滤波效果；
33)利用angle函数获得滤波后的互谱相位
Figure FDA0002566487720000027
N'为滤波后保留的互谱长度，
Figure FDA0002566487720000028
表示第k个频点的相位。
6.根据权利要求2所述一种基于互谱相位拟合的时延估计方法，其特征在于：所述步骤4)包括以下步骤：
41)根据互谱相位与频率的线性关系，将
Figure FDA0002566487720000029
除以
Figure FDA00025664877200000210
得到每个频点的时延值
Figure FDA00025664877200000211
p[k]表示第k个频点的时延值；
42)对
Figure FDA00025664877200000212
进行最小二乘线性拟合，确定两路接收信号
Figure FDA00025664877200000213
Figure FDA00025664877200000214
分数部分时延估计值
Figure FDA00025664877200000215
43)与步骤1)得到的整数部分时延m求和得到整体时延值Δn。
CN202010633065.5A 2020-07-02 2020-07-02 一种基于互谱相位拟合的时延估计方法 Withdrawn CN111856401A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010633065.5A CN111856401A (zh) 2020-07-02 2020-07-02 一种基于互谱相位拟合的时延估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010633065.5A CN111856401A (zh) 2020-07-02 2020-07-02 一种基于互谱相位拟合的时延估计方法

Publications (1)

Publication Number Publication Date
CN111856401A true CN111856401A (zh) 2020-10-30

Family

ID=73152157

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010633065.5A Withdrawn CN111856401A (zh) 2020-07-02 2020-07-02 一种基于互谱相位拟合的时延估计方法

Country Status (1)

Country Link
CN (1) CN111856401A (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112512116A (zh) * 2020-11-25 2021-03-16 华东师范大学 一种基于5gsrs信号的自适应二次互相关toa估计方法
CN113589054A (zh) * 2021-07-23 2021-11-02 中国电子科技集团公司第五十四研究所 一种测控天线组阵自动时延偏差估计与补偿装置
CN113777412A (zh) * 2021-08-11 2021-12-10 中电科思仪科技股份有限公司 一种提高天线方向图零深位置测试精度的方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104237871A (zh) * 2013-06-08 2014-12-24 中国科学院声学研究所 一种基于相位补偿的时延差估计方法
CN107644650A (zh) * 2017-09-29 2018-01-30 山东大学 一种基于渐进串行正交化盲源分离算法的改进声源定位方法及其实现系统
CN109669159A (zh) * 2019-02-21 2019-04-23 深圳市友杰智新科技有限公司 基于麦克风十字环阵列的声源定位跟踪装置及方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104237871A (zh) * 2013-06-08 2014-12-24 中国科学院声学研究所 一种基于相位补偿的时延差估计方法
CN107644650A (zh) * 2017-09-29 2018-01-30 山东大学 一种基于渐进串行正交化盲源分离算法的改进声源定位方法及其实现系统
CN109669159A (zh) * 2019-02-21 2019-04-23 深圳市友杰智新科技有限公司 基于麦克风十字环阵列的声源定位跟踪装置及方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
刘超: "基于麦克风阵列的声源定位算法研究", 《中国优秀硕士学位论文全文数据库_信息科技辑》 *
杨祥清 等: "基于麦克风阵列的三维声源定位算法及其实现", 《声学技术》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112512116A (zh) * 2020-11-25 2021-03-16 华东师范大学 一种基于5gsrs信号的自适应二次互相关toa估计方法
CN112512116B (zh) * 2020-11-25 2022-12-02 华东师范大学 一种基于5gsrs信号的自适应二次互相关toa估计方法
CN113589054A (zh) * 2021-07-23 2021-11-02 中国电子科技集团公司第五十四研究所 一种测控天线组阵自动时延偏差估计与补偿装置
CN113777412A (zh) * 2021-08-11 2021-12-10 中电科思仪科技股份有限公司 一种提高天线方向图零深位置测试精度的方法
CN113777412B (zh) * 2021-08-11 2024-03-19 中电科思仪科技股份有限公司 一种提高天线方向图零深位置测试精度的方法

Similar Documents

Publication Publication Date Title
CN111856401A (zh) 一种基于互谱相位拟合的时延估计方法
CN107644650B (zh) 一种基于渐进串行正交化盲源分离算法的改进声源定位方法及其实现系统
CN110837001B (zh) 一种电力系统中谐波和间谐波的分析方法与装置
Ferguson Time-delay estimation techniques applied to the acoustic detection of jet aircraft transits
CN111024209A (zh) 一种适用于矢量水听器的线谱检测方法
Josso et al. Source motion detection, estimation, and compensation for underwater acoustics inversion by wideband ambiguity lag-Doppler filtering
CN112926014A (zh) 基于rls与rssd的滚动轴承声信号多频带融合故障诊断方法
CN110196407B (zh) 一种基于频率预估的单矢量水听器信号来波方向估计方法
JP2007240605A (ja) 複素ウェーブレット変換を用いた音源分離方法、および音源分離システム
Niu et al. Mode separation with one hydrophone in shallow water: A sparse Bayesian learning approach based on phase speed
EP1480589A1 (en) Method and apparatus for removing noise from electronic signals
CN109901114B (zh) 一种适用于声源定位的时延估计方法
US20030128848A1 (en) Method and apparatus for removing noise from electronic signals
CN115015942B (zh) 一种自适应水下目标声学测速装置及方法
CN114813129B (zh) 基于wpe与emd的滚动轴承声信号故障诊断方法
CN112666520B (zh) 一种可调响应时频谱声源定位方法及系统
CN110426711B (zh) 一种基于极性零点检测的时延估计方法及系统
CN106356069B (zh) 一种信号处理方法和装置
CN113820004A (zh) 一种鲁棒的振动信号初始相位估计方法
CN109960843B (zh) 一种基于正交原理的多普勒频移数值仿真方法
CN116996137B (zh) 一种基于加权叠加的低信噪比宽带线性调频信号检测方法
Angelopoulos et al. Nonparametric spectral estimation-an overview
CN114822579B (zh) 一种基于一阶差分麦克风阵列的信号估计方法
CN112505627B (zh) 基于信号过零点信息的时延估计算法
Nakayama et al. Acoustic distance measurement based on phase interference using the cross-spectral method with adjacent microphones

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
WW01 Invention patent application withdrawn after publication

Application publication date: 20201030

WW01 Invention patent application withdrawn after publication