CN105072557B - 一种三维环绕声重放系统的扬声器环境自适应校准方法 - Google Patents

一种三维环绕声重放系统的扬声器环境自适应校准方法 Download PDF

Info

Publication number
CN105072557B
CN105072557B CN201510490058.3A CN201510490058A CN105072557B CN 105072557 B CN105072557 B CN 105072557B CN 201510490058 A CN201510490058 A CN 201510490058A CN 105072557 B CN105072557 B CN 105072557B
Authority
CN
China
Prior art keywords
transmission function
function
signal
sound
output signal
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
CN201510490058.3A
Other languages
English (en)
Other versions
CN105072557A (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.)
Peking University
Original Assignee
Peking 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 Peking University filed Critical Peking University
Priority to CN201510490058.3A priority Critical patent/CN105072557B/zh
Publication of CN105072557A publication Critical patent/CN105072557A/zh
Application granted granted Critical
Publication of CN105072557B publication Critical patent/CN105072557B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Stereophonic System (AREA)

Abstract

本发明公开了一种三维环绕声重放系统的扬声器环境自适应校准方法。本方法为:1)测量三维环绕声重放系统的听音环境内每一通道扬声器到听音位置的传递函数hi(n);2)计算每一传递函数hi(n)的逆函数gi(n);3)测量激励信号源的虚拟方位(θ,δ)以及各通道扬声器的摆放位置(θi,δi),依据三维环绕声算法计算各通道扬声器的理想输出信号ti;4)将各通道扬声器的理想输出信号ti与相应通道的逆函数gi(n)进行卷积计算,得到各通道扬声器的实际输出信号t'i,其中,i=1…M,M为扬声器总数,n代表时间。此方法可均衡由扬声器性能、非球面均匀分布、传输信道的不一致性。

Description

