CN106531145B - 基于m估计器的递归有源噪声控制方法 - Google Patents

基于m估计器的递归有源噪声控制方法 Download PDF

Info

Publication number
CN106531145B
CN106531145B CN201611079136.1A CN201611079136A CN106531145B CN 106531145 B CN106531145 B CN 106531145B CN 201611079136 A CN201611079136 A CN 201611079136A CN 106531145 B CN106531145 B CN 106531145B
Authority
CN
China
Prior art keywords
current time
noise
value
filter
estimator
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
CN201611079136.1A
Other languages
English (en)
Other versions
CN106531145A (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.)
Southwest Jiaotong University
Original Assignee
Southwest Jiaotong 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 Southwest Jiaotong University filed Critical Southwest Jiaotong University
Priority to CN201611079136.1A priority Critical patent/CN106531145B/zh
Publication of CN106531145A publication Critical patent/CN106531145A/zh
Application granted granted Critical
Publication of CN106531145B publication Critical patent/CN106531145B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10K—SOUND-PRODUCING DEVICES; METHODS OR DEVICES FOR PROTECTING AGAINST, OR FOR DAMPING, NOISE OR OTHER ACOUSTIC WAVES IN GENERAL; ACOUSTICS NOT OTHERWISE PROVIDED FOR
    • G10K11/00—Methods or devices for transmitting, conducting or directing sound in general; Methods or devices for protecting against, or for damping, noise or other acoustic waves in general
    • G10K11/16—Methods or devices for protecting against, or for damping, noise or other acoustic waves in general
    • G10K11/175—Methods or devices for protecting against, or for damping, noise or other acoustic waves in general using interference effects; Masking sound
    • G10K11/178—Methods or devices for protecting against, or for damping, noise or other acoustic waves in general using interference effects; Masking sound by electro-acoustically regenerating the original acoustic waves in anti-phase
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10K—SOUND-PRODUCING DEVICES; METHODS OR DEVICES FOR PROTECTING AGAINST, OR FOR DAMPING, NOISE OR OTHER ACOUSTIC WAVES IN GENERAL; ACOUSTICS NOT OTHERWISE PROVIDED FOR
    • G10K2210/00—Details of active noise control [ANC] covered by G10K11/178 but not provided for in any of its subgroups
    • G10K2210/30—Means
    • G10K2210/301—Computational
    • G10K2210/3023—Estimation of noise, e.g. on error signals
    • G10K2210/30232—Transfer functions, e.g. impulse response

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Acoustics & Sound (AREA)
  • Multimedia (AREA)
  • Soundproofing, Sound Blocking, And Sound Damping (AREA)

Abstract

基于M估计器的递归有源噪声控制方法，其步骤主要是：A、噪声输入向量的生成；B、权向量的生成；C、滤波器输出值的生成；D、残差信号阈值的确定，D1、置于噪声消除点的误差麦克风拾取出噪声消除点在当前时刻n的声音信号；并得到残差信号平方e2(n)，组成残差信号平方窗函数向量E(n)，E(n)＝[e2(n),e2(n‑1),...,e2(n‑Nw+1)]；D2、计算M估计器的期望方差D3、计算出M估计器残差信号的低限阈值ξ、中限阈值Δ1和高限阈值Δ2；E、误差信号限定值的计算，滤波器计算出M估计器的误差信号限定值ψ(n)；F、权向量更新；G、迭代，令n＝n+1,重复A、B、C、D、E、F的步骤,直至滤波结束。该方法对冲击噪声的消除效果好，收敛速度快，残留误差小。

Description

