CN105372492A - 基于三条dft复数谱线的信号频率测量方法 - Google Patents
基于三条dft复数谱线的信号频率测量方法 Download PDFInfo
- Publication number
- CN105372492A CN105372492A CN201410435495.0A CN201410435495A CN105372492A CN 105372492 A CN105372492 A CN 105372492A CN 201410435495 A CN201410435495 A CN 201410435495A CN 105372492 A CN105372492 A CN 105372492A
- Authority
- CN
- China
- Prior art keywords
- frequency
- spectral line
- signal
- discrete
- spectral lines
- 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
Landscapes
- Complex Calculations (AREA)
- Measurement Of Resistance Or Impedance (AREA)
- Measuring Frequencies, Analyzing Spectra (AREA)
Abstract
本发明涉及一种基于三条DFT复数谱线的信号频率测量方法,属于信号参数测量技术领域。本发明的特征在于其处理步骤包含:采样信号经过加窗处理后进行DFT变换,在对应频率范围内通过比较谱线幅值或者实部虚部绝对值之和搜索最高谱线及其相邻的两条谱线,基于三条谱线对应复数的实部虚部的运算求取中间参数,再通过求解方程、或反函数公式、或逼近多项式获得频偏参数,最终测量出信号频率。本发明在搜索最高谱线、以及计算中间参数时减少了乘法和开方计算,简化了信号频率测量算法的实现;同时采用相邻谱线取差的形式,抑制了其他信号旁瓣的影响,提高了测量准确度。
Description
技术领域
本发明涉及一种基于三条DFT复数谱线的信号频率测量方法,属于信号参数测量技术领域。
背景技术
当前,基于离散傅里叶变换DFT或其快速算法FFT分析频率信号的方法已经广泛使用。但是,DFT具有栏栅效应,即实际信号频率未必落在离散谱线上,由此需要采用插值算法估计实际信号的频率。2012年《中国电机工程学报》32卷16期上发表的“基于三谱线插值FFT 的电力谐波分析算法”文章中提出了对输入离散信号加窗傅里叶变换后,通过采用幅值最高谱线及其两边相邻的谱线,三条谱线共同插值测量信号频率的方法。如果三条谱线的离散频率序号分别对应k 1-1、k 1和k 1+1,三条谱线的幅值|Y 1 |=|Y(k 1-1)|、|Y 2|=|Y(k 1)|和|Y 3|=|Y(k 1+1)|满足:
,
其中,γ= k 0- k 1,k 0对应实际信号频率。当N较大时,上式可以简化表示为λ=g(γ),其反函数记为。该文进一步提出采用多项式逼近方法计算:
。
已有方法的不足在于基于谱线幅值进行最高谱线搜索。而且,获得谱线幅值需要计算实部和虚部的平方和、然后进行开方,该计算量较大。已有方法对中间参数γ的计算也基于谱线幅值实现,同时虽然采用了多条谱线,并未有效抑制其他频点信号的旁瓣干扰,测量准确度低。
发明内容
本发明的目的是在搜索最高谱线、以及计算中间参数λ时减少乘法和开方计算,减少运算量;同时,设计更合理的中间参数γ的计算公式,实现旁瓣干扰的抑制,提高测量准确度。
本发明为解决上述技术问题而提供一种基于三条DFT复数谱线的信号频率测量方法,该测量方法的步骤如下:
步骤(1):将采样率为FS 、采样点为连续截取的采样信号x(n),进行加窗处理得到加窗信号y(n),加窗处理公式为:
,
其中为N点的窗函数序列,n = 0:(N-1);
步骤(2):对加窗信号y(n)进行离散傅里叶DFT变换,得到离散频谱Y(k),其中离散频率序号k = 0:(N-1);
步骤(3):在设定频率范围所对应的离散频率序号范围[ks ,
ke ]内,搜索最大的谱线作为最高谱线,其离散频率序号为k 1,同时选择该谱线两侧、离散频率序号为k 1-1和k 1+1的相邻两条谱线,其中Re(Z)表示取复数Z的实部值,Im(Z)表示取复数Z的虚部值;
步骤(4):依据k 1-1、k 1和k 1+1对应的三条谱线Y 1 =Y(k 1-1)、Y 2 =Y(k 1)和Y 3 =Y(k 1+1)计算中间参数λ:
步骤(5):求解如下方程中的频偏参数γ:
,
其中,W(ω)是窗函数w(n)的离散时间傅里叶变换DTFT的结果,归一角频率ω=2π f/FS =2π k/N;
步骤(6):依据频偏参数γ计算被测信号频率f γ ,计算公式为:
。
进一步地,所述的步骤(3)搜索最高和次高谱线的处理为:在设定频率范围所对应的离散频率序号范围[ks ,
ke ]内,搜索|Y(k)|最大的谱线作为最高谱线,其离散频率序号为k 1,同时选择该谱线两侧、离散频率序号为k 1-1和k 1+1的相邻两条谱线。
进一步地,所述的步骤(5)采用中间参数实测值λ与频偏参数γ关系函数的反函数公式计算频偏参数γ,该反函数为:
。
进一步地,所述的步骤(5)采用中间参数实测值λ与频偏参数γ关系函数的反函数的逼近多项式公式计算频偏参数γ,该逼近多项式公式为:
,
其中,M是逼近多项式的最高次数,am (m=0:M)是多项式第m次项λ m 的系数。
本发明信号频率测量方法的设计原理是:假设一个频率为f 0、幅值为、初相位为的单一频率信号x(t),在经过了采样率为Fs的模数变换后得到如下形式的离散信号:
。
如果所加窗函数的时域形式为w(n),其离散时间傅里叶变换DTFT得到的连续频谱为W(ω),则忽略负频点-f 0处频峰的旁瓣影响,在正频点f 0附近的连续频谱函数可以表达为:
,
上式进行离散抽样,即可得到离散傅立叶变换DFT的表达式为:
,
其中,离散频率间隔为Δf=FS /N。
余弦窗函数是DFT最为常用的一类窗函数。对应余弦窗函数的统一时域形式为:
,
余弦窗w(n)的离散时间傅里叶变换DTFT结果为:
,
其中:
,
在信号DTFT频谱曲线的主瓣内,且当N较大时,近似有:
,
当时,上式取等号。依据常用窗函数系数,在主瓣-H<k<H内,其相邻两条谱线和的相位相差近似为π;而对应H<k<N/2的旁瓣内和接近同相位。由此,当谱线Y(k 1)和Y(k 2)的相位相差0或π时,有:
,
由此,不必计算谱线幅值,直接通过谱线实部和虚部的绝对值之和,即可实现最高谱线的搜索、以及中间参数λ的计算。
进一步地,通过对多数窗函数进一步进行频域处理所得到的新的窗函数,将进一步抑制其旁瓣,由此可以减小被测信号DFT之后负频点、谐波频点对待测信号强度的旁瓣干扰。由此,本发明以(Y(k 1)-Y(k 1-1))=Y 2-Y 1和(Y(k 1)-Y(k 1+1))=
Y 2-Y 3的比值为基础,识别实际信号频点位置。同时,为了减少开方运算,采用谱线的实部和虚部分别进行计算,即:
。
本发明方法基于上述原理设计,减少了乘法和开方运算,简化了信号频率测量算法实现。计算公式中采用相邻谱线取差的形式,进一步抑制其他信号旁瓣对测量的影响,提高了测量准确度。
附图说明
图1是本发明的基于三条DFT复数谱线的信号频率测量方法的计算流程图。
具体实施方式
下面结合附图1所示的计算流程图对本发明的两个具体实施方式作进一步的说明。这两个实施例应用于对50Hz附近频率信号进行测量。第一个具体实施方式采用哈宁(Hanning)窗,其具体步骤如下:
步骤(1):将采样率Fs=1500Hz、连续截取N=512点的信号x(n),进行加窗处理得到加窗信号y(n),加窗处理公式为:
,
其中w(n)选择N=512点的Hanning窗函数序列,即:
,n = 0:(N-1);
步骤(2):对加窗信号y(n)进行离散傅里叶DFT变换,得到离散频谱Y(k),其中离散频率序号k = 0:(N-1);
步骤(3):在离散频率序号范围[15,
23],即对应频率[43.945, 67.383]Hz的范围内,搜索最大的谱线作为最高谱线,其离散频率序号为k 1,同时选择该谱线两侧、离散频率序号为k 1-1和k 1+1的相邻两条谱线;
步骤(4):依据k 1-1、k 1和k 1+1对应的三条谱线Y 1 =Y(k 1-1)、Y 2 =Y(k 1)和Y 3 =Y(k 1+1)计算中间参数λ:
;
步骤(5):求解如下方程中的频偏参数γ:
;
由此,计算频偏参数γ的函数公式为:
,
步骤(6):依据频偏参数γ计算被测信号频率f γ ,计算公式为:
。
第二个具体实施方式采用布莱克曼(BlackMan)窗,其具体步骤如下:
步骤(1):将采样率Fs=1500Hz、连续截取 N=512点的信号x(n),进行加窗处理得到加窗信号y(n),加窗处理公式为:
,
其中w(n)选择N=512点的布莱克曼(BlackMan)窗函数序列,即:
,n = 0:(N-1);
步骤(2):对加窗信号y(n)进行离散傅里叶DFT变换,得到离散频谱Y(k),其中离散频率序号k = 0:(N-1);
步骤(3):在离散频率序号范围[15,
23],即对应频率[43.945, 67.383]Hz的范围内,搜索最大的谱线作为最高谱线,其离散频率序号为k 1,同时选择该谱线两侧、离散频率序号为k 1-1和k 1+1的相邻两条谱线;
步骤(4):依据k 1-1、k 1和k 1+1对应的三条谱线Y 1 =Y(k 1-1)、Y 2 =Y(k 1)和Y 3 =Y(k 1+1)计算中间参数λ:
;
步骤(5):采用中间参数实测值λ与频偏参数γ关系函数的反函数的逼近多项式公式计算频偏参数γ,采用最高为7次的逼近多项式公式,形式如下:
,
步骤(6):依据频偏参数γ计算被测信号频率f γ ,计算公式为:
。
依据第一个和第二个实施方式,分别输入相同的一组仿真测试数据,以验证两个实施例的计算结果。该输入信号x(n)是基波频率f 1为50.1
Hz、包含2至9次谐波的信号,具体形式为:
其中,基波和各次谐波的幅值分别是:1, 0.02, 0.1,
0.01, 0.05, 0.0, 0.02, 0.0, 0.01;初始相位分别是-23.1°, 115.6°, 59.3°, 52.4°, 123.8°, 161.8°, -31.8°, 119.9°, -63.7°。
第一个实施方式中,离散频率序号范围15至23的9条谱线计算结果依次是:0.152783905+j1.76229599,
4.70210906+j54.2254920, -10.9858392-j126.688437, 6.36734959+j73.4295187,
-0.221526634-j2.55345712, -0.0512202109-j0.589197696,
-0.0199885230-j0.228698046, -0.00996303987-j0.112669252。不论幅值比较还是采用实部虚部绝对值之和,均可得到最高谱线对应k 1=17,同时选择该谱线两侧、离散频率序号为k 1-1=16和k 1+1=18的相邻两条谱线。于是,中间参数计算结果为λ=0.05039996845,频偏参数γ=2∙λ=
0.1007999369,最终测量得到的信号频率为50.099999815Hz,测量误差-0.000000369%。
第二个实施方式中,离散频率序号范围15至23的9条谱线计算结果依次是:-0.63151373-j7.27542732,
4.87561219+j56.23242377, -9.38422261-j108.21320679, 6.10560557+j70.41558734,
-1.11464837-j12.84931264, 0.00734999+j0.0889938, -0.00022265+j0.00101658,
0.00048733+j0.00853355。不论幅值比较还是采用实部虚部绝对值之和,均可得到最高谱线对应k 1=17,同时选择该谱线两侧、离散频率序号为k 1-1=16和k 1+1=18的相邻两条谱线。于是,中间参数计算结果为λ=0.0413416339。带入逼近多项式后,得到频偏参数γ=0.1007992629,最终测量得到的信号频率为50.099997840Hz,测量误差-0.000004311%。
Claims (4)
1.一种基于三条DFT复数谱线的信号频率测量方法,其特征在于该信号频率测量方法包含如下步骤:
步骤(1):将采样率为FS 、采样点为连续截取的采样信号x(n),进行加窗处理得到加窗信号y(n),加窗处理公式为:
,
其中为N点的窗函数序列,n = 0:(N-1);
步骤(2):对加窗信号y(n)进行离散傅里叶DFT变换,得到离散频谱Y(k),其中离散频率序号k = 0:(N-1);
步骤(3):在设定频率范围所对应的离散频率序号范围[ks ,
ke ]内,搜索最大的谱线作为最高谱线,其离散频率序号为k 1,同时选择该谱线两侧、离散频率序号为k 1-1和k 1+1的相邻两条谱线,其中Re(Z)表示取复数Z的实部值,Im(Z)表示取复数Z的虚部值;
步骤(4):依据k 1-1、k 1和k 1+1对应的三条谱线Y 1 =Y(k 1-1)、Y 2 =Y(k 1)和Y 3 =Y(k 1+1)计算中间参数λ:
步骤(5):求解如下方程中的频偏参数γ:
,
其中,W(ω)是窗函数w(n)的离散时间傅里叶变换DTFT的结果,归一角频率ω=2πf/FS =2πk/N;
步骤(6):依据频偏参数γ计算被测信号频率f γ ,计算公式为:
。
2.根据权利要求1所述的基于三条DFT复数谱线的信号频率测量方法,其特征在于:所述的步骤(3)搜索最高和次高谱线的处理为:在设定频率范围所对应的离散频率序号范围[ks ,
ke ]内,搜索|Y(k)|最大的谱线作为最高谱线,其离散频率序号为k 1,同时选择该谱线两侧、离散频率序号为k 1-1和k 1+1的相邻两条谱线。
3.根据权利要求1所述的基于三条DFT复数谱线的信号频率测量方法,其特征在于:所述的步骤(5)采用中间参数实测值λ与频偏参数γ关系函数的反函数公式计算频偏参数γ,该反函数为:
。
4.根据权利要求1所述的基于三条DFT复数谱线的信号频率测量方法,其特征在于:所述的步骤(5)采用中间参数实测值λ与频偏参数γ关系函数的反函数的逼近多项式公式计算频偏参数γ,该逼近多项式公式为:
,
其中,M是逼近多项式的最高次数,am (m=0:M)是多项式第m次项λ m 的系数。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410435495.0A CN105372492B (zh) | 2014-08-31 | 2014-08-31 | 基于三条dft复数谱线的信号频率测量方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410435495.0A CN105372492B (zh) | 2014-08-31 | 2014-08-31 | 基于三条dft复数谱线的信号频率测量方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105372492A true CN105372492A (zh) | 2016-03-02 |
CN105372492B CN105372492B (zh) | 2018-06-22 |
Family
ID=55374862
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201410435495.0A Expired - Fee Related CN105372492B (zh) | 2014-08-31 | 2014-08-31 | 基于三条dft复数谱线的信号频率测量方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN105372492B (zh) |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105822289A (zh) * | 2016-03-25 | 2016-08-03 | 重庆科技学院 | 用于油井动液面检测的频率估算方法 |
CN106018956A (zh) * | 2016-08-10 | 2016-10-12 | 北京妙微科技有限公司 | 一种加窗谱线插值的电力系统频率计算方法 |
CN107064630A (zh) * | 2017-03-31 | 2017-08-18 | 许继集团有限公司 | 一种电力系统频率测量方法及装置 |
WO2017198111A1 (zh) * | 2016-05-16 | 2017-11-23 | 深圳市中兴微电子技术有限公司 | 一种频偏估计方法、装置及计算机存储介质 |
Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6329806B1 (en) * | 1994-06-01 | 2001-12-11 | Siemens Ag | Process for determining the harmonic oscillations of the fundamental component of an electrical signal |
FR2897441A1 (fr) * | 2006-02-10 | 2007-08-17 | Thales Sa | Recepteur numerique large bande de mesure de frequence |
CN101179549A (zh) * | 2007-12-14 | 2008-05-14 | 清华大学 | 采用三点加权插值算法的通信信号载频估计方法 |
CN101900761A (zh) * | 2009-11-05 | 2010-12-01 | 中国航天科技集团公司第五研究院第五一四研究所 | 一种高准确度非整周期采样谐波分析测量方法 |
CN102435844A (zh) * | 2011-11-01 | 2012-05-02 | 南京磐能电力科技股份有限公司 | 一种频率无关的正弦信号相量计算方法 |
CN102435845A (zh) * | 2011-10-19 | 2012-05-02 | 广西电网公司电力科学研究院 | 一种Nuttall自卷积窗加权傅里叶变换的谐波电能计量方法 |
CN102832620A (zh) * | 2012-08-31 | 2012-12-19 | 天津理工大学 | 基于加窗全相位fft的apf谐波检测与控制系统及方法 |
CN103576002A (zh) * | 2013-11-11 | 2014-02-12 | 华北电力大学(保定) | 一种容性绝缘设备介质损耗角的计算方法 |
CN103941090A (zh) * | 2014-04-22 | 2014-07-23 | 国家电网公司 | 基于谱线能量插值的谐波测量方法 |
-
2014
- 2014-08-31 CN CN201410435495.0A patent/CN105372492B/zh not_active Expired - Fee Related
Patent Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6329806B1 (en) * | 1994-06-01 | 2001-12-11 | Siemens Ag | Process for determining the harmonic oscillations of the fundamental component of an electrical signal |
FR2897441A1 (fr) * | 2006-02-10 | 2007-08-17 | Thales Sa | Recepteur numerique large bande de mesure de frequence |
CN101179549A (zh) * | 2007-12-14 | 2008-05-14 | 清华大学 | 采用三点加权插值算法的通信信号载频估计方法 |
CN101900761A (zh) * | 2009-11-05 | 2010-12-01 | 中国航天科技集团公司第五研究院第五一四研究所 | 一种高准确度非整周期采样谐波分析测量方法 |
CN102435845A (zh) * | 2011-10-19 | 2012-05-02 | 广西电网公司电力科学研究院 | 一种Nuttall自卷积窗加权傅里叶变换的谐波电能计量方法 |
CN102435844A (zh) * | 2011-11-01 | 2012-05-02 | 南京磐能电力科技股份有限公司 | 一种频率无关的正弦信号相量计算方法 |
CN102832620A (zh) * | 2012-08-31 | 2012-12-19 | 天津理工大学 | 基于加窗全相位fft的apf谐波检测与控制系统及方法 |
CN103576002A (zh) * | 2013-11-11 | 2014-02-12 | 华北电力大学(保定) | 一种容性绝缘设备介质损耗角的计算方法 |
CN103941090A (zh) * | 2014-04-22 | 2014-07-23 | 国家电网公司 | 基于谱线能量插值的谐波测量方法 |
Non-Patent Citations (1)
Title |
---|
王涛 等: "频差对方向图综合的影响及解决方法", 《通信技术》 * |
Cited By (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105822289A (zh) * | 2016-03-25 | 2016-08-03 | 重庆科技学院 | 用于油井动液面检测的频率估算方法 |
CN105822289B (zh) * | 2016-03-25 | 2019-03-19 | 重庆科技学院 | 用于油井动液面检测的频率估算方法 |
WO2017198111A1 (zh) * | 2016-05-16 | 2017-11-23 | 深圳市中兴微电子技术有限公司 | 一种频偏估计方法、装置及计算机存储介质 |
CN107395287A (zh) * | 2016-05-16 | 2017-11-24 | 深圳市中兴微电子技术有限公司 | 一种频偏估计方法和装置 |
CN107395287B (zh) * | 2016-05-16 | 2019-12-13 | 深圳市中兴微电子技术有限公司 | 一种频偏估计方法和装置 |
CN106018956A (zh) * | 2016-08-10 | 2016-10-12 | 北京妙微科技有限公司 | 一种加窗谱线插值的电力系统频率计算方法 |
CN106018956B (zh) * | 2016-08-10 | 2018-10-16 | 北京妙微科技有限公司 | 一种加窗谱线插值的电力系统频率计算方法 |
CN107064630A (zh) * | 2017-03-31 | 2017-08-18 | 许继集团有限公司 | 一种电力系统频率测量方法及装置 |
CN107064630B (zh) * | 2017-03-31 | 2019-11-12 | 许继集团有限公司 | 一种电力系统频率测量方法及装置 |
Also Published As
Publication number | Publication date |
---|---|
CN105372492B (zh) | 2018-06-22 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN102435844B (zh) | 一种频率无关的正弦信号相量计算方法 | |
CN103308804B (zh) | 基于快速k-s变换电能质量扰动信号时频参数提取方法 | |
CN104897960B (zh) | 基于加窗四谱线插值fft的谐波快速分析方法及系统 | |
CN107247182B (zh) | 一种基于量测相量数据的间谐波分量还原方法 | |
CN103399203B (zh) | 一种基于复合迭代算法的谐波参数高精度估计方法 | |
CN109030941A (zh) | Hanning自乘卷积窗FFT三谱线插值谐波分析方法 | |
CN108414833B (zh) | 一种信号分量频率的精确估计方法 | |
CN103308766A (zh) | 一种基于凯撒自卷积窗双谱线插值fft谐波分析方法及其装置 | |
CN104597321B (zh) | 基于四条离散傅里叶复数谱线的信号频率测量方法及装置 | |
CN105137180B (zh) | 基于六项余弦窗四谱线插值的高精度谐波分析方法 | |
CN103353550A (zh) | 一种测量电力系统信号频率及谐波参数的方法 | |
CN102818930B (zh) | 一种高精度快速计算电力谐波参数的方法 | |
CN106018956B (zh) | 一种加窗谱线插值的电力系统频率计算方法 | |
CN107271774B (zh) | 一种基于频谱泄漏校正算法的apf谐波检测方法 | |
CN104391178A (zh) | 一种基于Nuttall窗的时移相位差稳态谐波信号校正方法 | |
CN110837001A (zh) | 一种电力系统中谐波和间谐波的分析方法与装置 | |
CN105372492A (zh) | 基于三条dft复数谱线的信号频率测量方法 | |
CN101216512A (zh) | 一种非正弦周期信号实时高精度检测方法 | |
CN109541312A (zh) | 一种新能源汇集地区次同步谐波检测方法 | |
CN104297740A (zh) | 基于相位分析的雷达目标多普勒谱估计方法 | |
CN105486921A (zh) | 凯撒三阶互卷积窗三谱线插值的谐波与间谐波检测方法 | |
CN103091545A (zh) | 一种频率无关的正弦信号相量半波计算方法 | |
CN103941090A (zh) | 基于谱线能量插值的谐波测量方法 | |
CN106546949A (zh) | 一种基于频率预估计的双阵元正弦信号来波方向估计方法 | |
Li et al. | Frequency estimation based on modulation FFT and MUSIC algorithm |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
TA01 | Transfer of patent application right |
Effective date of registration: 20180419 Address after: 213161 A1002, 3, Chuang Kong port, Changzhou science and Education Town, 18 Wujin Road, Wujin, Changzhou. Applicant after: CHANGZHOU HAOYUN INDUSTRIAL CONTROL TECHNOLOGY Co.,Ltd. Address before: 100094, new building, building 4, 7 Hui Feng Road, Beijing, Haidian District, 407 Applicant before: Shengji Upgo (Beijing) Technology Co.,Ltd. |
|
TA01 | Transfer of patent application right | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20180622 |
|
CF01 | Termination of patent right due to non-payment of annual fee |