一种三维环绕声重放系统的扬声器环境自适应校准方法
技术领域
本发明属于三维环绕声重放技术领域,本发明提出了一种扬声器环境自适应校准方法,解决了三维环绕声重放系统对扬声器摆放位置固定,幅频特性一致的要求。
背景技术
3D多媒体的时代已经到来,3D音视频系统也正迅速走向电影院,家庭影院,及手持终端设备,成为全球各大电子制造商的新焦点。环绕声重放技术在3D音视频系统占有重要的地位,目前主流的技术包括VBAP(Vector Base Amplitude Panning)、Ambisonics,WFS(Wave Field Synthetize)。其中Ambisonics方法是1973年由牛津大学的Michael Gerzon提出的(参考:Gerzon M.“Periphony:With-Height Sound Reproduction,”Journal oftheAudio Engineering Society,vol.21(1),pp.2-10,1973),主要是通过基于球谐函数对原始声场的分解与重建来控制虚拟声源的方位。以Ambisonics声重放系统为例,基于Ambisonics声重放系统技术特点是编解码分离,在编码阶段,根据虚拟声源的方向得到各球谐基函数的投影值;在声场重放阶段,根据重发扬声器的数量、方位和编码环节得到的投影值,得到不同通道信号的输出增益,把此增益输送给对应的扬声器重发,达到在扬声器阵列中心位置处重建源声场的目的。
尽管Ambisonics方法编解码分离的方案给3D声音录音和重放带来了很大的优势,但在走向市场的道路中却遇到困难,原因之一就是Ambisonics是基于声场重构方法通过复杂的数学计算得到的,它假设复现声场为自由声场,播放设备各通道幅频特性一致,扬声器大致均匀分布在一个以听者为中心的球面上。而这些条件在影院尤其是家庭影院等实际应用中很难满足,导致基于Ambisonics重建的声场出现较大重构误差,无法满足听音需求。
发明内容
针对现有技术中存在的技术问题,本发明的目的在于提供一种应用于三维环绕声重放系统的环境自适应扬声器校准方法,可均衡非自由声场,非一致性信道,以及扬声器非球面分布带来的影响。
本发明的具体思想是,首先测出听音环境中各通道扬声器到听音位置的冲击响应函数(冲击响应函数记录了播放系统及扬声器的幅频特性,声场中扬声器位置到预先设定的听音位置的混响特性),其次求出冲激响应函数的逆系统,最后将计算得到的理想三维环绕声重放系统中各扬声器应播放的音频信号卷积上此逆系统的冲激响应,可均衡非自由声场,非一致性信道,以及扬声器非球面分布带来的影响,解决三维环绕声重放系统中虚拟声源定位不准的问题。
图1为本发明流程框图,分为四个步骤,
1、测量各通道扬声器到听音位置的传递函数h1…M(n),M为扬声器总数。
2、计算传递函数hi(n)的逆函数gi(n),使i=1…M。式中δ(n)为单位冲击函数,n代表时间,n=0时,该函数值为1,其他均为0。
3、已知信号s,测量信号源的虚拟方位(θ,δ)以及各扬声器的摆放位置(θi,δi),依据三维环绕声算法(比如Ambisonics算法),计算各扬声器理想输出信号ti,i=1…M。
4、把各通道扬声器理想的输出信号ti与相应通道的逆函数gi(n)进行卷积得到各扬声器的实际输出信号t'i,i=1…M。
与现有技术相比,本发明的积极效果为:
本发明在听者的位置放置一个麦克风,用MLS序列测量各扬声器的传递函数,再求其逆函数,最后,把逆函数作用于基于三维环绕声技术求出的各音箱播放信号,此方法可均衡由扬声器性能、非球面均匀分布、传输信道的不一致性,以及室内混响带来的影响,通过对理想系统、未均衡系统以及均衡后系统的球谐函数分解系数进行比较可得,均衡后系统优于未均衡系统。
附图说明
图1为本发明方法流程图;
图2为传递函数测量系统图;
图3为基于MLS序列测量传递函数的流程图;
图4为扬声器摆放位置图;
图5为二维空间三阶球谐函数分解系数图;其中,
(a)为理想系统系数图,(b)为未均衡系统系数图,(c)为均衡后系统系数图,
(d)为理想系统系数图,(e)为未均衡系统系数图,(f)为均衡后系统系数图,
(g)为理想系统系数图,(h)为未均衡系统系数图,(i)为均衡后系统系数图,
(j)为理想系统系数图,(k)为未均衡系统系数图,(l)为均衡后系统系数图;
图6为三种系统的空间系数相关图;其中,
(a)为理想系统的相关系数,(b)未均衡系统的相关系数,(c)为均衡后系统的相关系数;
图7为空间指向图;其中,
(a)为理想系统水平角10度,(b)为未均衡系统水平角10度,(c)为均衡后系统水平角10度,(d)为理想系统水平角20度,(e)为未均衡系统水平角20度,(f)均衡化系统水平角20度。
具体实施方式
下面结合附图对本发明进行进一步详细描述,本发明的流程如图1所示。
步骤1:各通道扬声器到听音位置传递函数测量:
传递函数测量方法的激励信号包括最大长度伪随机序列(MLS:Maximum LengthSequence),Gelay码及扫频信号等,本发明中采用最大长度伪随机序列(MLS)作为激励信号,测量听音环境内各扬声器到听音位置的传递函数。具体方法为:
声源到达听音者的位置在音量不是很大时可近似满足线性时不变的条件,我们可以将这个过程看成是一个线性时不变系统,任何系统都有自己的传输特性,即传递函数。我们把输入激励信号用x(n)表示,测试系统的传递函数用hi(n)表示,输出信号用yi(n)表示,则常用的测量问题总是通过测量x(n)和yi(n),来求解二者之间的关系hi(n),如图2所示。其中,yi(n)为hi(n)与x(n)的卷积。在测量传递函数时,测试信号x(n)作为输入信号进入系统,通过测量最后得到的输出信号yi(n)。根据二者关系求得整个系统的传递函数,即为我们想得到的传递函数。
理想的测量传递函数的激励信号应满足的特点:激励信号的序列可重复再生;激励信号是确定信号,激励信号是宽带信号,具有最大的信噪比SNR,具有最小的非线性时变误差。可选用MLS序列,Golay码和扫频信号。本发明采用MLS序列,但不限于MLS序列。
首先,假设测量系统是线性时不变系统。长度为L的MLS序列x(n)的圆自相关函数为:
以上性质说明,x(n)的自相关函数若除以序列长度,则近似于一个脉冲信号。利用MLS序列的以上性质,我们可以推导出MLS序列测量传递函数的原理:
其中,h′i(n)为实际测量得到的传递函数,hi(n)为待测系统的传递函数,x(n)为输入系统的MLS序列,y(n)为系统测量得到的输出信号。
系统得到的输出信号y(n)与输入的MLS序列x(n)计算互相关后,除以序列长度L。当序列长度L足够时,近似于一个脉冲信号。这样的方法得到的h′i(n)近似为hi(n)。在测量时需要注意,MLS序列的长度L需要大于室内环境的混响时间,否则x(n)与互相关中的圆卷积和线性卷积并不相等,会出现时间混叠的问题。测试流程如图3。
所以利用MLS序列测量传递函数的具体过程如下:
1)生成N阶长度为L=2N-1的MLS序列x(n);
2)测量得到系统的输出信号y(n);
3)计算x(n)与y(n)的互相关,并除以序列长度L得到h′i(n);
4)多次测量取平均值。
步骤2:传递函数的逆系统计算
步骤1测量的传递函数包含了播放系统及扬声器的幅频特性,声场中扬声器位置到听音位置的混响特性,本发明第二步采用最小二乘方法,但不限于最小二乘方法,估计步骤1中得到传递函数的逆系统的传递函数,来均衡非自由声场,非一致性信道,以及扬声器非球面分布带来的影响。具体方法描述为,
测得的传递函数hi(n),则逆函数gi(n)应满足
式中T为系统延迟,如果直接对此式求解,解得的gi(n)会出现不稳定的情况,因此本发明采用最小二乘算法,在最小均方意义上解得gi(n)的稳定近似解。
引入w(n)为随机白噪声序列,得到序列f(n)满足下式,
式中hi(n)为传递函数,gi(n)为传递函数hi(n)的逆函数。
理想情况下f(n)=w(n-T),将上式展开:
式中d为逆滤波器的长度,n的取值限定值为e,应大于等于d,保证方程超定。
转化为矩阵形式,
F=BiGi (7)
其中:
其中d为逆系统阶数,e为白噪声长度。
目标即是解出逆系统的冲激响应G,在此式中,问题已转化为解超定方程了,可以采用最小二乘法求解。
步骤3:基于Ambisonics的三维环绕声算法
根据平面声波可以在球坐标系中用球谐函数无限展开的原理,振幅为s的平面波传播到内部没有声源的有限空间Ω的波动方程也可以写为:
用球谐函数分解声场
结合得
简化为
由于具有特定的与声源方位相关的空间指向特性,因此可认为,平面波S的空间信息可被空间采样并保留下来,这个过程也被称为Ambisonic系统的编码过程。也被称为原始Ambisonic信号。当然,在实施过程中不可能进行无限阶球谐函数取样,只能到M阶截断。此时3维情况下,对平面波S进行M阶采样后共有(M+1)2个原始Ambisonic信号。而2维情况下,对平面波S进行M阶采样后共2M+1个原始Ambisonic信号。
根据球谐函数重构声场的理论,用Ambisonic系统对声场进行重构时,根据重发扬声器的数量和方位,将由编码环节得到的原始Ambisonic信号以不同的比例作线性混合,然后输送给对应的扬声器重发,最终达到在扬声器阵列中心位置处重建源声场的目的。确定各原始Ambisonic信号的混合比例的过程就是Ambisonic系统的解码过程。假设方位角为的第j个扬声器的重发信号为tj,则按照Ambisonic系统重构声场的原理,二维情况时tj应满足以下方程式:
简化表达
Y·=A (16)
g可以用以下方法解出
t=pinv(g)·A=(YT·Y)-1·YT·A (17)
当然,扬声器数量越多,重发效果越好,因此二维情况时,重发扬声器的数量推荐采用k≥2M+2.
步骤4:卷积计算输出
步骤3求出了个通道理想的音频输出信号,但此信号经扬声器播放时,必然受到信道及房间环境的影响,为均衡此影响,把此信号卷积上步骤2求出的逆函数g(n),即可得到输出到各扬声器的音频信号,
客观评测实验
本发明采用听音位置复现的声场与理想声场的相似程度来评价所采用的方法,根据声场的球谐函数分解原理得到空间任意一点的声压的分解公式为
虚拟声源最直接的效果即虚拟声源与预设方位相差多少,而虚拟声源的方位则体现在球谐函数的各阶系数上,因此本发明主要用复现的声场与理想声场各自的球谐函数系数的相似度来作为一个客观评价指标。
参与评价的系统有三种,1.理想情况下的声场;2.未加均衡处理情况下的声场;3.采用本发明方法均衡处理后的声场。其各自球谐函数的系数为:
1.理想情况下的声场
2.未加均衡处理情况下的声场,因为扬声器幅频特性不一致,到达中心点距离不一致,传播路径不同,导致扬声器信号到达听者位置时,相当于经过了一个滤波器,此滤波器既是步骤1中测得的传递函数。因而,计算各扬声器的球谐函数分解系数也要乘以滤波器系数
3.采用本发明方法均衡处理后的声场,相当于在上系统2基础上再乘以步骤3求出的逆系统的传递函数,
系统评价环境如图4所示,共用11个扬声器。
采用二维空间三维球谐函数分解,分别按式(14,19,20)计算三个系统的各阶球谐函数分解系数,结果如图5所示,图中显示,均衡后的系统空间参数分布更贴近理想系统,优于理想系统。为了评价球谐系数,本发明采用复现的声场与理想声场各自的球谐函数系数的相关系数来作为一个客观评价指标,其具体定义为:
设三种评价的系统(1.理想系统;2.未加均衡处系统;3.采用本发明方法均衡处理后的系统)各自球谐函数转为向量为,
1.理想系统
r1(f)=A1·A1/(A1|*|A1|) (23)
2.未加均衡系统,
r2(f)=A2(f)·A1/(|A1|*|A2(f)|) (24)
3.采用本发明方法均衡处理后的系统
r3(f)=A3(f)·A1/(|A1|*|A3(f)|) (25)
图6为三种系统的相关系数图,左图为理想系统的相关系数r1,中图未均衡系统的相关系数r2,右图为均衡后系统的相关系数r3。从图中可以看出,均衡后系统的空间系数相关图优于未均衡系统。
空间增益图可以说明了所实现系统空间方位的分辨率和混淆度,本发明采用空间指向图来进一步评价所实现系统的性能。三种系统的空间指向图计算公式如下,
1.理想系统
2.未加均衡处理情况下的声场,
3.采用本发明方法均衡处理后的声场
num这里取3。
空间指向里选则水平角为0°,20°情况,计算结果如图7所示,上图为水平角10度,下图为水平角20度,左图为理想系统,中图为未均衡系统,右图为均衡后系统。从图中可以看出,均衡后的系统指向图优于未均衡系统。
综上,我们采用了三种对空间参数的评价标准,均衡后的系统均优于未经处理的实际系统,本发明方法有效。本发明基于各通道扬声器到听音者的传递函数,基于一定的优化准则法估计其逆函数,利用此逆函数均衡听音环境及信道对各通道声音在幅值和相位上的影响,客观实验结果表明了本发明的有效性。