基于M估计器的递归有源噪声控制方法
技术领域
本发明属于自适应有源噪声控制技术领域。
背景技术
噪音控制一直以来是一个重要的课题。随着科技的发展，噪声污染的问题越来越严重。很多电子设备及工业装置都会产生噪声。这些噪声不仅可能对通信或自动控制系统造成干扰，还会影响人的身心健康。噪声控制技术主要分为无源控制技术和有源(主动)控制技术。无源控制技术，主要采用消声器、振动隔离、阻尼减振等方法，不仅造价高，且对于低频噪声往往没有效果。有源控制技术，采用主动噪声控制(Active noise control,ANC)系统，利用两列频率相同、相位差固定的声波，叠加后会产生相加性或相消性干涉，从而使声能得到增强或减弱，利用声波相消性干涉来消除噪声。此类方法，有着有效的去噪能力(尤其是低频噪声)，且成本低；近年来得到了广泛的关注。
在一个有源(主动)噪声控制系统中，主麦克风置于噪声源处用于拾取噪声源的原始噪声(初级噪声)，控制器通常是一个自适应滤波器，自适应滤波器对噪声信号滤波后通过扬声器输出一个反向(相)噪声(次级噪声)；误差麦克风置于噪声消除点(扬声器处)，用于拾取噪声消除点的残差信号(初级噪声与次级噪声抵消后的误差信号)，该残差信号用于主动噪声控制器的反馈修正。从噪声源到噪声消除点的路径称为主级通路，噪声源经控制器到噪声消除点(扬声器)的路径称为次级通路。
自适应有源噪声控制方法中，滤波-x最小均方(FxLMS)算法应用最为广泛，对于初级噪声信号为非冲击噪声的环境下，有着良好的鲁棒性，且结构简单。但是，当初级噪声为α稳定分布的冲击噪声时，FxLMS算法会出现不收敛的情况。由于FxLMS算法采用均方误差，而α稳定分布的噪声方差是无穷的，这样会导致控制器更新时参数发散。为避免这一缺陷，一些改进方法被相继提出，较成熟的方法有以下两种：
(1)滤波p范数(FxLMP)方法
参考文献1“Leahy R,Zhou Z,Hsu YC.Adaptive filtering of stableprocesses for active attenuation of impulsive noise.in:Proceedings of the1995International Conference on Acoustics,Speech,and Signal Processing,vol.5；1995.p.2983-2986.)该方法是将分数阶误差策略，即p范数(1<p<2)引入有源噪声控制系统中，使得当强噪声环境下，能够把滤波器限定在一个范围内，避免其参数发散，从而提高了算法在冲击噪声下的稳定性，也提高了算法收敛速度。但是实验结果表面该方法某些强冲击收敛速度仍然不够快。
(2)鲁棒FxLMS(RFxLMS)方法
参考文献2“George NV,Panda G.A robust filtered-s LMS algorithm fornonlinear active noise control.Appl.Acoust.2012；73(8):836–41.”。该文献提出了一种鲁棒的FxLMS消除方法。该方法采用一种基于对数的误差信号更新，而这种对数变换可以在一定程度上压缩α稳定分布噪声的幅值，从而限定控制器的参数，提高了算法在大部分情况下的稳定性，但是该算法具有较大的残差。
发明内容
本发明的发明目的就是提供一种基于M估计器的递归有源噪声控制方法，该方法对冲击噪声的消除效果好，收敛速度快，残留误差小。
本发明实现其发明目的所采用的技术方案是，基于M估计器的递归有源噪声控制方法，其步骤如下：
1、基于M估计器的递归有源噪声控制方法，其步骤如下：
A、噪声输入向量的生成
置于噪声源的主麦克风拾取出噪声源当前时刻n产生的初级噪声值u(n)；将当前时刻n和前L-1个时刻的L个初级噪声值u(n),u(n-1)...,u(n-L+1)，构成滤波器当前时刻n的噪声输入向量U(n)，U(n)＝[u(n),u(n-1)...,u(n-L+1)]T；其中L＝128，是滤波器的抽头数，上标T代表转置；
B、权向量的生成
滤波器生成当前时刻n和前L-1个时刻的权系数w(n),w(n-1)…w(n-L+1)，并将这L个权系数构成当前时刻n的权向量W(n)，W(n)＝[w(n),w(n-1),...,w(n-L+1)]；当前时刻n<129时，W(n)＝0；
C、滤波器输出值的生成
滤波器根据当前时刻n的权向量W(n)对当前时刻n的噪声输入向量U(n)进行滤波得到当前时刻n输出到扬声器的滤波器输出值 其中符号*代表卷积运算，s表示噪声源经滤波器到扬声器的次级通路固有的脉冲响应；
D、残差信号阈值的确定
D1、置于噪声消除点的误差麦克风拾取出噪声消除点在当前时刻n的声音信号，即当前时刻n的残差信号e(n)；并得到当前时刻n的残差信号平方e2(n)，将当前时刻n和前Nw-1个时刻的残差信号平方e2(n),e2(n-1),...,e2(n-Nw+1)，组成残差信号平方窗函数向量E(n)，E(n)＝[e2(n),e2(n-1),...,e2(n-Nw+1)]，其中，Nw是窗函数长度，取值为100；
D2、计算M估计器的期望方差
其中，C为计算纠正因子，C＝1.483{1+5/(Nw-1)}；ζ为缩放参数，其取值为0.9～1；med为取中位数运算；
D3、计算出M估计器残差信号的低限阈值ξ、中限阈值Δ1和高限阈值Δ2，
E、误差信号限定值的计算
滤波器根据当前时刻n的残差信号e(n)，计算出M估计器当前时刻n的误差信号限定值ψ(n)：
其中，sign表示符号函数，当e(n)>0时，其值为1；当e(n)<0时，其值为-1；当e(n)＝0时，其值为0；
F、权向量更新
滤波器生成当前时刻n的增益向量Φ(n)，
式中，λ表示遗忘因子，其取值为0.9～1；F(n)为当前时刻n的递归矩阵，其大小为L×L，当前时刻n<129时，F(n)＝0；下一时刻(n+1)的递归矩阵F(n+1)由公式F(n+1)＝λ-1F(n)-λ-1Φ(n)(s*U(n))TF(n)得到；
随后，滤波器算出下一时刻(n+1)的权向量W(n+1)，W(n+1)＝W(n)+Φ(n)e(n)
G、迭代
令n取n+1的值,重复A、B、C、D、E、F的步骤,直至滤波结束。
与现有技术相比，本发明的有益效果是：
本发明通过对残差信号进行递归M估计器运算，M估计器是对误差信号的鲁棒估计，当模型假设发生某种微小变化时，相应的估计器统计推断也只有微小改变。当误差信号限定值是有界函数时，基于M估计器的算法就会更加稳定，其原因是即使出现异常值，也不会对算法中滤波器权向量的更新产生无法控制的影响：当残差信号e(n)位于误差信号的低位阈值(下限)附近时，外界的信号的变化只会引起滤波器的参数微小变化；当残差信号e(n)位于误差信号中位阈值附近时，滤波器参数的变化与外界环境的变化成正相关；当残差信号e(n)位于误差信号的高限阈值附近时，滤波器参数的变化不再受外界环境的影响。因此，本发明对噪声信号有着相较于传统算法有着更好的信息采集、追踪能力；尤其对含有各种高阶信号的冲击噪声有着非常稳定的去噪性能，且收敛速度快。
仿真实验证明，本发明的收敛速度快，较现有方法其收敛速度提高了大约50％；去噪性能好，残差功率小，其平均噪声残留较现有方法，降低了约5dB。
下面结合附图和具体实施方式对本发明进行详细说明
附图说明
图1是本发明方法和文献1、文献2方法的仿真实验的平均噪声残留对比图。
图2是本发明方法和文献1、文献2方法的仿真实验的残留噪声频谱密度对比图。
下面结合附图和具体实施方式对本发明进行详细说明。
具体实施方式
实施例
本发明的一种具体实施方式是，基于M估计器的递归有源噪声控制方法，其步骤如下：
1、基于M估计器的递归有源噪声控制方法，其步骤如下：
A、噪声输入向量的生成
置于噪声源的主麦克风拾取出噪声源当前时刻n产生的初级噪声值u(n)；将当前时刻n和前L-1个时刻的L个初级噪声值u(n),u(n-1)...,u(n-L+1)，构成滤波器当前时刻n的噪声输入向量U(n)，U(n)＝[u(n),u(n-1)...,u(n-L+1)]T；其中L＝128，是滤波器的抽头数，上标T代表转置；
B、权向量的生成
滤波器生成当前时刻n和前L-1个时刻的权系数w(n),w(n-1)…w(n-L+1)，并将这L个权系数构成当前时刻n的权向量W(n)，W(n)＝[w(n),w(n-1),...,w(n-L+1)]；当前时刻n<129时，W(n)＝0；
C、滤波器输出值的生成
滤波器根据当前时刻n的权向量W(n)对当前时刻n的噪声输入向量U(n)进行滤波得到当前时刻n输出到扬声器的滤波器输出值 其中符号*代表卷积运算，s表示噪声源经滤波器到扬声器的次级通路固有的脉冲响应；
D、残差信号阈值的确定
D1、置于噪声消除点的误差麦克风拾取出噪声消除点在当前时刻n的声音信号，即当前时刻n的残差信号e(n)；并得到当前时刻n的残差信号平方e2(n)，将当前时刻n和前Nw-1个时刻的残差信号平方e2(n),e2(n-1),...,e2(n-Nw+1)，组成残差信号平方窗函数向量E(n)，E(n)＝[e2(n),e2(n-1),...,e2(n-Nw+1)]，其中，Nw是窗函数长度，取值为100；
D2、计算M估计器的期望方差
其中，C为计算纠正因子，C＝1.483{1+5/(Nw-1)}；ζ为缩放参数，其取值为0.9～1；med为取中位数运算；
D3、计算出M估计器残差信号的低限阈值ξ、中限阈值Δ1和高限阈值Δ2，
E、误差信号限定值的计算
滤波器根据当前时刻n的残差信号e(n)，计算出M估计器当前时刻n的误差信号限定值ψ(n)：
其中，sign表示符号函数，当e(n)>0时，其值为1；当e(n)<0时，其值为-1；当e(n)＝0时，其值为0；
F、权向量更新
滤波器生成当前时刻n的增益向量Φ(n)，
式中，λ表示遗忘因子，其取值为0.9～1；F(n)为当前时刻n的递归矩阵，其大小为L×L，当前时刻n<129时，F(n)＝0；下一时刻(n+1)的递归矩阵F(n+1)由公式F(n+1)＝λ-1F(n)-λ-1Φ(n)(s*U(n))TF(n)得到；
随后，滤波器算出下一时刻(n+1)的权向量W(n+1)，W(n+1)＝W(n)+Φ(n)e(n)
G、迭代
令n取n+1的值,重复A、B、C、D、E、F的步骤,直至滤波结束。
仿真实验：
为了验证本发明的有效性，进行了仿真实验，并与文献1-2的算法进行了对比。
仿真实验的主级通路和次级通路均采用高阶FIR滤波器建模。有源噪声控制器的阶数设定为128阶。初级噪声分别为一个α稳定分布的冲击噪声，α取值越大噪声越弱。ζ＝0.999。
图1是本发明方法和文献1、文献2方法的仿真实验的平均噪声残留(滤波后的残差信号与未经滤波的初级噪声传输到消除点的信号之比)对比图。由图1可知，本发明的收敛速度快，在大约10000时刻收敛，而文献1滤波器参数发散；文献2在大约在15000时刻收敛；本发明的收敛速度较效果好的文献2的方法，提高了大约50％。且平均噪声残留较文献1，文献2的方法，降低了约5dB。
图2是本发明方法和文献1、文献2方法的仿真实验的残留噪声频谱密度对比图。由图2可知，本发明比文献1、2收敛速度更快，残差功率更小。

