CN108735231B - 主旋律音高序列估计方法 - Google Patents

主旋律音高序列估计方法 Download PDF

Info

Publication number
CN108735231B
CN108735231B CN201810390549.4A CN201810390549A CN108735231B CN 108735231 B CN108735231 B CN 108735231B CN 201810390549 A CN201810390549 A CN 201810390549A CN 108735231 B CN108735231 B CN 108735231B
Authority
CN
China
Prior art keywords
pitch
main melody
frame
melody
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.)
Expired - Fee Related
Application number
CN201810390549.4A
Other languages
English (en)
Other versions
CN108735231A (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.)
Dalian Minzu University
Original Assignee
Dalian Minzu 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 Dalian Minzu University filed Critical Dalian Minzu University
Priority to CN201810390549.4A priority Critical patent/CN108735231B/zh
Publication of CN108735231A publication Critical patent/CN108735231A/zh
Application granted granted Critical
Publication of CN108735231B publication Critical patent/CN108735231B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G10MUSICAL INSTRUMENTS; ACOUSTICS
    • G10LSPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L25/00Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00
    • G10L25/48Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 specially adapted for particular use
    • G10L25/51Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 specially adapted for particular use for comparison or discrimination
    • G10L25/54Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 specially adapted for particular use for comparison or discrimination for retrieval
    • GPHYSICS
    • G10MUSICAL INSTRUMENTS; ACOUSTICS
    • G10LSPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L25/00Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00
    • G10L25/90Pitch determination of speech signals
    • GPHYSICS
    • G10MUSICAL INSTRUMENTS; ACOUSTICS
    • G10LSPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L25/00Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00
    • G10L25/90Pitch determination of speech signals
    • G10L2025/906Pitch tracking

Landscapes

  • Engineering & Computer Science (AREA)
  • Computational Linguistics (AREA)
  • Signal Processing (AREA)
  • Health & Medical Sciences (AREA)
  • Audiology, Speech & Language Pathology (AREA)
  • Human Computer Interaction (AREA)
  • Physics & Mathematics (AREA)
  • Acoustics & Sound (AREA)
  • Multimedia (AREA)
  • Auxiliary Devices For Music (AREA)

Abstract

主旋律音高序列估计方法,属于音乐信息检索领域,为了解决降低短时强伴奏的影响及八度误差的问题,要点是将每帧主旋律音高序列平滑估计值的上下八度范围作为该帧的主旋律音高的搜索范围,效果是将每个频率点的显著度值定义为各次谐波分量幅度的加权和,进一步降低了八度误差。

Description