Claims (7)

1.一种三维环绕声重放系统的扬声器环境自适应校准方法,其步骤为:
1)测量三维环绕声重放系统的听音环境内每一通道扬声器到听音位置的传递函数hi(n);
2)计算每一传递函数hi(n)的逆函数gi(n);
3)测量激励信号源的虚拟方位(θ,δ)以及各通道扬声器的摆放位置(θi,δi),依据三维环绕声算法计算各通道扬声器的理想输出信号ti
4)将各通道扬声器的理想输出信号ti与相应通道的逆函数gi(n)进行卷积计算,得到各通道扬声器的实际输出信号t'i,其中,i=1…M,M为扬声器总数,n代表时间;
其中,采用最小二乘算法计算出所述逆函数gi(n),其方法为:首先引入一随机白噪声序列w(n),得到一输出序列然后根据f(n)=w(n-T)对输出序列进行展开并转换为一矩阵F=BiGi;其中,T为系统延迟,d为逆系统阶数,e为白噪声阶数,然后采用最小二乘法对矩阵F进行求解,得到所述逆函数gi(n)。
2.如权利要求1所述的方法,其特征在于,测量所述传递函数的激励信号为MLS序列。
3.如权利要求2所述的方法,其特征在于,获取所述传递函数hi(n)的方法为:
31)生成N阶长度为L=2N-1的MLS序列x(n);
32)根据输入序列x(n),测量得到系统的输出信号yi(n);
33)计算序列x(n)与信号yi(n)的互相关,并除以序列长度L得到本次测量的传递函数h′i(n);
34)重复步骤31)~33),将多次测量的传递函数h′i(n)取平均值,得到所述传递函数hi(n)。
4.如权利要求2或3所述的方法,其特征在于,所述MLS序列的长度L大于听音环境的混响时间。
5.如权利要求1所述的方法,其特征在于,所述扬声器的数量k≥2M+2,其中,M为三维环绕声重放系统编码过程中球谐函数展开的截断阶数。
6.如权利要求1所述的方法,其特征在于,测量所述传递函数的激励信号为Golay码。
7.如权利要求1所述的方法,其特征在于,测量所述传递函数的激励信号为扫频信号。
CN201510490058.3A 2015-08-11 2015-08-11 一种三维环绕声重放系统的扬声器环境自适应校准方法 Active CN105072557B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510490058.3A CN105072557B (zh) 2015-08-11 2015-08-11 一种三维环绕声重放系统的扬声器环境自适应校准方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510490058.3A CN105072557B (zh) 2015-08-11 2015-08-11 一种三维环绕声重放系统的扬声器环境自适应校准方法