Claims (1)

1.基于M估计器的递归有源噪声控制方法，其步骤如下：
A、噪声输入向量的生成
置于噪声源的主麦克风拾取出噪声源当前时刻n产生的初级噪声值u(n)；将当前时刻n和前L-1个时刻的L个初级噪声值u(n),u(n-1)...,u(n-L+1)，构成滤波器当前时刻n的噪声输入向量U(n)，U(n)＝[u(n),u(n-1)...,u(n-L+1)]T；其中L＝128，是滤波器的抽头数，上标T代表转置；
B、权向量的生成
滤波器生成当前时刻n和前L-1个时刻的权系数w(n),w(n-1)…w(n-L+1)，并将这L个权系数构成当前时刻n的权向量W(n)，W(n)＝[w(n),w(n-1),...,w(n-L+1)]；当前时刻n<129时，W(n)＝0；
C、滤波器输出值的生成
滤波器根据当前时刻n的权向量W(n)对当前时刻n的噪声输入向量U(n)进行滤波得到当前时刻n输出到扬声器的滤波器输出值 其中符号*代表卷积运算，s表示噪声源经滤波器到扬声器的次级通路固有的脉冲响应；
D、残差信号阈值的确定
D1、置于噪声消除点的误差麦克风拾取出噪声消除点在当前时刻n的声音信号，即当前时刻n的残差信号e(n)；并得到当前时刻n的残差信号平方e2(n)，将当前时刻n和前Nw-1个时刻的残差信号平方e2(n),e2(n-1),...,e2(n-Nw+1)，组成残差信号平方窗函数向量E(n)，E(n)＝[e2(n),e2(n-1),...,e2(n-Nw+1)]，其中，Nw是窗函数长度，取值为100；
D2、计算M估计器的期望方差
其中，C为计算纠正因子，C＝1.483{1+5/(Nw-1)}；ζ为缩放参数，其取值为0.9～1；med为取中位数运算；
D3、计算出M估计器残差信号的低限阈值ξ、中限阈值Δ1和高限阈值Δ2，
E、误差信号限定值的计算
滤波器根据当前时刻n的残差信号e(n)，计算出M估计器当前时刻n的误差信号限定值ψ(n)：
其中，sign表示符号函数，当e(n)>0时，其值为1；当e(n)<0时，其值为-1；当e(n)＝0时，其值为0；
F、权向量更新
滤波器生成当前时刻n的增益向量Φ(n)，
式中，λ表示遗忘因子，其取值为0.9～1；F(n)为当前时刻n的递归矩阵，其大小为L×L，当前时刻n<129时，F(n)＝0；下一时刻(n+1)的递归矩阵F(n+1)由公式F(n+1)＝λ-1F(n)-λ-1Φ(n)(s*U(n))TF(n)得到；
随后，滤波器算出下一时刻(n+1)的权向量W(n+1)，W(n+1)＝W(n)+Φ(n)e(n)
G、迭代
令n取n+1的值,重复A、B、C、D、E、F的步骤,直至滤波结束。
CN201611079136.1A 2016-11-30 2016-11-30 基于m估计器的递归有源噪声控制方法 Active CN106531145B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201611079136.1A CN106531145B (zh) 2016-11-30 2016-11-30 基于m估计器的递归有源噪声控制方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201611079136.1A CN106531145B (zh) 2016-11-30 2016-11-30 基于m估计器的递归有源噪声控制方法