主旋律音高序列估计方法
技术领域
本发明属于音乐信息检索领域,涉及一种主旋律音高序列估计方法。
背景技术
音乐信息检索领域中的主旋律提取即从音乐信号中估计出主旋律音高序列,其主要依据是音乐中主旋律成分的显著性和时序连续性。传统的主旋律提取方法大多把主旋律提取问题分解成多音高估计和旋律轮廓构建两个子问题,导致无法同时兼顾两个特性,在多音高估计阶段没有考虑到连续帧音高的高度相关性。Jo等提出采用贝叶斯滤波框架建模主旋律提取问题,并借助粒子滤波方法估计参数,但该方法在粒子滤波之前仍然需要进行多音高估计,且具有较多的参数需要估计(S.Jo,C.D.Yoo,and A.Doucet,Melody trackingbased on sequential Bayesian model,IEEE Journal of Selected Topics in SignalProcessing,2011,5(6):1216–1227.)。
此外,音乐信号非常复杂,具有典型非平稳性的本质特点,即某些短时强伴奏音功率可能超过主旋律成分,且主旋律音高频率分布较广,包含几个八度范围,导致现有方案在主旋律提取中来自于伴奏音的错误音高和八度误差广泛存在。
发明内容
为了解决降低短时强伴奏的影响及八度误差的问题,本发明提出如下方案:一种主旋律音高序列估计方法,包括如下步骤:
(1)将每帧主旋律音高序列平滑估计值的上下八度范围作为该帧的主旋律音高的搜索范围;
(2)计算每帧主旋律音高搜索范围内每个频率点的显著度值S′t(f):
Figure BDA0001643334790000021
其中Nh=10是最大谐波次数,|Xt(hf)|为频率f的第h次谐波幅度;
(3)对每帧的S′t(f)进行归一化,得到St(f):
Figure BDA0001643334790000022
(4)利用动态规划算法得到最终的主旋律音高序列值,动态规划的初始条件为:
D(1,ft,j)=S1(fj) (11)
其中D(t,ft,j)是第t帧取ft,j作为旋律音高的累积代价函数,S1(fj)由公式(10)得到,是第1帧的显著度函数在第j个频点的取值,则求解主旋律音高序列最优解的递归公式为:
D(t,ft,j)=St(fj)+max{D(t-1,ft-1,k)-λd(ft-1,k,ft,j)},t=1,...,Nfrm (12)
其中ft,j为第t帧主旋律音高的可能值,ft-1,k为第t-1帧主旋律音高估计值,d(ft-1,k,ft,j)为ft-1,k和ft,j的半音差,λ为惩罚因子;
不断迭代公式(12)得到最终主旋律音高序列估计。
作为技术方案的补充,λ=0.05。
本发明采用粒子滤波得到主旋律音高序列的粗略估计,然后将其作平滑处理,用以限定动态规划的主旋律音高搜索范围,再由动态规划在限定的范围内获得最终主旋律音高序列估计,降低了搜索范围外的强伴奏影响;本发明在动态规划阶段引入惩罚因子进一步减少了短时强伴奏的影响。
本发明在粒子滤波的似然函数中,引入谐波平滑性因子降低了八度误差;本发明的动态规划阶段在主旋律音高搜索范围内将每个频率点的显著度值定义为各次谐波分量幅度的加权和,进一步降低了八度误差。
具体实施方式
本公开涉及一种音乐主旋律提取方法,用于解决传统主旋律提取方法无法兼顾显著性和时序连续性的问题,同时用于降低音乐信号非平稳性导致的自于伴奏音的错误音高和八度误差。
其方案如下:在粒子滤波框架下,随机掷粒子,将粒子滤波的似然函数建模为各次谐波幅度平方和与谐波平滑性因子的乘积,用logistic分布拟合转移概率函数,由粒子滤波的预测和更新方程逐帧递归完成主旋律音高序列粗略估计,然后对其进行平滑处理,取平滑处理后的每帧音高的上、下八度范围作为该帧的主旋律音高的搜索范围,并于被限定的范围内执行动态规划算法以得到最终主旋律音高序列估计。
该方法将显著性和时序连续性融合到同一框架中;采用粒子滤波和动态规划相结合的两阶段方案,限定最终主旋律音高序列的可能范围,降低了主旋律音高估计中的短时强伴奏音影响;采用粒子滤波中的谐波平滑性因子和动态规划中加权求和的显著度函数,减少了主旋律音高估计中的八度误差。
该方法具体包括如下步骤:
S1、计算音乐信号的常Q变换;
计算音乐音频信号的常Q变换方法:对非平稳的音乐信号进行分帧,然后利用常Q变换计算按照对数分布的各频率点的幅度值,每八度范围具有36个按对数分布的频率点。
S2、获得基于粒子滤波的主旋律音高序列的粗略估计;
获得基于粒子滤波的主旋律音高序列粗略估计的方法:
初始化:
Figure BDA0001643334790000041
其中
Figure BDA0001643334790000042
表示第i个粒子的初始音高
Figure BDA0001643334790000043
服从[110Hz,1200Hz]范围内均匀分布,Np=300为粒子数,
Figure BDA0001643334790000044
表示第i个粒子权重。
迭代求解:
(1)按转移概率,预测第t帧旋律音高频率:
Figure BDA0001643334790000045
根据实际音乐中音高转移概率统计分析结果,将音高转移概率建模为logistic分布,其概率密度函数为:
Figure BDA0001643334790000046
其中x=f0,t/f0,t-1,f0,t为第t帧旋律音高频率,且μ=1.00003,s=0.0055045。
(2)计算粒子权重
Figure BDA0001643334790000047
Figure BDA0001643334790000048
其中P和S分别代表了第t帧第i个粒子的功率和谐波平滑性因子,分别定义为:
Figure BDA0001643334790000051
Figure BDA0001643334790000052
其中Am,t为第t帧第i个粒子音高频率f0,t所对应的第m次谐波分量的幅度,H为最大谐波次数。
(3)归一化粒子权重
Figure BDA0001643334790000053
Figure BDA0001643334790000054
(4)计算有效粒子数量
Figure BDA0001643334790000055
Figure BDA0001643334790000056
Figure BDA0001643334790000057
重新初始化粒子音高
Figure BDA0001643334790000058
跳到步骤(1);否则,跳到步骤(5),其中fmin和fmax分别是主旋律音高的频率下限和上限。
(5)计算当前旋律音高估计
Figure BDA0001643334790000059
Figure BDA00016433347900000510
其中
Figure BDA00016433347900000511
是权值最大的各粒子音高频率的平均值。
(6)对粒子按照归一化权重进行重采样,并将所有重采样后的粒子权重都置为1/Np
(7)令t=t+1,若t≤Nfrm(Nfrm为音频的总帧数),跳到步骤(1);否则,结束迭代过程。
结束(1)-(7)的步骤后,由
Figure BDA0001643334790000061
构成主旋律音高序列粗略估计。
S3、获得主旋律音高序列平滑估计序列;
获得主旋律音高序列平滑估计序列的方法:对前后100毫秒时间内的主旋律音高序列粗略估计结果取平均,作为当前帧的主旋律音高序列平滑估计值,各个帧的平滑估计值构成主旋律音高序列平滑估计序列。
S4、获得最终主旋律音高序列估计;
获得最终主旋律音高序列估计的方法:
(1)将每帧主旋律音高序列平滑估计值的上下八度范围作为该帧的主旋律音高的搜索范围。
(2)计算每帧主旋律音高搜索范围内每个频率点的显著度值S′t(f):
Figure BDA0001643334790000062
其中Nh=10是最大谐波次数,|Xt(hf)|为频率f的第h次谐波幅度。
(3)对每帧的S′t(f)进行归一化,得到St(f):
Figure BDA0001643334790000063
(4)利用动态规划算法得到最终的主旋律音高序列值,动态规划的初始条件为:
D(1,ft,j)=S1(fj) (11)
其中D(t,ft,j)是第t帧取ft,j作为旋律音高的累积代价函数,S1(fj)由公式(10)得到,是第1帧的显著度函数在第j个频点的取值。则求解主旋律音高序列最优解的递归公式为:
D(t,ft,j)=St(fj)+max{D(t-1,ft-1,k)-λd(ft-1,k,ft,j)},t=1,...,Nfrm (12)
其中ft,j为第t帧主旋律音高的可能值,ft-1,k为第t-1帧主旋律音高估计值,d(ft-1,k,ft,j)为ft-1,k和ft,j的半音差,λ=0.05为惩罚因子。其中,ft,j、ft-1,k中的j、k使用不同的符号表示,是由于在迭代到第t帧时,第t-1帧主旋律音高估计值ft-1,k已经得到并为确定值,而第t帧主旋律音高的可能值ft,j还没确定,因而以示区别。
不断迭代公式(12)得到最终主旋律音高序列估计。
本发明将旋律的显著性和平滑性融合到粒子滤波框架中,显著性通过粒子滤波中似然函数的各次谐波幅度平方和方式体现,而平滑性则在粒子滤波中将转移概率拟合成logistic分布的方式实现,这样该框架同时完成多音高估计和旋律轮廓构建。为了减少短时强伴奏影响,本发明通过两个途径解决:(1)采用粒子滤波得到主旋律音高序列的粗略估计,然后将其作平滑处理,用以限定动态规划的主旋律音高搜索范围,再由动态规划在限定的范围内获得最终主旋律音高序列估计,降低了搜索范围外的强伴奏影响;(2)在动态规划阶段引入惩罚因子进一步减少短时强伴奏的影响。为了减少八度误差,本发明通过两个方式解决:(1)在粒子滤波的似然函数中,引入谐波平滑性因子降低八度误差;(2)动态规划阶段在主旋律音高搜索范围内将每个频率点的显著度值定义为各次谐波分量幅度的加权和,进一步降低八度误差。因此,本发明既将显著性和平滑性同时融合到了粒子滤波这一框架中,又降低了短时强伴奏的影响和八度误差。
也即由上述方案,本实施例所述方案具有如下有益效果:将显著性和平滑性融合到粒子滤波框架中,可同时完成多音高估计和旋律轮廓构建;通过粒子滤波和动态规划相结合的方案完成最终的主旋律音高序列估计,并在动态规划中限定主旋律音高搜索范围,还引入惩罚因子,降低了短时强伴奏的影响;通过在粒子滤波中引入谐波平滑性因子和动态规划中加权求和的显著度函数减少了主旋律音高估计中的八度误差。
本实施例所述的音乐主旋律提取方法,其与本公开背景技术中述及的Jo等公开的方案相比较,区别如下:
1.本发明方法和Jo等的方法的似然函数和转移概率密度都不一样,具有不同的音高显著度描述方法和转移概率密度表达式,尤其是概率密度,本发明的方法是在实际音乐音频的基础上作统计分析,将音高转移概率拟合为logistic分布,而Jo等的方法是用常规的高斯密度拟合转移概率;本方法还在似然函数中引入了谐波平滑性因子用以减少八度误差;
2.Jo等的方法先用多音高估计得到各帧的主旋律音高候选,然后采用粒子滤波的方法得到最终的主旋律音高序列。在某些主旋律音符结尾部分,功率可能小于某些强伴奏的场合,Jo等的方法仍无法排除这些错误音高,准确性较差;
3.两个方法的流程不同,本发明的方法是先粒子滤波再动态规划,Jo等的方法是先多音高估计再粒子滤波;
4.本方法具有较少的参数,参数调整容易;而Jo等的方法具有较多的参数,调整困难。
由上述,音乐具有很强的非平稳特性,强伴奏音可能在某些短时间内超过主旋律功率,比如,在某些音符结尾部分,可能主旋律功率小于有些伴奏;此外,八度误差也经常存在,在某个音符演奏过程中,估计音高会因为各次谐波幅度比例的变化而在不同的八度之间切换。因此,若能事先确定主旋律的大致范围,将能滤除范围外的伴奏音高和八度误差。
Jo等的方法先用多音高估计得到各帧的旋律音高候选,然后采用粒子滤波的方法得到最终的主旋律音高序列。该方法仍不能解决上述问题。本发明提出的方法首先用粒子滤波估计出旋律音高,同样也具有上述问题,但是粒子滤波方法估计的结果大部分都是准的,因此对估计序列作平滑处理就削弱了错误估计的影响,仍能确定旋律音高的大致范围,在被限制的范围内再利用动态规划作处理,就能滤除范围外的短时强伴奏音高和八度误差。
以上所述,仅为本发明创造较佳的具体实施方式,但本发明创造的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本发明创造披露的技术范围内,根据本发明创造的技术方案及其发明构思加以等同替换或改变,都应涵盖在本发明创造的保护范围之内。

