TWI489298B - 一種分析多尺度之固有熵値的方法 - Google Patents

一種分析多尺度之固有熵値的方法 Download PDF

Info

Publication number
TWI489298B
TWI489298B TW103100339A TW103100339A TWI489298B TW I489298 B TWI489298 B TW I489298B TW 103100339 A TW103100339 A TW 103100339A TW 103100339 A TW103100339 A TW 103100339A TW I489298 B TWI489298 B TW I489298B
Authority
TW
Taiwan
Prior art keywords
intrinsic
entropy
time series
time
nth
Prior art date
Application number
TW103100339A
Other languages
English (en)
Other versions
TW201528003A (zh
Inventor
葉家榮
黃鍔
羅孟宗
彭仲康
Original Assignee
國立中央大學
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 國立中央大學 filed Critical 國立中央大學
Priority to TW103100339A priority Critical patent/TWI489298B/zh
Priority to US14/245,431 priority patent/US20150193376A1/en
Application granted granted Critical
Publication of TWI489298B publication Critical patent/TWI489298B/zh
Publication of TW201528003A publication Critical patent/TW201528003A/zh

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Algebra (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Databases & Information Systems (AREA)
  • Information Retrieval, Db Structures And Fs Structures Therefor (AREA)
  • Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)

Description

一種分析多尺度之固有熵值的方法
本發明係有關於一種時間序列訊號之分析方法，以多個固有時間尺度的亂度特徵，來表現一個動態系統在不同時間尺度的動態特性。
熵值是一個描述動態系統複雜度的物理量。在控制論、機率論、數論、天體物理、生命科學等的重要領域，經常利用熵值所表現之特徵值或變化量代表該領域所定義之系統的動態特性，雖然傳統熵值的評估方法，如近似熵值和取樣熵值，能夠以一個熵值表現複雜系統之整體亂度特性，然而熵值並非只是複雜系統之宏觀特徵，基於一個複雜系統涵蓋多個不同時間尺度的調控機制，其對應於不同時間尺度所表現之動態特性亦會有差異；由於現有之評估方法並無法由單一熵值判讀出複雜系統不同尺度的動態特性，故相關之技術人員亦積極發展適切之解析方法。
而目前之相關領域的研究，如美國專利案US61/195,894以及US12/411,539所述，係有利用經驗模態分解法(Empirical Mode Decomposition,EMD)與希爾伯特-黃頻譜分析(Hilbert Spectral Analysis,HSA)的組合，以其方法之分析轉換提供一時間-頻率分析；其中，經驗模態分解法可 將一時間序列訊號中，分解出複數個固有模態函數，其實施之過程如下所述：(1)找出時間序列訊號中之所有局部極值，並利用三次曲線(cubic spline)連接所有的局部最大值(最小值)作為上包絡線(下包絡線)；(2)將數據與兩包絡線的平均包絡線相減，以取得第一分量；並(3)重複步驟1和步驟2，直到上下包絡線對稱於時間軸。
簡而言之，經驗模態分解法是一種自適性的模態分解法，按 照時間序列訊號的動態本質由高頻的擾動到低頻的擾動逐步拆解出該時間序列訊號的複數個固有模態函數。
本發明進一步為利用經驗模態分解法以及創新的數值分析方 法可量化出一時間序列訊號於不同時間尺度之熵值指標，則相關領域之技術人員係可利用本發明所提供之方法，由複雜系統所表現之特定時間尺度之熵值差異，作為評估系統狀態之細部動態特徵的參考數值。
有鑑於此，本發明提供一種分析一系統之固有熵值的方法，用以解析非線性及非穩態之時間序列訊號的動態特性表現，該方法之特徵在於具有下列步驟：
步驟A. 取得該系統之時間序列訊號，該系統係為一非線性及非穩態之時間序列(nonlinear and non-stationary time series)的動態系統，並利用一非線性及非穩態之模態分解方法將該時間序列訊號拆解為複數個固有模態函數(Intrinsic Mode Function,IMF)，該非線性及非穩態之模態分解方法係為一經驗模態分解法。
該些固有模態函數係為一時間序列函數，以該些固有模態函數的平均週期代表固有模態函數之固有時間尺度，且該些固有模態函數之平均週期由小至大依序為一第1固有時間尺度、一第2固有時間尺度、...、一第n-1固有時間尺度至一第n固有時間尺度，則該些固有模態函數之組合，係可表現該系統之時間序列訊號為一結合不同時間尺度的非線性以及非穩態之時間序列。
步驟B. 選取一第1時間序列，該第1時間序列為該第1固有時間尺度之固有模態函數，並以複數個粗粒化尺度粗粒化該第1時間序列，該些粗粒化尺度係為該系統之複雜度分佈，得到一第1粗粒化時間序列組。
步驟C. 於步驟B之後，利用一熵值分析方法計算該第1粗粒化時間序列組，該熵值分析方法係為一樣本熵演算法(sample entropy)，得到複數個第1粗粒化時間序列組之熵值，並取該些熵值之最大值為一第1固有時間尺度之固有熵值。
步驟D. 於步驟C之後，選取一第n時間序列，該第n時間序列為該第1固有時間尺度到該第n固有時間尺度之固有模態函數的組成，並取該第n時間序列之標準差代入該熵值分析方法，得到一第n粗粒化時間序列組以及複數個第n粗粒化時間序列組之熵值。
步驟E. 於步驟D之後，將該些第n-1粗粒化時間序列組以及該第n粗粒化時間序列組之熵值相減，得到複數個粗粒化尺度之熵值差，並取該些熵值差之最大者為一第n固有時間尺度之固有熵值。
步驟F. 重覆步驟D至步驟E所敘之內容，選取複數個時間序列，直至得到一第2固有時間尺度、一第3固有時間尺度、...、一第n-1固有 時間尺度以及一第n固有時間尺度之固有熵值。
步驟G. 產出一固有熵值組或一固有熵值特徵之圖式，該固有熵值組係包含該系統之時間序列訊號的該些固有時間尺度與該些固有時間尺度之固有熵值，進一步將該固有熵值組或該固有熵值特徵之圖式與一資料庫進行比對。
關於本發明之優點與精神，以及更詳細的實施方式可以藉由以下的實施方式以及所附圖式得到進一步的瞭解。
S100至S114‧‧‧分析固有熵值的方法之步驟
S0‧‧‧原始時間尺度
IMF1、IMF2至IMFn‧‧‧固有時間尺度
S1‧‧‧第1時間序列
F1、F2至Fn‧‧‧粗粒化尺度
S1F‧‧‧第1粗粒化序列組
S1F1、S1F2至S1Fn‧‧‧第1粗粒化序列
S2‧‧‧第2時間序列
第一圖顯示本發明一實施例中分析固有熵值之方法流程圖；第二圖顯示本發明一實施例中拆解時間序列訊號為複數個時間尺度固有模態函數之方法示意圖；第三圖顯示本發明一實施例中取得第1時間序列的粗粒化時間序列組之方法示意圖；第四圖顯示本發明一實施例中組成第2時間序列之方法示意圖；第五圖顯示人體心跳訊號於固有模態函數之時間尺度以及粗粒化尺度之解析；以及第六圖顯示人體心跳訊號於特定之固有熵值以及其數值之解析。
為讓 鈞局貴審查委員及習於此技術人士，對本發明之功效完全了解，茲配合圖示及圖號，就本發明較佳之實施例說明如下：
請參照第一圖所示，其係本發明之分析固有亂度值的方法之實施步驟，本發明所提之方法，用以解析複雜系統之動態特性表現，該方法之特徵在於具有下列步驟：
先如步驟S100所示，取得一系統之一時間序列訊號，該時間序列訊號係具有一原始時間尺度；最佳的，該系統係為一非穩態以及非線性之時間序列的動態系統；在一實施例中，本發明係可取得一生理系統、一工程系統或一環境系統等之時間序列訊號，本發明不以此為限。
接著，如步驟S102所示，利用一非線性及非穩態之模態分解方法將該時間序列訊號拆解為複數個固有模態函數，其中各別固有模態函數具有一平均週期，該平均週期由小至大依序為一第1固有時間尺度、一第2固有時間尺度、...、一第n-1固有時間尺度至一第n固有時間尺度；最佳的，該非線性及非穩態之模態分解方法係為一經驗模態分解法。
在一實施例中，請參照第二圖所示，係取得一人體生理系統之一時間訊號，該時間訊號為是一600秒內的個體心跳數值，其600秒係為該時間訊號之原始時間尺度S0；接著，利用經驗模態分解法拆解為複數個固有模態函數，以10秒為該些固有模態函數之平均週期的最小值，其10秒係為第1固有時間尺度IMF1；次以30秒為該些固有模態函數之平均週期的第二最小值，其30秒係為第2固有時間尺度IMF2；而後以600秒為該些固有模態函數之平均週期的最大值，其600秒係為第n固有時間尺度IMFn，則固有 模態函數之數量，本發明並不以此為限。
接著，如步驟S104所示，選取一第1時間序列，該第1時間序列為一第1固有時間尺度之固有模態函數，並以複數個粗粒化尺度粗粒化該第1時間序列，得到一第1粗粒化時間序列組；在一實施例中，該些粗粒化尺度係為該時間訊號的組合頻率。
接著，如步驟S106所示，利用一熵值分析方法計算該第1粗粒化時間序列組，得到複數個第1粗粒化時間序列組之熵值，並取該些熵值之最大值為一第1固有時間尺度之固有熵值；最佳的，其中該熵值分析方法係為一樣本熵演算法。
在一實施例中，請參照第三圖所示，係選取以第1固有時間尺度IMF1之固有模態函數為第1時間序列S1，並以複數個粗粒化尺度(F1、F2、...、Fn)粗粒化該第1時間序列S1，得到一第1粗粒化時間序列組S1F；接著，經由樣本熵演算法得到複數個第1粗粒化時間序列組(S1F1、S1F2、...、S1Fn)之熵值，若其中S1F2之熵值為該第1粗粒化時間序列組S1F之最大值，則S1F2序列之熵值為第1固有時間尺度IMF1之固有熵值，則粗粒化尺度之數量，本發明不並以此為限。
接著，如步驟S108所示，選取一第n時間序列，該第n時間序列為該第1固有時間尺度到該第n固有時間尺度之固有模態函數的組成，並取該第n時間序列之標準差代入該熵值分析方法，得到一第n粗粒化時間序列組以及複數個第n粗粒化時間序列組之熵值。
接著，如步驟S110所示，將該些第n-1粗粒化時間序列組以及該第n粗粒化時間序列組之熵值相減，得到複數個粗粒化尺度之熵值差， 並取該些熵值差之最大者為一第n固有時間尺度之固有熵值。
在一實施例中，請參照第四圖所示，係選取以第1固有時間 尺度IMF1到該第2固有時間尺度IMF2之固有模態函數的組成為第2時間序列S2，並取該第2時間序列S2之標準差代入樣本熵演算法，得到一第2粗粒化時間序列組S2F以及複數個第2粗粒化時間序列組(S2F1、S2F2、...、S2Fn)之熵值，則固有模態函數組成之數量以及粗粒化尺度之數量，本發明不並以此為限。
進一步將該些第2粗粒化時間序列組(S2F1、S2F2、...、S2Fn) 以及該第1粗粒化時間序列組(S1F1、S1F2、...、S1Fn)之熵值相減，得到複數個粗粒化尺度(F1、F2、...、Fn)之熵值差(S2D1、S2D2、...、S2Dn)，若其中粗粒化尺度F5具有最大熵值差，則該熵值S2D5為一第2固有時間尺度IMF2之固有熵值。
在一實施例中，係選取以第1固有時間尺度IMF1到該第3固 有時間尺度IMF3之固有模態函數的組成為第3時間序列S3，並取該第3時間序列S3之標準差代入樣本熵演算法，得到一第3粗粒化時間序列組S3F以及複數個第n粗粒化時間序列組(S3F1、S3F2、...、S3Fn)之熵值。
進一步將該些第3粗粒化時間序列組(S3F1、S3F2、...、S3Fn) 以及該第2粗粒化時間序列組(S2F1、S2F2、...、S2Fn)之熵值相減，得到複數個粗粒化尺度(F1、F2、...、Fn)之熵值差(S3D1、S3D2、...、S3Dn)，若其中粗粒化尺度F7具有最大熵值差，則該熵值S3D7為一第3固有時間尺度IMF3之固有熵值。
接著，如步驟S112所示，重覆選取複數個時間序列，直至得 到該第2固有時間尺度IMF2、該第3固有時間尺度IMF3、...、該第n-1固有時間尺度IMFn-1以及該第n固有時間尺度IMFn之固有熵值。
而後，如步驟S114所示，產出一固有熵值組，係包含該時間 序列訊號之該些固有時間尺度與該些固有時間尺度之固有熵值，將該固有熵值組與一資料庫之對比數值組進行一比對，最佳地，亦產出一固有熵值特徵之圖式。
在一實施例中，請參照第五圖與第六圖所示，係揭露本發明 之分析固有亂熵值的方法所解析之人體心跳訊號的實施態樣，統計141件心跳時間訊號的複數固有熵值之組合，其中44件係為36.39±9.4歲的健康年輕個體之案例，28件係為66.2±3.7歲的健康年長個體之案例，22件係以輕度充血性心臟衰竭(Congestive Heart Failure,CHF)個體為CHF 1-2之案例，22件係以重度充血性心臟衰竭為CHF 3-4之案例，另有25件係取得自心房纖維顫動(Atrial Fibrillation,AF)個體之案例，由於人體心跳訊號對應於不同生理狀況以及病理特徵會表現出不同的特徵值，如第五圖a至e所示之圖形，係以x軸表示固有模態函數之時間尺度與y軸表示固有模態函數之粗粒化尺度，為說明本發明之方法解析該些個體案例之心跳時間訊號，所取得之固有熵值特徵之圖式，比較心臟疾病患者與正常人之固有熵值特徵之圖式，係可發現心臟疾病患者於各尺度數值之表現均有下降的趨勢，由其在大尺度時，該數值下降最為明顯，醫療人員可將訊號特徵與資料庫作比對，判斷出正確之疾病表現；又如第六圖所示之圖表，係以x軸表示特定之固有熵值與y軸表示個體案例之固有熵值數值，為說明因個體之生理狀況或病理特徵的不同，於固有熵值的表現上亦會有相異處，醫療人員可由固有熵值的表 現，建立準確之心臟疾病評估的依據。
上述實施例僅為說明本發明之原理及其功效，並非限制本發明。因此習於此技術之人士對上述實施例進行修改及變化仍不脫本發明之精神。本發明之權利範圍應如後述之申請專利範圍所列。
S100-S114‧‧‧分析固有熵值的方法之步驟

Claims (6)

  1. 一種分析一系統之固有熵值的方法，步驟為：A取得該系統之時間序列訊號，並利用一非線性及非穩態之模態分解方法將該時間序列訊號拆解為複數個固有模態函數，其中該些固有模態函數之平均週期由小至大依序為一第1固有時間尺度、一第2固有時間尺度直至一第n-1固有時間尺度和一第n固有時間尺度；B選取一第1時間序列，該第1時間序列為該第1固有時間尺度之固有模態函數，並以複數個粗粒化尺度粗粒化該第1時間序列，得到一第1粗粒化時間序列組；C於步驟B之後，利用一熵值分析方法計算該第1粗粒化時間序列組，得到複數個第1粗粒化時間序列組之熵值，並取該些熵值之最大值為一第1固有時間尺度之固有熵值；D於步驟C之後，選取一第n時間序列，該第n時間序列為該第1固有時間尺度到該第n固有時間尺度之固有模態函數的組成，並取該第n時間序列之標準差代入該熵值分析方法，得到一第n粗粒化時間序列組以及複數個第n粗粒化時間序列組之熵值；E於步驟D之後，將該些第n-1粗粒化時間序列組以及該第n粗粒化時間序列組之熵值相減，得到複數個粗粒化尺度之熵值差，並取該些熵值差之最大者為一第n固有時間尺度之固有熵值；F重覆步驟D至步驟E所敘之內容，選取複數個時間序列，直至得到一第2固有時間尺度、一第3固有時間尺度直至一第n-1固有時間尺度和一第n固有時間尺度之固有熵值； G於步驟F之後，產出一固有熵值組，係包含該時間序列訊號之該些固有時間尺度與該些固有時間尺度之固有熵值。
  2. 如申請範圍第1項所敘之分析方法，於步驟G更包含產出一固有熵值特徵之圖式。
  3. 如申請範圍第2項所敘之分析方法，於步驟G更包含將該固有熵值組或該固有熵值特徵之圖式與一資料庫進行比對。
  4. 如申請範圍第1項所敘之分析方法，其中該系統係為一非穩態以及非線性之時間序列的動態系統。
  5. 如申請範圍第1項所敘之分析方法，其中該非線性及非穩態之模態分解方法係為一經驗模態分解法。
  6. 如申請範圍第1項所敘之分析方法，其中該熵值分析方法係為一樣本熵演算法。
TW103100339A 2014-01-06 2014-01-06 一種分析多尺度之固有熵値的方法 TWI489298B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
TW103100339A TWI489298B (zh) 2014-01-06 2014-01-06 一種分析多尺度之固有熵値的方法
US14/245,431 US20150193376A1 (en) 2014-01-06 2014-04-04 Method of multi-scales intrinsic entropy analysis

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
TW103100339A TWI489298B (zh) 2014-01-06 2014-01-06 一種分析多尺度之固有熵値的方法

Publications (2)

Publication Number Publication Date
TWI489298B true TWI489298B (zh) 2015-06-21
TW201528003A TW201528003A (zh) 2015-07-16

Family

ID=53495299

Family Applications (1)

Application Number Title Priority Date Filing Date
TW103100339A TWI489298B (zh) 2014-01-06 2014-01-06 一種分析多尺度之固有熵値的方法

Country Status (2)

Country Link
US (1) US20150193376A1 (zh)
TW (1) TWI489298B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9916332B2 (en) * 2015-07-09 2018-03-13 Entit Software Llc Dataset chart scaling
CN111447872A (zh) * 2017-12-11 2020-07-24 艾德脑科技股份有限公司 生理信号分析装置与方法
CN108594161B (zh) * 2018-05-03 2020-06-19 国网重庆市电力公司电力科学研究院 一种电能表内异物声音信号降噪方法、系统
CN109580222A (zh) * 2018-12-04 2019-04-05 河北科技大学 基于变分模态分解-传递熵的轴承退化状态识别预测方法
CN110974259B (zh) * 2019-11-05 2022-06-24 华南师范大学 基于均值粗粒化的脑电情感识别方法、系统和存储介质
CN112515686B (zh) * 2020-11-30 2022-12-30 中国科学院空天信息创新研究院 一种脑电数据处理方法、装置以及计算机可读存储介质
CN113624710A (zh) * 2021-09-16 2021-11-09 首都师范大学 一种液体物质成分变化的太赫兹信息识别方法及其应用
CN114202921B (zh) * 2021-12-09 2023-08-04 洛阳师范学院 交通流量的多尺度符号动态熵分析方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090275853A1 (en) * 2008-04-30 2009-11-05 The General Electric Company Method and apparatus for monitoring physiological state of a subject
US7882167B2 (en) * 2005-02-18 2011-02-01 Beth Israel Deaconess Medical Center Complexity-based dynamical analysis of a network
CN102670190A (zh) * 2012-05-21 2012-09-19 南京大学 基于心率变异性非线性特性的充盈性心衰自动诊断方法
TW201300081A (zh) * 2011-06-17 2013-01-01 Ind Tech Res Inst 生理參數指標運算系統、方法、記錄媒體及電腦程式產品

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6738734B1 (en) * 1996-08-12 2004-05-18 The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration Empirical mode decomposition apparatus, method and article of manufacture for analyzing biological signals and performing curve fitting
WO2006089086A2 (en) * 2005-02-18 2006-08-24 Beth Israel Deaconess Medical Center, Inc. Complexity-based dynamical assay for assessing the toxicity and efficacy of pharmaceutical and other therapeutic interventions

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7882167B2 (en) * 2005-02-18 2011-02-01 Beth Israel Deaconess Medical Center Complexity-based dynamical analysis of a network
US20090275853A1 (en) * 2008-04-30 2009-11-05 The General Electric Company Method and apparatus for monitoring physiological state of a subject
TW201300081A (zh) * 2011-06-17 2013-01-01 Ind Tech Res Inst 生理參數指標運算系統、方法、記錄媒體及電腦程式產品
CN102670190A (zh) * 2012-05-21 2012-09-19 南京大学 基于心率变异性非线性特性的充盈性心衰自动诊断方法

Also Published As

Publication number Publication date
TW201528003A (zh) 2015-07-16
US20150193376A1 (en) 2015-07-09

Similar Documents

Publication Publication Date Title
TWI489298B (zh) 一種分析多尺度之固有熵値的方法
Henriques et al. Nonlinear methods most applied to heart-rate time series: a review
Azami et al. Amplitude-and fluctuation-based dispersion entropy
Thuraisingham et al. On multiscale entropy analysis for physiological data
Gómez-Herrero et al. Assessing coupling dynamics from an ensemble of time series
Costa et al. Multiscale analysis of heart rate dynamics: entropy and time irreversibility measures
Zhao et al. Determination of sample entropy and fuzzy measure entropy parameters for distinguishing congestive heart failure from normal sinus rhythm subjects
Amarantidis et al. Interpretation of entropy algorithms in the context of biomedical signal analysis and their application to EEG analysis in epilepsy
Sharma et al. Automated detection of hypertension using physiological signals: a review
Culbreth et al. Entropic approach to the detection of crucial events
Billeci et al. A combined independent source separation and quality index optimization method for fetal ECG extraction from abdominal maternal leads
Chien et al. Automatic sleep-arousal detection with single-lead EEG using stacking ensemble learning
Byeon et al. Pre-configured deep convolutional neural networks with various time-frequency representations for biometrics from ECG signals
Triwiyanto et al. Embedded machine learning using a multi-thread algorithm on a Raspberry Pi platform to improve prosthetic hand performance
Alqudah An enhanced method for real-time modelling of cardiac related biosignals using Gaussian mixtures
Zhao et al. Evaluation of a single-channel EEG-based sleep staging algorithm
Afghah et al. Game theoretic approach for systematic feature selection; application in false alarm detection in intensive care units
Jara et al. Using complexity–entropy planes to detect Parkinson’s disease from short segments of haemodynamic signals
Hua et al. A visual approach for the SARS (severe acute respiratory syndrome) outbreak data analysis
Chou et al. Effects of ecg data length on heart rate variability among young healthy adults
McLachlan et al. Pitch and plasticity: insights from the pitch matching of chords by musicians with absolute and relative pitch
Bajić et al. On quantization errors in approximate and sample entropy
Elgendi Eventogram: a visual representation of main events in biomedical signals
Roland Motion artifact suppression for insulated EMG to control myoelectric prostheses
Casson An analog circuit approximation of the discrete wavelet transform for ultra low power signal processing in wearable sensor nodes

Legal Events

Date Code Title Description
MM4A Annulment or lapse of patent due to non-payment of fees