Publications (2)

Publication Number Publication Date
CN106531145A CN106531145A (zh) 2017-03-22
CN106531145B true CN106531145B (zh) 2019-05-17

Family

ID=58355238

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201611079136.1A Active CN106531145B (zh) 2016-11-30 2016-11-30 基于m估计器的递归有源噪声控制方法

Country Status (1)

Country Link
CN (1) CN106531145B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109147753B (zh) * 2018-07-24 2022-03-25 西南交通大学 基于误差平方与误差平方对数之差最小的凸组合降噪方法
CN111785243B (zh) * 2020-06-10 2024-02-09 厦门市派美特科技有限公司 一种耳机有源降噪方法及系统
CN112087178A (zh) * 2020-09-14 2020-12-15 四川长虹电器股份有限公司 一种永磁同步电机的参数识别方法
CN113295753A (zh) * 2021-04-29 2021-08-24 山东非金属材料研究所 一种pH计电压误差的双FIR测量方法
CN114488224B (zh) * 2021-12-24 2023-04-07 西南交通大学 一种用于卫星集中式自主导航的自适应滤波方法
CN117451951B (zh) * 2023-12-22 2024-03-19 微瑞精仪(厦门)科技有限公司 一种高适应性单分子测量仪

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1886104A (zh) * 2003-11-26 2006-12-27 加利福尼亚大学董事会 有源噪声控制方法以及包含前馈和反馈控制器的设备
CN104683614A (zh) * 2015-03-24 2015-06-03 西南交通大学 基于m估计的记忆成比例仿射投影自适应回声消除方法
CN105976806A (zh) * 2016-04-26 2016-09-28 西南交通大学 基于最大熵的有源噪声控制方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8824666B2 (en) * 2009-03-09 2014-09-02 Empire Technology Development Llc Noise cancellation for phone conversation

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1886104A (zh) * 2003-11-26 2006-12-27 加利福尼亚大学董事会 有源噪声控制方法以及包含前馈和反馈控制器的设备
CN104683614A (zh) * 2015-03-24 2015-06-03 西南交通大学 基于m估计的记忆成比例仿射投影自适应回声消除方法
CN105976806A (zh) * 2016-04-26 2016-09-28 西南交通大学 基于最大熵的有源噪声控制方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于M估计的混合噪声主动控制算法;田亚雄;《佳木斯大学学报（自然科学版）》;20160915;第34卷(第5期);791-794

