CN108134604A - 一种基于直接法求解瞬时频率的毛刺定位与消除方法 - Google Patents
一种基于直接法求解瞬时频率的毛刺定位与消除方法 Download PDFInfo
- Publication number
- CN108134604A CN108134604A CN201711344919.2A CN201711344919A CN108134604A CN 108134604 A CN108134604 A CN 108134604A CN 201711344919 A CN201711344919 A CN 201711344919A CN 108134604 A CN108134604 A CN 108134604A
- Authority
- CN
- China
- Prior art keywords
- instantaneous
- signal
- phase
- local
- instantaneous frequency
- 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
- 238000000034 method Methods 0.000 title claims abstract description 47
- 230000003121 nonmonotonic effect Effects 0.000 claims abstract description 19
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 16
- 238000005070 sampling Methods 0.000 claims description 6
- 230000006870 function Effects 0.000 description 18
- 230000009466 transformation Effects 0.000 description 8
- 238000005516 engineering process Methods 0.000 description 3
- 238000001914 filtration Methods 0.000 description 3
- 230000009471 action Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 235000013399 edible fruits Nutrition 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 230000000877 morphologic effect Effects 0.000 description 2
- QHGUCRYDKWKLMG-UHFFFAOYSA-N octopamine Chemical compound NCC(O)C1=CC=C(O)C=C1 QHGUCRYDKWKLMG-UHFFFAOYSA-N 0.000 description 2
- 229960001576 octopamine Drugs 0.000 description 2
- 241000208340 Araliaceae Species 0.000 description 1
- 206010021703 Indifference Diseases 0.000 description 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 1
- 235000003140 Panax quinquefolius Nutrition 0.000 description 1
- 241000219000 Populus Species 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000000295 complement effect Effects 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 238000011089 mechanical engineering Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 230000008447 perception Effects 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000012731 temporal analysis Methods 0.000 description 1
- 238000000700 time series analysis Methods 0.000 description 1
- 230000001052 transient effect Effects 0.000 description 1
- 238000002077 transients frequency Methods 0.000 description 1
Classifications
-
- H—ELECTRICITY
- H03—ELECTRONIC CIRCUITRY
- H03L—AUTOMATIC CONTROL, STARTING, SYNCHRONISATION OR STABILISATION OF GENERATORS OF ELECTRONIC OSCILLATIONS OR PULSES
- H03L7/00—Automatic control of frequency or phase; Synchronisation
- H03L7/06—Automatic control of frequency or phase; Synchronisation using a reference signal applied to a frequency- or phase-locked loop
- H03L7/08—Details of the phase-locked loop
- H03L7/085—Details of the phase-locked loop concerning mainly the frequency- or phase-detection arrangement including the filtering or amplification of its output signal
-
- H—ELECTRICITY
- H04—ELECTRIC COMMUNICATION TECHNIQUE
- H04L—TRANSMISSION OF DIGITAL INFORMATION, e.g. TELEGRAPHIC COMMUNICATION
- H04L27/00—Modulated-carrier systems
- H04L27/32—Carrier systems characterised by combinations of two or more of the types covered by groups H04L27/02, H04L27/10, H04L27/18 or H04L27/26
- H04L27/34—Amplitude- and phase-modulated carrier systems, e.g. quadrature-amplitude modulated carrier systems
- H04L27/36—Modulator circuits; Transmitter circuits
Landscapes
- Engineering & Computer Science (AREA)
- Computer Networks & Wireless Communication (AREA)
- Signal Processing (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
- Complex Calculations (AREA)
Abstract
一种基于直接法求解瞬时频率的毛刺定位与消除方法，先使用局部均值分解方法将信号分解为乘积函数分量，从而得到乘积函数的调频信号，将调频信号的所有局部极大值点置为1，所有的局部极小值点置为‑1；再将调频信号取反余弦，得到非单调的瞬时相位，对非单调的瞬时相位进行展开，得到单调递增的瞬时相位；然后对单调递增的瞬时相位进行前向差分，得到调频信号的带有毛刺的瞬时频率；再确定带有毛刺的瞬时频率的毛刺位置；然后将带有毛刺的瞬时频率在毛刺位置序列处的数据删除，并使用三次样条插补出毛刺位置序列处的瞬时频率值，最终得到无毛刺的瞬时频率；本发明具有无参数，不依赖使用者经验的优点，具有局部处理的特点。
Description
技术领域
本发明属于信号处理技术领域，具体涉及一种基于直接法求解瞬时频率的毛刺定位与消除方法。
背景技术
瞬时频率是调幅-调频(AM-FM)信号的重要调制参数，对于探究非平稳、非线性过程的详细机制具有重要意义。目前，瞬时频率的求解思路是先将信号分解成为一系列单分量信号，然后对这些单分量信号进行解调。比较常用的包括Hilbert-Huang变换(NordenE.Huang.etc.The empirical mode decomposition and the Hilbert spectrum fornonlinear and non-stationary time series analysis[J].Proceedings of the royalsociety A:Mathematical,Physical and Engineering Sciences.1998:pp903-995.)，经验AM-FM分解(Norden E.Huang et al,on instantaneous frequency,Advances inAdaptive Data Analysis 1(2),2009:pp177-229.)，局部均值分解(Jonathan S.Smith,The local mean decomposition and its application to EEG perception data,J.R.Soc.Interface 2(5),2005:pp443–454.)三种方式。Hilbert-Huang变换是首先使用经验模式分解将信号分解为一系列本征模态函数，然后使用希尔伯特变换求取本征模态函数的瞬时频率；经验AM-FM分解和Hilbert-Huang变换一样首先将信号分解为本征模态函数，然后将本征模式函数分解为调频部分的包络函数和调幅部分的纯调频信号，然后使用直接正交法对调频部分的纯调频信号求解瞬时频率；局部均值分解是将信号分解成为一系列的乘积函数，每一乘积函数是由调幅部分的包络函数和调频部分的纯调频信号组成，使用直接法求取乘积函数的纯调频信号的瞬时频率。
直接法是一种性能比较好的纯调频信号求解方法，相较于Hilbert变换、Teager能量法等其他求解方法，具有无端点效应、求解速度快、不会出现负频率的优点。但是，在对纯调频信号进行反余弦时，由于反余弦函数是定义在区间[-1,1]上的单调函数，而纯调频信号则是非单调函数，那么对纯调频信号求反余弦后，得到的相位函数就是非单调的。根据瞬时频率的定义，瞬时频率应当是相位的导数，且瞬时频率应当都是正的，于是，将非单调的相位函数展开为单调递增的瞬时相位是获取瞬时频率的关键一步。文献(任达千,杨世锡,吴昭同,等.信号瞬时频率直接计算方法与Hilbert变换及Teager能量法比较[J].机械工程学报,2013,49(9):42-48.)推导出了直接法的相位展开方法。直接法存在的一个重要问题是，其得到的瞬时频率曲线在极值点附近存在较大误差，形成一系列毛刺状突变点(以下简称毛刺)，给直接法的应用带来了一定的局限。
为了克服使用直接法求取的瞬时频率在信号极值点位置附近产生的毛刺，许多学者提出了一些方案。因为毛刺往往对应于IF曲线中的高频成分，因此有学者(任达千,杨世锡,吴昭同,等.信号瞬时频率直接计方法与Hilbert变换及Teager能量法比较[J].机械工程学报,2013,49(9):42-48.)使用滑动平均法对得到的IF进行平滑，但是未给出包括平滑次数、平滑步长等平滑参数的选择方式；亦有学者使用形态学滤波法(唐贵基，王晓龙.LMD的LabVIEW实现及其在齿轮故障诊断中的应用[J].中国测试，2014.40(1):p101-105.)来对直接法求取的瞬时频率中的畸变进行消除，同样未给出使用的结构元素、操作类型等参数，且这些方法都是对整个瞬时频率序列进行的无差别处理，在处理毛刺位置的同时，也对非毛刺位置的正常数值产生了影响，目前未见可以精确定位这些毛刺点并将之定点消除的方法。
发明内容
为了克服上述现有技术的缺点，本发明的目的在于提供一种基于直接法求解瞬时频率的毛刺定位与消除方法，可以精确定位毛刺点并将之定点消除，然后使用三次样条插值补全被消除的数据，使得所求的瞬时频率结果保证在非极值点位置处数值不被改变的同时，较为精准的求出毛刺位置处的瞬时频率值。
为了达到上述目的，本发明采取的技术方案是：
一种基于直接法求解瞬时频率的毛刺定位与消除方法，包括以下步骤：
1)使用局部均值分解方法将原始信号x(t)分解为乘积函数分量，从而得到乘积函数的调频信号s(t)，将调频信号s(t)所有的局部极大值点MaxEi置为1，所有的局部极小值点MinEi置为-1，且满足-1≤s(t)≤1；
2)将调频信号s(t)取反余弦，得到非单调的瞬时相位
3)对非单调的瞬时相位进行展开，得到单调递增的瞬时相位其展开方法为：将非单调的瞬时相位的所有局部极值点进行编号；如果第一个局部极值点是局部极大值点，那么，如果局部极值点的编号i为奇数，第i极值点与第i+1极值点之间的相位重新定义为如果i为偶数，第i与i+1个极值点之间的相位重新定义为如果第一个局部极值点是局部极小值点，那么，如果局部极值点的编号i为奇数，第i极值点与第i+1极值点之间的相位重新定义为如果i为偶数，第i与i+1个极值点之间的相位重新定义为最终得到单调递增的瞬时相位
4)对单调递增的瞬时相位进行前向差分，设调频信号s(t)的采样频率为fs，得到调频信号s(t)的带有毛刺的瞬时频率f0(t)：
上式中diff表示前向差分运算；
5)确定带有毛刺的瞬时频率f0(t)的毛刺位置：确定非单调的瞬时相位的局部极值点E，针对第i个局部极值点为Ei，将Ei位置及其两侧各2个采样点位置共计5个点按幅值的绝对值大小进行降序排序，分别标识为D1～D5，其对应采样时刻分别为
设如果Z＜0，那么，若则将与存储至毛刺位置序列T中，若则将与存储至毛刺位置序列T中；如果Z＞0，若则将与存储至毛刺位置序列T中，若则将与存储至毛刺位置序列T中；最终得到毛刺位置序列T；
6)将步骤4)得到的带有毛刺的瞬时频率f0(t)在步骤5)得到的毛刺位置序列T处的数据删除，并使用三次样条插补出毛刺位置序列T处的瞬时频率值，最终得到无毛刺的瞬时频率f(t)。
本发明的有益效果为：
(1)本发明与现有的毛刺去除方法相比，本发明在确定毛刺位置点后，只剔除瞬时频率中的毛刺点，具有局部处理的特性；而滑动平均、形态学率滤波等传统方式是全局性的处理方式，经过传统方式之后，原本的非毛刺位置的瞬时频率数值也被改变；
(2)传统的滑动平均法需要设定步长、滑动平均次数，形态学滤波需要设定结构元素与操作类型，这些参数的选择都有赖于使用者的经验积累，而本发明提出的方法是无参数的，不依赖经验的，因此效果是较为客观的。
附图说明
图1是本发明方法流程图。
图2是实施例待分解的原始信号x(t)的波形图。
图3是实施例使用局部均值分解对原始信号x(t)进行分解得到的分解结果。
图4是实施例原始信号x(t)第1个乘积函数的调频部分经过局部极大值点置1，局部极小值点置-1的操作之后得到的纯调频信号s(t)的波形图。
图5是实施例将纯调频信号s(t)取反余弦，得到非单调的瞬时相位
图6是实施例将非单调的瞬时相位进行展开，得到单调递增的瞬时相位
图7是实施例纯调频信号s(t)的带有毛刺的瞬时频率f0(t)。
图8是实施例将单调递增的瞬时相位中的毛刺定位出来并使用“*”标识。
图9为实施例无毛刺的瞬时频率f(t)。
图10是使用Hilbert变换求解的调频信号s(t)的瞬时频率。
具体实施方式
以下结合附图和实施例作对本发明进一步的详细说明。
参照图1，一种基于直接法求解瞬时频率的毛刺定位与消除方法，包括以下步骤：
1)产生一个仿真信号作为原始信号x(t)，公式为：
x(t)＝(1+0.5cos(10πt))cos(100πet+2cos(10πt))+(1+0.3sin(5πt))sin(30πt)+sin(10πt)
设定采样频率为1024Hz，t∈[0,1]，其波形如图2所示；将原始信号x(t)使用局部均值分解方法分解为4个乘积函数(Product Function,PF)和1个余项R，结果如图3所示；本实施例取第1个PF的调频信号s1(t)为分析对象，按照直接法规则，将调频信号s1(t)的全部局部极大值点置为1，将调频信号s1(t)的全部局部极小值点置为-1，得到纯调频信号s(t)，其波形如图4所示；
2)将调频信号s(t)取反余弦，得到非单调的瞬时相位如图5所示；
3)对非单调的瞬时相位进行展开，得到单调递增的瞬时相位其展开方法为：将非单调的瞬时相位的所有局部极值点进行编号；如果第一个局部极值点是局部极大值点，那么，如果局部极值点的编号i为奇数，第i极值点与第i+1极值点之间的相位重新定义为如果i为偶数，第i与i+1个极值点之间的相位重新定义为如果第一个局部极值点是局部极小值点，那么，如果局部极值点的编号i为奇数，第i极值点与第i+1极值点之间的相位重新定义为如果i为偶数，第i与i+1个极值点之间的相位重新定义为最终得到单调递增的瞬时相位本实施例的单调递增的瞬时相位如图6所示；
4)对单调递增的瞬时相位进行前向差分，调频信号s(t)的采样频率为fs＝1024Hz，得到调频信号s(t)的带有毛刺的瞬时频率f0(t)，如图7所示；
上式中diff表示前向差分运算；
5)确定带有毛刺的瞬时频率f0(t)的毛刺位置：确定非单调的瞬时相位的局部极值点E，针对第i个局部极值点为Ei，将Ei位置及其两侧各2个采样点位置共计5个点按幅值的绝对值大小进行降序排序，分别标识为D1～D5，其对应采样时刻分别为
设如果Z＜0，那么，若则将与存储至毛刺位置序列T中，若则将与存储至毛刺位置序列T中；如果Z＞0，若则将与存储至毛刺位置序列T中，若则将与存储至毛刺位置序列T中；最终得到毛刺位置序列T，使用“*”在f0(t)中将毛刺点位置标出，如图8所示；
6)将步骤4)得到的带有毛刺的瞬时频率f0(t)在步骤5)得到的毛刺位置序列T处的数据删除，并使用三次样条插补出毛刺位置序列T处的瞬时频率值，最终可以得到无毛刺的瞬时频率f(t)，如图9所示；对比图9和图8，发现图8中杂乱的毛刺被悉数去除。
为说明本方法求解的客观性和准确性，Hilbert变换是目前学术界较为常用的瞬时频率求解方法，因此在本实施例中，使用Hilbert变换对纯调频信号s(t)进行求解，得到其瞬时频率如图10所示；对比图9和图10，发现两者形状差异不大，说明本发明并没有波及到非毛刺位置数值，实现了对调频信号s(t)的正确解调。
Claims (1)
1.一种基于直接法求解瞬时频率的毛刺定位与消除方法，其特征在于，包括以下步骤：
1)使用局部均值分解方法将原始信号x(t)分解为乘积函数分量，从而得到乘积函数的调频信号s(t)，将调频信号s(t)所有的局部极大值点MaxEi置为1，所有的局部极小值点MinEi置为-1，且满足-1≤s(t)≤1；
2)将调频信号s(t)取反余弦，得到非单调的瞬时相位
3)对非单调的瞬时相位进行展开，得到单调递增的瞬时相位其展开方法为：将非单调的瞬时相位的所有局部极值点进行编号；如果第一个局部极值点是局部极大值点，那么，如果局部极值点的编号i为奇数，第i极值点与第i+1极值点之间的相位重新定义为如果i为偶数，第i与i+1个极值点之间的相位重新定义为如果第一个局部极值点是局部极小值点，那么，如果局部极值点的编号i为奇数，第i极值点与第i+1极值点之间的相位重新定义为如果i为偶数，第i与i+1个极值点之间的相位重新定义为最终得到单调递增的瞬时相位
4)对单调递增的瞬时相位进行前向差分，设调频信号s(t)的采样频率为fs，得到调频信号s(t)的带有毛刺的瞬时频率f0(t)：
上式中diff表示前向差分运算；
5)确定带有毛刺的瞬时频率f0(t)的毛刺位置：确定非单调的瞬时相位的局部极值点E，针对第i个局部极值点为Ei，将Ei位置及其两侧各2个采样点位置共计5个点按幅值的绝对值大小进行降序排序，分别标识为D1～D5，其对应采样时刻分别为
设如果Z＜0，那么，若则将与存储至毛刺位置序列T中，若则将与存储至毛刺位置序列T中；如果Z＞0，若则将与存储至毛刺位置序列T中，若则将与存储至毛刺位置序列T中；最终得到毛刺位置序列T；
6)将步骤4)得到的带有毛刺的瞬时频率f0(t)在步骤5)得到的毛刺位置序列T处的数据删除，并使用三次样条插补出毛刺位置序列T处的瞬时频率值，最终得到无毛刺的瞬时频率f(t)。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711344919.2A CN108134604B (zh) | 2017-12-15 | 2017-12-15 | 一种基于直接法求解瞬时频率的毛刺定位与消除方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711344919.2A CN108134604B (zh) | 2017-12-15 | 2017-12-15 | 一种基于直接法求解瞬时频率的毛刺定位与消除方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108134604A true CN108134604A (zh) | 2018-06-08 |
| CN108134604B CN108134604B (zh) | 2021-01-15 |
Family
ID=62389489
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201711344919.2A Active CN108134604B (zh) | 2017-12-15 | 2017-12-15 | 一种基于直接法求解瞬时频率的毛刺定位与消除方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108134604B (zh) |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6990436B1 (en) * | 2003-11-28 | 2006-01-24 | The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration | Computing frequency by using generalized zero-crossing applied to intrinsic mode functions |
| CN101869477A (zh) * | 2010-05-14 | 2010-10-27 | 北京工业大学 | 一种自适应脑电信号中眼电伪迹的自动去除方法 |
| CN102508120A (zh) * | 2011-10-18 | 2012-06-20 | 中国民航大学 | 利用局部均值分解实现飞机导线微弱故障诊断与定位方法 |
| CN102904847A (zh) * | 2012-09-17 | 2013-01-30 | 重庆大学 | 一种双重消噪的基于瞬时信息的信号调制识别方法 |
| CN106932641A (zh) * | 2016-10-24 | 2017-07-07 | 南京理工大学 | 一种基于改进hht变换的电网故障诊断方法 |
-
2017
- 2017-12-15 CN CN201711344919.2A patent/CN108134604B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6990436B1 (en) * | 2003-11-28 | 2006-01-24 | The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration | Computing frequency by using generalized zero-crossing applied to intrinsic mode functions |
| CN101869477A (zh) * | 2010-05-14 | 2010-10-27 | 北京工业大学 | 一种自适应脑电信号中眼电伪迹的自动去除方法 |
| CN102508120A (zh) * | 2011-10-18 | 2012-06-20 | 中国民航大学 | 利用局部均值分解实现飞机导线微弱故障诊断与定位方法 |
| CN102904847A (zh) * | 2012-09-17 | 2013-01-30 | 重庆大学 | 一种双重消噪的基于瞬时信息的信号调制识别方法 |
| CN106932641A (zh) * | 2016-10-24 | 2017-07-07 | 南京理工大学 | 一种基于改进hht变换的电网故障诊断方法 |
Non-Patent Citations (3)
| Title |
|---|
| QING LI: "Revision of bearing fault characteristic spectrum using LMD and interpolation correction algorithm", 《9TH INTERNATIONAL CONFERENCE ON DIGITAL ENTERPRISE TECHNOLOGY- DET2016 – "INTELLIGENT MANUFACTURING IN》 * |
| 任达千: "信号瞬时频率直接计算法与Hilbert变换及", 《机械工程学报》 * |
| 任达千: "基于LMD 的信号瞬时频率求取方法及实验", 《浙江大学学报》 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108134604B (zh) | 2021-01-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Lin | Power harmonics and interharmonics measurement using recursive group-harmonic power minimizing algorithm | |
| Zhao et al. | Power spectrum estimation of the welch method based on imagery EEG | |
| CN106405230B (zh) | 频率测量方法和装置 | |
| CN103926462B (zh) | 一种电力系统谐波快速分析方法及运行装置 | |
| CN104898427B (zh) | 电机伺服系统的分段滤波迭代学习控制方法 | |
| CN106610450B (zh) | 一种单相电能计量芯片 | |
| CN109241915B (zh) | 基于振动信号平稳非平稳性判别及特征甄别的智能电厂泵机故障诊断方法 | |
| CN105137373B (zh) | 一种指数信号的去噪方法 | |
| CN101320060A (zh) | 快速相位计 | |
| CN106226590A (zh) | 一种电力系统同步相量测量方法 | |
| CN109001532A (zh) | 基于fpga的模拟信号fft实现方法及其电路 | |
| CN101701985A (zh) | 定频变点电网谐波检测方法及其测量仪 | |
| CN108134604A (zh) | 一种基于直接法求解瞬时频率的毛刺定位与消除方法 | |
| CN108090270B (zh) | 一种基于形态学滤波和盲源分离的暂态振荡参数识别方法 | |
| DE112013001232T5 (de) | Verfahren und System zum Charakterisieren einer Frequenzübersetzungsvorrichtung | |
| CN205453724U (zh) | 一种多功能通用混沌电路系统 | |
| Zhou et al. | Research on rough and detailed combination spectrum analysis of electromagnetic signals based on complex modulation ZoomFFT | |
| Osipov et al. | Algorithms of packet wavelet transform for power determination under nonsinusoidal modes | |
| CN110286268A (zh) | 基于频率波动的电力波形畸变评估方法 | |
| Pigazo et al. | Accurate and computationally efficient implementation of the IEEE 1459-2000 standard in three-phase three-wire power systems | |
| CN204302410U (zh) | 一种智能波形参数自动测试分析仪 | |
| CN108957130A (zh) | 一种基于瞬时无功功率理论的电力系统谐波检测方法 | |
| CN107863963B (zh) | 一种适用于非连续可导模拟信号的采样与重构方法 | |
| CN108932373A (zh) | 一种基于虚拟仪器技术的振荡器相位噪声分析方法 | |
| CN106771678B (zh) | 一种基于希尔伯特-黄变换和专家系统的鉴相方法及装置 |
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 |