Claims (2)

1.一种主旋律音高序列估计方法,其特征在于,包括如下步骤:
(1)将每帧主旋律音高序列平滑估计值的上下八度范围作为该帧的主旋律音高的搜索范围,其中计算平滑估计值的方法为:在粒子滤波框架下,随机掷粒子,将粒子滤波的似然函数建模为各次谐波幅度平方和与谐波平滑性因子的乘积,用logistic分布拟合转移概率函数,由粒子滤波的预测和更新方程逐帧递归完成主旋律音高序列粗略估计,对前后100毫秒时间内的主旋律音高序列粗略估计结果取平均,作为当前帧的主旋律音高序列平滑估计值;
(2)计算每帧主旋律音高搜索范围内每个频率点的显著度值S′t(f):
Figure FDA0003155062900000011
其中Nh=10是最大谐波次数,|Xt(hf)|为频率f的第h次谐波幅度;
(3)对每帧的S′t(f)进行归一化,得到St(f):
Figure FDA0003155062900000012
(4)利用动态规划算法得到最终的主旋律音高序列值,动态规划的初始条件为:
D(1,ft,j)=S1(fj) (11)
其中D(t,ft,j)是第t帧取ft,j作为旋律音高的累积代价函数,S1(fj)由公式(10)得到,是第1帧的显著度函数在第j个频点的取值,则求解主旋律音高序列最优解的递归公式为:
D(t,ft,j)=St(fj)+max{D(t-1,ft-1,k)-λd(ft-1,k,ft,j)},t=1,...,Nfrm (12)
其中ft,j为第t帧主旋律音高的可能值,ft-1,k为第t-1帧主旋律音高估计值,d(ft-1,k,ft,j)为ft-1,k和ft,j的半音差,λ为惩罚因子;
不断迭代公式(12)得到最终主旋律音高序列估计。
2.如权利要求1所述的主旋律音高序列估计方法,其特征在于λ=0.05。
CN201810390549.4A 2018-04-27 2018-04-27 主旋律音高序列估计方法 Expired - Fee Related CN108735231B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810390549.4A CN108735231B (zh) 2018-04-27 2018-04-27 主旋律音高序列估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810390549.4A CN108735231B (zh) 2018-04-27 2018-04-27 主旋律音高序列估计方法