Also Published As

Publication number Publication date
CN106531145A (zh) 2017-03-22

Similar Documents

Publication Publication Date Title
CN106531145B (zh) 基于m估计器的递归有源噪声控制方法
Lu et al. A survey on active noise control in the past decade—Part I: Linear systems
Thenua et al. Simulation and performance analysis of adaptive filter in noise cancellation
CN105976806B (zh) 基于最大熵的有源噪声控制方法
CN107333017A (zh) 一种重用权系数指数型误差自适应回声消除方法
Xiao et al. A new efficient filtered-x affine projection sign algorithm for active control of impulsive noise
CN104035332B (zh) 基于m‑估计的冲击类噪声有源控制方法
CN109119061A (zh) 一种基于梯度矩阵的有源噪声控制方法
CN113078884B (zh) 添加非线性拟合的自适应算法
Gomathi et al. Variable step size for improving convergence of FxLMS algorithm
Kar et al. An improved filtered-x least mean square algorithm for acoustic noise suppression
CN113299260A (zh) 一种基于emfnl滤波器在线建模次级通道的有源降噪方法
CN111193497A (zh) 一种基于emfnl滤波器的次级通道建模方法
Hassanpour et al. An efficient online secondary path estimation for feedback active noise control systems
Avalos et al. Multichannel Filtered-X Error Coded Affine Projection-Like Algorithm with Evolving Order
CN113851104A (zh) 含次级通道在线辨识的反馈型主动噪声控制系统及方法
Mehmood et al. A comprehensive review of adaptive noise cancellation techniques in the internet of things
Ramadan et al. An adaptive noise canceller using error nonlinearities in the lms adaptation
Akhtar et al. Robust adaptive algorithm for active noise control of impulsive noise
Zhou et al. Robust generalized hyperbolic secant algorithm for nonlinear active noise control
Lu et al. A recursive least M-estimate adaptive controller for active impulsive noise control
Mohapatra et al. A sigmoid function based feedback filtered-X-LMS algorithm with improved offline modelling
Dhubkarya et al. Comparative Performance Analysis of Adaptive Algorithms for Simulation & Hardware Implementation of an ECG Signal
Foo et al. Active noise cancellation headset
Rodríguez et al. Convex Combination of FXECAP–FXECLMS Algorithms for Active Noise Control

Legal Events

Date Code Title Description
C06 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