Publications (2)

Publication Number Publication Date
CN105072557A CN105072557A (zh) 2015-11-18
CN105072557B true CN105072557B (zh) 2017-04-19

Family

ID=54501814

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510490058.3A Active CN105072557B (zh) 2015-08-11 2015-08-11 一种三维环绕声重放系统的扬声器环境自适应校准方法

Country Status (1)

Country Link
CN (1) CN105072557B (zh)

Families Citing this family (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106303843B (zh) * 2016-07-29 2018-04-03 北京工业大学 一种多区域不同语音声源的2.5d重放方法
CN107219512B (zh) * 2017-03-29 2020-05-22 北京大学 一种基于声传递函数的声源定位方法
CN107147975B (zh) * 2017-04-26 2019-05-14 北京大学 一种面向不规则扬声器摆放的Ambisonics匹配投影解码方法
US11044552B2 (en) * 2017-08-31 2021-06-22 Harman International Industries, Incorporated Acoustic radiation control method and system
CN109756683B (zh) * 2017-11-02 2024-06-04 深圳市裂石影音科技有限公司 全景音视频录制方法、装置、存储介质和计算机设备
CN108012214B (zh) * 2017-11-08 2019-05-10 西北工业大学 基于广义极小极大凹罚函数的声场重构方法
CN108419174B (zh) * 2018-01-24 2020-05-22 北京大学 一种基于扬声器阵列的虚拟听觉环境可听化实现方法及系统
CN109379694B (zh) * 2018-11-01 2020-08-18 华南理工大学 一种多通路三维空间环绕声的虚拟重放方法
CN111526455A (zh) * 2020-05-21 2020-08-11 菁音电子科技(上海)有限公司 车载音响的校正增强方法及系统
CN113314129B (zh) * 2021-04-30 2022-08-05 北京大学 一种适应环境的声场重放空间解码方法
CN114283833A (zh) * 2021-12-24 2022-04-05 北京达佳互联信息技术有限公司 语音增强模型训练方法、语音增强方法、相关设备及介质

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6798889B1 (en) * 1999-11-12 2004-09-28 Creative Technology Ltd. Method and apparatus for multi-channel sound system calibration
CN104041074A (zh) * 2011-11-11 2014-09-10 汤姆逊许可公司 处理用于产生声场的高保真度立体声响复制表示的刚性球上的球形麦克风阵列的信号的方法和装置
CN104240695A (zh) * 2014-08-29 2014-12-24 华南理工大学 一种优化的基于耳机重放的虚拟声合成方法
CN104581604A (zh) * 2013-10-17 2015-04-29 奥迪康有限公司 再现声学声场的方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005286828A (ja) * 2004-03-30 2005-10-13 Victor Co Of Japan Ltd オーディオ再生装置
JP2012010011A (ja) * 2010-06-23 2012-01-12 Kyoto Univ 境界音場制御の原理を用いた音場再現システム
JP5668765B2 (ja) * 2013-01-11 2015-02-12 株式会社デンソー 車載音響装置

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6798889B1 (en) * 1999-11-12 2004-09-28 Creative Technology Ltd. Method and apparatus for multi-channel sound system calibration
CN104041074A (zh) * 2011-11-11 2014-09-10 汤姆逊许可公司 处理用于产生声场的高保真度立体声响复制表示的刚性球上的球形麦克风阵列的信号的方法和装置
CN104581604A (zh) * 2013-10-17 2015-04-29 奥迪康有限公司 再现声学声场的方法
CN104240695A (zh) * 2014-08-29 2014-12-24 华南理工大学 一种优化的基于耳机重放的虚拟声合成方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于HRTF的虚拟声源定位系统研究以及ARM平台实现;杨飞然;《万方数据 信号与信息处理》;20101125;第12-16页 *

Also Published As

Publication number Publication date
CN105072557A (zh) 2015-11-18

Similar Documents

Publication Publication Date Title
CN105072557B (zh) 一种三维环绕声重放系统的扬声器环境自适应校准方法
CN102907116B (zh) 用于测量多个扬声器和麦克风阵列的设备和方法
TWI275314B (en) System and method for automatic room acoustic correction in multi-channel audio environments
CN103583054B (zh) 用于产生音频输出信号的装置和方法
Gupta et al. Three-dimensional sound field reproduction using multiple circular loudspeaker arrays
TWI512720B (zh) 用以產生多個參數式音訊串流之裝置及方法和用以產生多個揚聲器信號之裝置及方法
US9398393B2 (en) Aural proxies and directionally-varying reverberation for interactive sound propagation in virtual environments
Amengual Garí et al. Optimizations of the spatial decomposition method for binaural reproduction
CN107147975A (zh) 一种面向不规则扬声器摆放的Ambisonics匹配投影解码方法
Frank Source width of frontal phantom sources: Perception, measurement, and modeling
CN1849844B (zh) 确定声场的表示的系统和方法
Borra et al. Soundfield reconstruction in reverberant environments using higher-order microphones and impulse response measurements
Hacihabiboğlu et al. Panoramic recording and reproduction of multichannel audio using a circular microphone array
Zhang et al. Distance-dependent modeling of head-related transfer functions
CN113314129B (zh) 一种适应环境的声场重放空间解码方法
Ge et al. Partially matching projection decoding method evaluation under different playback conditions
Chaman et al. Multipath-enabled private audio with noise
Wendt et al. Perceptual and room acoustical evaluation of a computational efficient binaural room impulse response simulation method
Liu et al. Cocktails, but no party: Multipath-enabled private audio
Diepold et al. HRTF Measurements with Recorded Reference Signal
Asakura et al. Development of a simulation system of the audiovisual environment using the headphone and head-mounted display
US11381927B2 (en) System and method for spatial processing of soundfield signals
Song et al. Simulation of realistic background noise using multiple loudspeakers
Giurda et al. Evaluation of an ILD-based hearing device algorithm using Virtual Sound Environments
Dickins et al. Validation of a Practical Spatial Soundfield Reproduction System Using a Directional Microphone

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