Publications (2)

Publication Number Publication Date
CN108735231A CN108735231A (zh) 2018-11-02
CN108735231B true CN108735231B (zh) 2021-11-12

Family

ID=63940045

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810390549.4A Expired - Fee Related CN108735231B (zh) 2018-04-27 2018-04-27 主旋律音高序列估计方法

Country Status (1)

Country Link
CN (1) CN108735231B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109979488B (zh) * 2019-03-14 2021-05-14 浙江大学 基于重音分析的人声转乐谱系统
CN112735365A (zh) * 2020-12-24 2021-04-30 重庆悦派对信息科技有限公司 一种基于概率模型的主旋律自动提取算法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20060019096A (ko) * 2004-08-26 2006-03-03 주식회사 케이티 허밍 기반의 음원 질의/검색 시스템 및 그 방법
CN103412886A (zh) * 2013-07-18 2013-11-27 北京航空航天大学 一种基于音高序列的音乐旋律匹配方法
CN103823867A (zh) * 2014-02-26 2014-05-28 深圳大学 一种基于音符建模的哼唱式音乐检索方法及系统
CN107256710A (zh) * 2017-08-01 2017-10-17 中国农业大学 一种基于动态时间伸缩算法的哼唱旋律识别方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20060019096A (ko) * 2004-08-26 2006-03-03 주식회사 케이티 허밍 기반의 음원 질의/검색 시스템 및 그 방법
CN103412886A (zh) * 2013-07-18 2013-11-27 北京航空航天大学 一种基于音高序列的音乐旋律匹配方法
CN103823867A (zh) * 2014-02-26 2014-05-28 深圳大学 一种基于音符建模的哼唱式音乐检索方法及系统
CN107256710A (zh) * 2017-08-01 2017-10-17 中国农业大学 一种基于动态时间伸缩算法的哼唱旋律识别方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Melody Extraction From Polyphonic Music Using Particle Filter and Dynamic Programming;Weiwei Zhang et al.;《 IEEE/ACM Transactions on Audio, Speech, and Language Processing》;20180930;第26卷(第9期);第1621-1629页 *
Melody Tracking Based on Sequential Bayesian Model;Seokhwan Jo et al.;《 IEEE Journal of Selected Topics in Signal Processing》;20111031;第5卷(第6期);第1216-1223页 *
Probabilistic model for main melody extraction using Constant-Q transform;Benoit Fuentes et al.;《2012 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)》;20120831;第5357-5360页 *
结合改进欧几里得算法和动态规划的音乐主旋律提取;张维维等;《信号处理》;20180825;第1009-1013页 *

Also Published As

Publication number Publication date
CN108735231A (zh) 2018-11-02

Similar Documents

Publication Publication Date Title
WO2020177371A1 (zh) 一种用于数字助听器的环境自适应神经网络降噪方法、系统及存储介质
CN111261183B (zh) 一种语音去噪的方法及装置
US20050060153A1 (en) Method and appratus for speech characterization
CN104685562B (zh) 用于从嘈杂输入信号中重构目标信号的方法和设备
US9754608B2 (en) Noise estimation apparatus, noise estimation method, noise estimation program, and recording medium
CN110070880B (zh) 用于分类的联合统计模型的建立方法及应用方法
WO2005066927A1 (ja) 多重音信号解析方法
Pertusa et al. Multiple fundamental frequency estimation using Gaussian smoothness
Shahnaz et al. Pitch estimation based on a harmonic sinusoidal autocorrelation model and a time-domain matching scheme
DE602004000716T2 (de) Rauschadaptierungssystem für ein Sprachmodell, Verfahren zur Rauschadaptierung und Programm zur Rauschadaptierung für Spracherkennung
EP1891624A2 (en) Multi-sensory speech enhancement using a speech-state model
CN108735231B (zh) 主旋律音高序列估计方法
CN108647596B (zh) 音乐主旋律提取方法
CN109192200A (zh) 一种语音识别方法
JP4512848B2 (ja) 雑音抑圧装置及び音声認識システム
JP5344251B2 (ja) 雑音除去システム、雑音除去方法および雑音除去プログラム
JP5351856B2 (ja) 音源パラメータ推定装置と音源分離装置とそれらの方法と、プログラムと記憶媒体
CN108628982B (zh) 基于粒子滤波的主旋律音高序列的粗略估计方法
CN108536871B (zh) 粒子滤波并限定动态规划搜索范围的音乐主旋律提取方法及装置
CN108595648B (zh) 音乐主旋律提取系统
Dionelis et al. Modulation-domain Kalman filtering for monaural blind speech denoising and dereverberation
JP4673828B2 (ja) 音声信号区間推定装置、その方法、そのプログラム及び記録媒体
Stouten et al. Joint removal of additive and convolutional noise with model-based feature enhancement
CN108573698B (zh) 一种基于性别融合信息的语音降噪方法
Pandiaraj et al. A confidence measure based—Score fusion technique to integrate MFCC and pitch for speaker verification

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
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20211112