JP2697957B2 - 時間により変化するスペクトル分析に関する方法および装置 - Google Patents

時間により変化するスペクトル分析に関する方法および装置

Info

Publication number
JP2697957B2
JP2697957B2 JP5516728A JP51672893A JP2697957B2 JP 2697957 B2 JP2697957 B2 JP 2697957B2 JP 5516728 A JP5516728 A JP 5516728A JP 51672893 A JP51672893 A JP 51672893A JP 2697957 B2 JP2697957 B2 JP 2697957B2
Authority
JP
Japan
Prior art keywords
spectrum
coefficient
signal
time
processor
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 - Lifetime
Application number
JP5516728A
Other languages
English (en)
Other versions
JPH07505949A (ja
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 ナショナル インストルメンツ コーポレイション
Publication of JPH07505949A publication Critical patent/JPH07505949A/ja
Application granted granted Critical
Publication of JP2697957B2 publication Critical patent/JP2697957B2/ja
Anticipated expiration legal-status Critical
Expired - Lifetime legal-status Critical Current

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01RMEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R23/00Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
    • G01R23/16Spectrum analysis; Fourier analysis

Landscapes

  • Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Complex Calculations (AREA)
  • Measurement Of Resistance Or Impedance (AREA)
  • Analysing Materials By The Use Of Radiation (AREA)
  • Radar Systems Or Details Thereof (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Description

【発明の詳細な説明】 限定された著作権放棄 この特許文書の開示の部分には、著作権保護の申し立
てをした素材が含まれている。著作権者は、特許の文書
や開示の複写物が、米国特許商標庁の出願書類や記録に
表れるように、だれにより行われようと異議を申し立て
ない。しかし、その他の権利はすべて留保している。
発明の背景 発明の分野 この発明は、デジタル信号処理に関するものであり、
特に、時間により変化する周波数成分を有する信号の特
性を定めるための方法および装置に関するものである。
関連する技術の説明 デジタル信号処理のための典型的なシステムは、入力
アナログ信号をデジタル・サンプル信号系列への変換、
デジタル・サンプル信号系列をデジタル信号に対して配
列やブロックレベルで計算を行うデジタル信号処理エン
ジンへの供給および計算結果に対する分析を含んでい
る。この分野は、音声処理、レーダーやソナーの画像認
識、水中音波探査器その他の技術に急速に拡大してい
る。
離散フーリェ変換(Discrete Fourier Transform DF
T)と高速フーリェ変換(Fast Fourier Transform FF
T)によるDFTの実現は、デジタル信号処理システムにお
ける基本的な処理技術である。これらの技術は、処理さ
れている信号に対して、時間により周波数成分がどのよ
うに変化しているのかを特徴づけることができないとい
う制限がある。
従来の時間により周波数成分が変化する信号を分析す
るシステムにおいては、所謂短時間フーリェ変換(STF
T)が用いられていた。このアルゴリズムは、入力信号
の多数の短いウィンドのフーリエ変換を計算することに
基づいている。各短い時間のウィンドごとに変換して結
合することにより、入力信号の時間により変化するスペ
クトルが求まる。しかし、STFT技術を用いて求めた時間
により変化するスペクトルは、よい分解能を有していな
い。また、STFT技術は、計算量が多い。
Gabor変換は別のデジタル信号処理の手法で、信号の
合成された時間−周波数表現を生成する(D.Gabor“The
ory of Communication"J.IEE(London),Vol.93 No.II
I,November,1946,pp.429−457を参照)。この応用は、
しかしながら、主として2重直交補助ウィンド関数γを
選択することが難しいために、制限されていた。最近、
有限サイクリック離散Gabor変換のための補助関数γを
設計するための枠組みが開発された(Wexler et al.“D
iscrete Gabor Expansions",Signal Processing,Vol.21
No.3,November,1990,pp.207−221.参照)。しかし、沢
山の信号処理の応用に対して、有限サイクリックGabor
変換は適当ではない。有限離散Gabor変換の一般的解法
が存在しないため、このデジタル処理技術の有用性が制
限されている。そして、時間的に変化するエネルギ・ス
ペクトルの作成に適用することに成功していない(Qna
n,et al.“Wigner Distribution Decomposition and Cr
oss−term Interference Cencellation"UMBC Technical
Report No.EER−91−1,University of Maryland,Janua
ry,1991を参照)。
発明の概要 本発明は、時間的に変化する周波数成分で特徴づけら
れた入力信号の時間的に変化するスペクトルを生成する
信号分析装置を提供するものである。この発明は、ま
た、時間で変化するスペクトルを基にして信号を分析す
る方法である。
したがって、第1の特徴として、本発明は、入力信号
を表すデジタル信号系列を生成する変換器を有する信号
分析装置を含む。デジタル信号系列は、直交的離散Gabo
r変換係数Cm,nをその信号系列に応答して計算し、この
係数に応答して入力信号エレルギの時間で変化するスペ
クトル(ここではGaborスペクトルと記述している)を
計算する第1のプロセッサに供給される。最後に、第2
のプロセッサはこのスペクトルをその後の分析や表示の
ために処理する。時間で変化するスペクトルが有用な一
つの特別な分析は、入力信号を別々の構成要素に分解す
ることである。別々の構成要素は色々な技術を用いて各
々分析される。
第1のプロセッサは入力バッファを有しており、バッ
ファ制御装置により、入力されるデジタル信号系列を受
信し、係数の計算のために適当なフォーマットで供給さ
れる。同時に、プロセッサは出力バッファを有してお
り、それには時間で変化するスペクトルを表すマトリッ
クスの値を、プロセッサ・エンジンにより、格納してい
る。
デジタル信号プロセッサは、ソフトウェアで制御され
た汎用コンピュータかハードウェアで実現されて、係数
と時間で変化するスペクトルを計算するようにプログラ
ムされている。そのルーティンは、次のステップを基に
している。
入力信号を表すデジタル信号系列を取り入れるステッ
プ、 デジタル信号系列をサンプリングして、長さLの複数
のウィンドウを定義し、複数のウィンドウの各々は、複
数のデジタル信号を含み、隣接のウィンドウに関連する
ΔMの長さのデジタル信号でシフトするステップ、 サンプリング間隔ΔNとΔMの直交的離散Gabor変換
係数Cm,nを、n=0からN−1まで各ウィンドウmに
対してと、複数のウィンドウm=0からM−1まで、Δ
NΔMがLより小さい値に対して、計算するステップ、 係数Cm,nに応答して、入力信号のエネルギの時間で
変化するスペクトルを計算するステップ、 入力信号を時間で変化するスペクトルに応答して分析
するステップ 直交的離散Gabor変換の係数は、非周期的な局所離散
ウィンドウ関数hおよびhと同様でhとγは重ならない
ような離散補助関数γを用いて計算される。時間で変化
するスペクトルは、クロス項を削除した(cross−term
deleted)Wigner−Ville分布を用いて計算される。
本発明の別の特徴は、離散ウィンドウ関数は、約0.5
(L/12)以下の分散σを有するガウス関数(Lは計
算で使用したウィンドウの長さ)である。
また本発明の別の特徴は、スペクトルが、係数を基に
したエネルギ分布関数で使用した予め計算した項目を有
するルックアップ・テーブルを使用して計算することで
ある。それぞれの係数に対して、スペクトル上の係数の
周囲の更新領域は、係数の大きさに対応して定められて
いる。エネルギ分布は、各係数の更新領域内で計算され
る。予め計算された項目のルックアップ・テーブルは、
最大の大きさの係数の更新領域をカバーするのに十分な
サイズを有するように計算される。この様に、定められ
た大きさのルックアップ・テーブルは、与えられたGabo
rスペクトルエンジンに対して使用することができる。
本発明の信号分析装置および方法は、音声処理、音や
電磁波の通信、音や電磁波の監視、レーダー、ソナー、
およびその他の信号処理の環境で遭遇する時間で変化す
る周波数成分で特徴づけられる信号に対して特に有用で
ある。例えば、複雑に信号のグループの中で、ある周波
数スペクトルを有する信号の始まりと終わりの時間を特
定する能力は、即ち他の複数の人が存在する部屋内で特
定の人がしゃべり始めるときに必要であるが、本発明の
Gaborスペクトル分布装置や方法により大幅に改良され
る。
本発明の他の特徴と効果は、明細書、図面、特許請求
の範囲をよく読むことで理解される。
図面の簡単な説明 第1図は、本発明の信号分析装置のブロック図であ
る。
第2A図および第2B図は、本発明のGaborスペクトルを
計算する方法をしめすフローチャートである。
第3図は、Gaborスペクトルの計算で使用されるデジ
タル信号のフレームの構成を示す図である。
第4図は、Gaborスペクトルの計算で使用されるサン
プルを行う分布を示す説明図である。
第5図は、本発明によるGaborスペクトルのエネルギ
分布に計算に使用する更新領域の例を示す図である。
第6図は、本発明のGaborスペクトルを生成する信号
分析装置のより詳細なブロック図である。
第7A図および第7B図は、それぞれ周波数ホッピング信
号および雑音のあるホッピング信号で、時間で変化する
周波数成分を有する入力信号を示した図である。
第8A図および第8B図は、それぞれ本発明を第7A図およ
び第7B図に示した入力信号に適用して生成されるGabor
スペクトルを示すグラフの図である。
第9図は、本発明の信号分析装置の応用に用いられる
局所化したウィンドウ関数hと離散補助関数γを示す図
である。
発明を実施するための最良の形態 本発明の実施例を図面を用いて詳細に説明する。
第1図は、本発明の信号分析装置を説明するブロック
図である。分析装置は、符号10で参照されるが、アナロ
グ/デジタル変換器11、デジタル信号プロセッサ12、デ
ジタル信号プロセッサに接続された制御部13、Gaborス
ペクトル分析装置14、そして表示装置15を有する。入力
信号s(t)が変換器や他の入力源から線16に供給され
る。アナログ/デジタル変換器11は入力信号s(t)を
表すデジタル値の系列s(i)を発生する。系列s
(i)はデジタル信号プロセッサ12に供給される。デジ
タル信号プロセッサは、デジタル信号処理(DSP)機能
を実行するようにプログラムされた汎用コンピュータを
使用することで、またはこの目的に適合した特別のデジ
タル信号処理(DSP)を行うエンジンで実施される。デ
ジタル信号プロセッサ12は離散Gabor変換(DGT)を、入
力信号s(i)のGabor係数を計算し、この係数に応答
しクロス項を削除した(cross−term deleted)Wigner
−Ville分布(PWDV)を基にして時間および周波数のエ
ネルギ分布を計算することにより実行する。エネルギ分
布GS(i,s)はデジタル信号プロセッサの出力から線17
上に供給される。エネルギ分布GS(i,s)はGaborペクト
ル分析装置14で分析され、表示装置15に供給するために
エネルギ分布の画像表示を生成する。また、Gaborスペ
クトル分析装置14は、他の分析機能、例えば、音声処
理、伝送、画像認識等で使用するために、時間で入力信
号の成分を分解すること等を行うことができる。スペク
トル分析装置14は、プログラムされた汎用コンピュータ
やこの目的のために設計された専用装置である。分析装
置14はデジタル信号プロセッサ12と同じプロセッサ上に
実現してもよい。
制御ブロック13は、1実施例では、付表Aのようなコ
ンピュータ・プログラムで構成されている。
付表Aに示されているプログラムは、LABVIEW(Natio
nal Instrumentsの商標)計測ソフトウェアシステムの
ために書かれたソースコード・ルーティンの一例であ
る。明らかに、広範囲の種類の他のソースコードで実装
することができる。
付表Aは、GaborSpectrum.c,GaborMagnitude.c、そし
てautoWV.cを含む3つのルーティンで構成されている。
GaborSpectrum.cはGabor係数の計算のためのデータを設
定する。ユーザは、スペクトルを計算するためのオーバ
ーサンプリング率,分散,およびウィンドウの長さのよ
うなパラメータを選択してルーティンを初期化する。パ
ラメータに関するより詳細な説明が下記に記載されてい
る。ルーティンGaborSpectrum.cはルーティンGaborMagn
itude.cを呼び出す。そのルーティンは、複素Gabor係数
を計算し、Gaborスペクトルを計算するのに使用する係
数の大きさを出力として供給する。ルーティンautoWV.c
は、係数に応答してエネルギ分布を計算し、Gaborスペ
クトルを出力する。
基本的には、Gaborスペクトル・アルゴリズムは、次
の1−8のステップで構成されている。
ステップ1: ウィンドウの長さLを選択する。ある実装では、値L
は常に自乗である。例えば、L=64,128,256,512,…等
である。Lを定めた後、周波数分解能が式(1)により
定められる。
ここでfsは変換器のサンプリング周波数である。
L=128またはL=256が音声処理に応用する場合に適
当である。
ステップ2: オーバーサンプリング率を選択する。オーバーサンプ
リング率は、式(2)により定められる。
ここで、ΔMは時間領域の増分、ΔNは周波数領域の
増分である。
(mΔM,nΔN)は、Gabor変換で見出されるガウス関
数の座標である。ここで、mは0から(SL/ΔM−1)
まで、SLはサンプル長、そして、nは0から(L/ΔN−
1)までである。図4はGabor変換アルゴリズムにおけ
る係数Cm,nの位置を示す図である。オーバーサンプリ
ング率として4が選択される。しかし、他の値も可能で
ある。オーバーサンプリング率が大きくなると、計算量
が増加する。オーバーサンプリング率の4が、出力スペ
クトルの分解能と必要な計算量との間の妥協として良い
ことが証明されている。
ステップ3: 式(3)となるように、ΔM,ΔNを選択する。
付表Aの実装では、L=128の場合、ΔM=8および
ΔN=4を選択した。しかし、他の選択も同様に適して
いる。ΔMとΔNの選択で、Gaborスペクトルの時間軸
上および周波数軸上のそれぞれの分解能を定まる。
ステップ4: L,ΔM,ΔNが与えられると、ウィンドウ関数の分散σ
は(4)の式で定められる。
これは、ウィンドウ関数hに対する局所化したガウス
関数を確保する。分散を大きくするとよりよい周波数領
域の分解能が得られる。分散を小さくするとよりよい時
間領域の分解能が得られる。
ステップ5: Lとσが与えられると、ガウス・ウィンドウ関数h
[L]が与えられる。そして、離散補助関数γ[L]
が、付表Bに示されているように、γが最小自乗誤差に
関してhと同様に定れられる。最適なシステムでは、γ
とhは示された領域の制限内で重ならない。補助関数γ
[L]とウィンドウ関数h[L]の例が第9図に示され
ている。
ステップ1−5は、準備的なステップである。実際の
スペクトル計算はステップ6以下で行われる。
ステップ6: 第3図は、どのようにして連続的な信号が離散化する
のか、そしてどのようにしてL点フレームの離散信号s
(i)の系列に適応するのかを示している。あるフレー
ムのL点が求められると、L点に対するGabor係数の大
きさを以下に従って計算する。
*は共役複素数を示す。
mΔMの値は入力データの使用され方を示す。次の例
は、この過程を説明するためである。
m=−L/ΔM+1,s(−L+ΔM),s(−L+Δ+
1),…,s(ΔM−1)が用いられる。
m=0,s(0),s(1)…,s(L−1)が用いられる; m=1,s(ΔM),s(ΔM+2),…,s(ΔM+L−
1)が用いられる; m=2,s(2ΔM),s(2ΔM+2),…,s(2ΔM+
L−1)が用いられる; … 一般的に、 s(mΔM),s(smΔM+1),…,s(mΔM+L−
1)が用いられる。
各々、係数Cm,nがn=0,1,・・・N−1の場合に対
して計算される。
示されているように、この計算は、FFTの順序の形式
を有している。
ステップ7: 係数Cm,nがs(i)に対して計算された後、信号の
スペクトルがクロス項を削除した(cross−term delete
d)Wigner−Ville分布を用いて計算される。Gabor数C
m,nは疎に位置しているので、スペクトルの完全な再構
築には、係数Cm,n間の点の補間が必要となる。Wigner
−Ville分布のクロス項の干渉(interference)を削除
後、信号のエネルギは以下の式で計算される。
ステップ8: (mΔM,nΔN)における各Gabor係数|Cm,n|に対し
て、スペクトル値GS8(i,k)がステップ7の式で計算さ
れ、前のGS8(i,k)と累算される。即ち、第5図に示さ
れているように、点(i,k)におけるエネルギは、係数C
22,C21,C31,C32そしてその他の分布の累算された寄与で
ある。計算機言語の表現を用いて、その累算は以下のよ
うに表される。
式(8)における定数4は正規化係数で、特定の応用
により変化するものである。
上記の項は、(mΔM,nΔN)におけるGabor係数から
離れた点(i,k)では、急速に減少するので、このステ
ップにおける反復は、点(mΔM,nΔN)の周囲の小さ
い領域に限定される。正確な領域の大きさは、応用に合
わせて調整される。また、上記指数の項は毎回計算する
必要はない。この項は予め計算され、ルックアップ・テ
ーブルに格納される。
スペクトル計算上の重要な問題は、いつ計算を停止す
るか即ち更新領域の大きさである。Cm,nの隣接する点
へのエネルギ寄与は、指数関数的に減少するので、(m
ΔM,nΔN)から遠く離れた点を計算する必要はない。
このため、各Gabor係数Cm,nに対して更新領域を定める
必要がある。更新領域は、第5図に示されるように、ス
ペクトルが計算される必要がある楕円として定義され
る。
i′=i−mΔM、k′=k−nΔNとすると、式
(9)は以下のように表される。
これは、第5図に60のような2次元の楕円である。
任意の係数に対する更新領域は、以下のように定義さ
れる。
ここで、|Cm,n|maxは、全てのGabor係数の最大の大き
さである。
したがって、各係数に対する更新領域は、以下ので表
される。
最大の更新領域は、|Cm,n|2=|Cm,n|max 2の場合であ
る。このとき、以下の式が導き出せる。
このように、i′maxとj′maxは、ルックアップ・テ
ーブルのサイズを与える。正のインデックスと値のみを
格納する必要がある。これは、式(B.2)はi′とk′
に対して偶関数であるからである。
一旦、テーブルが作成されると、式(8)計算はすぐ
できる。対応する項目をテーブル内で探し、その値と|C
m,n|2とを掛ける。計算上のコストが高い指数関数の計
算がこの様にして避けることができる。
第2A図と第2B図は、好ましい実施例に従ってGaborス
ペクトルを計算する基本的な手法が示されている。
プロセッサをまず初期化するルーティンが含まれる
(ブロック50)。この初期化は、上記に概略したステッ
プ1−5を実行することが含まれる。
次に、Gabor係数が第2A図のブロック51,52,53,54、そ
して55で実現されているアルゴリズムに従って計算され
る。
特に、インデックスiが−L+ΔMに設定され、イン
デックスmがL/ΔM+1に設定され、Mがゼロに設定さ
れる(ブロック51)。
次に、ルーティンは次のL点を取り入れ、入力関数s
(i)をiがゼロより小さい値の場合にゼロに等しいと
設定する(ブロック52)。
次のステップは、L点直交的離散Gabor係数をインデ
ックスn=0からN−1まで計算する。これは、定数に
等しいmに対する係数の行(column)を計算する(ブロ
ック53)。次のステップで、mは増分され、インデック
スiはi+ΔMと等しく設定すれる(ブロック54)。
次に、ルーティンは係数が全てのフレームに対して計
算されたかを、ループのインデックスMとループで処理
されるサンプルされた系列中のフレームの数とが等しい
かを検査することにより判断する。ルーティンが終了し
ていない場合、ルーティンはブロック55Aに分岐し、こ
こでMを増分してからブロック52に戻る。信号系列が終
了している場合は、アルゴリズムは各係数のエネルギ分
布を第2B図に示されているルーティンに従って計算する
ことに進む。
各係数の周囲のエネルギ分布に基づいて、スペクトル
を計算するために、アルゴリズムは第2A図に示されてい
るルーティンにより生成される係数を受け取る(ブロッ
ク56)。次のステップは、インデックスmを−L/ΔM+
1に設定する(ブロック57)。次は、インデックスnを
ゼロに設定する(ブロック58)。
次に、係数Cm,nに対する更新領域を係数の大きさを
基にして見出す(ブロック59)。
係数に対応する更新領域内の各点(i,k)に対して、
エネルギ分布即ちスペクトルを計算する(ブロック6
0)。次に、インデックスnを増分する(ブロック6
1)。次のステップで、インデックスnは値Nに対して
テストされ、行中の全ての係数の分布が計算されたを決
める。nがNに等しくない場合、アルゴリズムはブロッ
ク59に戻る。nがNに等しい場合、そのときは行は完了
しており、アルゴリズムはブロック63に進み、インデッ
クスmを増分する。次のステップ(ブロック64)で、ア
ルゴリズムはインデックスmが最大値Mに達したかを決
める。達していない場合、そのときアルゴリズムは次の
係数のフレームを処理するためにブロック58にループす
る。mがMに等しい場合、そのときアルゴリズムは完了
し、エネルギ分布値はGaborスペクトルを出力するため
に使用される(ブロック65)。最後に、ルーティンは停
止する(ブロック66)。
第6図は、本発明によるGaborスペクトル分析装置を
示している。Gaborスペクトル分析装置は、時間により
変化する周波数スペクトルにより特徴付けられる入力信
号を分析するのに用いられる装置である。第6図に示さ
れている装置は、典型的な音声処理アルゴリズムで遭遇
するような0〜20kHの範囲の信号に対して特に有用であ
る。しかし、装置は様々な周波数領域で動作するように
設計することができる。このため、同じ基本的な装置を
様々な周波数範囲を有する広い範囲の入力信号に対して
適用できる。
システムは、折返し除去フィルタとA/D変換器を有す
る、変換器/フィルタ100を含んでいる。折返し除去フ
ィルタは様々なデータ収集システムの構成要素である。
折返し除去フィルタの目的は、連続している信号をデジ
タル化する場合に折返しの問題がないように、不要な高
い周波数成分を除去することである。サンプリング理論
によれば、装置の周波数範囲の最高周波数はfとする
と、折返し除去フィルタは、f以上の周波数を実質的に
取り除くように設計される。完全なローパス・フィルタ
を設計することができないため、折返し除去フィルタ
は、fより少し高いカットオフ周波数であるように設計
されることが多い。サンプリング周波数(連続の信号が
デジタル化されるレート)は、少なくても最高周波数の
2倍でなければならない。第6図に示されている分析装
置は、最高速度は20kHであるが、最小サンプリング・レ
ートは40kHである。安全を見越して、48kHのオーディオ
・ボードのサンプリングが使用され、これにより0〜20
kH間には折返しがないことを補償している。
変換器/フィルタ100の出力は、線101に信号s(i)
として供給される。記憶装置102は、入力信号系列を後
で使用するのに記憶するために使用されている。獲得さ
れ、記憶されている信号に対しては、変換器/フィルタ
100はこの応用には必要がない。記憶装置102は、アナロ
グ信号を表しているデジタル信号系列の供給源として作
動する。
線101上の離散データは入力バッファ103へ供給され
る。入力バッファは入力バッファ制御装置により制御さ
れる。入力バッファ103の目的は、この後の処理に対し
て十分な点(L)を蓄積することである。各点Lのブロ
ックはフレームと呼ばれる。
入力バッファ制御装置は、入力バッファを制御する。
その主要な機能は入力バッファを保持し、ルーティンの
その後の処理に対して、オーバーラッピング・データ・
アドレス等を供給することである。選択された時間領域
の増分ΔMに依存して、入力バッファ制御装置は、ΔM
点を現在のフレームに追加し最初のΔM点を取り除い
て、新しいフレームを形成する。
L点データフレームはバッファ103の出力をして線105
上に得られる。線105上のフレームは、離散Gabor変換エ
ンジン106に供給される。エンジン106は、第2A図に示さ
れるアルゴリズムに従ってGabor変換を実行する。エン
ジン106は、一次元の入力信号s(i)をインデックス
(i,k)に対する二次元の時間−周波数領域のGabor係数
に展開する。
Gabor変換エンジン106の出力は、フレームmのインデ
ックスnに対する周波数増分の信号系列に対する係数C
m,nの配列である。この係数の配列は、レジスタ配列の
様な係数バッファ107に供給される。係数バッファ107か
ら、係数の配列が擬似Wigner−Ville分布エンジン108に
供給される。エンジン108は、第2B図に示されているル
ーティンに従ってクロス項削除干渉(cross term delet
ed interference)により擬似Wigner−Ville分布を実行
する。出力は各係数Cm,n周囲のエネルギのWigner分布
で定義されたスペクトルである。係数の配列に対するス
ペクトルの集まりが、Gaborスペクトル記録GS(i,k)で
ある。スペクトル記録はスペクトル記録バッファ109
へ、Gaborスペクトル記録バッファ制御装置110の制御で
供給される。Gaborスペクトル記録バッファ制御装置は
スペクトル記録バッファ109を保持し、制御する。
擬似Wigner−Ville分布エンジンに接続されて、ルッ
クアップ・テーブルがある。ルックアップ・テーブル
は、上記した分布関数で用いられる項目を記憶してい
る。ルックアップ・テーブルを使用して、アルゴリズム
をリアルタイムに近いGaborスペクトル記録の計算に適
しているようにするために、これらの項目の計算がオフ
ラインで計算される。
スペクトル記録バッファ制御装置110は、バッファ109
のデータを、バッファが一杯になると、一時に出力す
る。例えば、データが線112上から表示装置113に供給さ
れ、Gaborスペクトル記録の画像表現が生成され、表示
される。その他に、線112上のデータは、記憶装置114に
後で処理するために供給される。
第6図に示される実施例では、線112上のデータは分
析エンジン115に供給されている。線101上の入力信号系
列s(i)もまた分析装置の入力として供給される。分
析エンジンは、線112上のGaborスペクトル記録データに
応答して入力信号s(i)を区分することができる。そ
して区分された部分の入力信号を独立して分析すること
ができる。
表示されたGaborスペクトル記録の応用や、分析エン
ジン115へのスペクトル記録により生成されるデータ
は、スペクトル分析装置,音声認識システム,音声分析
および合成システム,音声検出システム,水中音響検知
装置,レーダ装置,自動車電話通信,そして電話トーン
検知システムが含まれる。
第7A図と第7B図は、本発明の付表Aのようなルーティ
ンを使用して分析される入力信号を示している。第7A図
において、周波数ホッピング信号が示され、これは離散
時間ウィンドウにおける4つの離散周波数成分がある。
第7B図には、信号/雑音比0dBにおける雑音のあるホッ
ピング信号が示されている。
第8A図は、第7A図の周波数ホッピング信号に対するGa
borスペクトルを示しており、等高線はスペクトル上の
時間と周波数の位置における増加するエネルギ値を示し
ている。このように、第8A図に示すグラフが三次元で表
されると、等高線を有する各領域は、紙の上にピークと
して表される。
ホッピング信号の4つの周波数成分の各々は、Gabor
スペクトル内の時間で局所化している。このような明確
な成分間の局所化と分離のため、Gaborスペクトルは時
間で変化する周波数成分を有する信号を分析するのに特
に有用である。
第8B図は、第7B図に示す雑音のある周波数ホッピング
信号の4重のオーバーサンプリング率によるGaborペク
トルを示している。このように、雑音が相当含まれてい
る信号でも、Gaborスペクトルは周波数成分を比較的よ
く局所化する。Gaborスペクトルは、時間で変化する入
力信号のスペクトルにおいて従来技術よりよいと思われ
る。
上記の発明の好ましい実施例の記述は、説明や記述の
目的のために行われた。これは、包括的とすること、ま
たは発明を開示した詳細な形態に限定することを意図し
たものではない。多数の修正や変形がこの分野の知識を
有するものにとって明白である。実施例は、この発明の
原理と実際の応用を説明するのに最適であるから選択さ
れた。このため、この発明の色々な実施例と特定の意図
した用途に適した改造が、この分野の当業者にとって可
能である。この発明の範囲は特許請求の範囲で定義され
たものおよびその均等物である。
付表B 補助関数γの計算 Hが階数pの満たされた列を有すると仮定すると、QR
分解を用いることにより、 が得られる。ここでQは直交行列があり、Rは上三角行
列である。Hの最初の列がhT=[h(0),h(1),
…,h(L−1)]である。
式(A.1)式は |h,…|=|r1,1q1,r1,2q1+r2,2q2,…|となる。そ
して、 h=r1,1q1 (A.2) ここでq1はQの第1行である。式(A.1)を Hγ=μ (A,3) に代入して、 を得る。式(A.4)において、xは、 x=(RT-1μ (A.5) により定まり、ベクトルyは任意である。QTQ=1であ
るので、 いくつかの注を順に述べる。
(a)h=γ1,1q1であるので、ウィンドウ関数h
(i)は行列Qxの範囲内である。このため Qy Th=0 (A.7) となる。
(b)複直交補助ウィンドウ関数γは2つの直交ベク
トルQxx+Qyyの和である。このため、 ‖γ‖=‖x‖+‖y‖ (A.8) である。
(c)Qxx=Qy(RT-1μは、複直交の関係で定ま
り、Qyyは、特定の制限に依存している。Qyy=0のと
き、解はγ=Qxx=γminで、最小のエネルギーを有す
る。
Γを最小にすると、γは と等しくなる。
関係式 γTh=xTQx Th+yTQy Th=(Qxx)Th (A.11) と式(A.8)により、式(A.10)は と書くことができる。
明らかに、ξの最大値は、‖y‖がゼロ・ベクトルの
とき生じる。
y=0を式(A.6)に代入すると γ=Qx x=Qx(RT-1μ (A.13) が得られる。
対応する最小自乗誤差は である。式(A.13)は、hに最も類似した最適なγを計
算した解である。解は、‖y‖=0である必要があり、
式(A.8)はγが最小エネルギーを含むことを示唆して
いる。γは最小エネルギー条件 γ=HT(HHT-1μ (A.15) を満足する。

Claims (32)

    (57)【特許請求の範囲】
  1. 【請求項1】入力信号を表すデジタル信号系列の供給源
    と、 前記供給源に接続し、前記系列に応答して直交的離散Ga
    bor変換の係数Cm,nを計算し、前記係数を使用して入力
    信号エネルギの時間変化するスペクトルを計算するプロ
    セッサと を有し、前記スペクトルは前記入力信号を分析するのに
    使用することを特徴とする信号分析装置。
  2. 【請求項2】前記プロセッサは、係数を生成する第1の
    手段と、スペクトルを計算する第2の手段とを有するこ
    とを特徴とする請求の範囲第1項記載の信号分析装置。
  3. 【請求項3】前記プロセッサは、スペクトルの計算に使
    用する項目のテーブルを格納する記憶手段を含むことを
    特徴とする請求の範囲第2項記載の信号分析装置。
  4. 【請求項4】前記プロセッサは、クロス項を削除した
    (cross−term deleted)Wigner−Ville分布を計算する
    ことを特徴とする請求の範囲第1項記載の信号分析装
    置。
  5. 【請求項5】前記プロセッサは、局所化したガウス離散
    ウィンドウ関数hおよびhと同様な離散補助関数γを基
    に、直交的離散Gabor変換を行うことにより、係数Cm,n
    計算する計算エンジンと、各係数Cm,nに対するエネル
    ギ分布関数を計算する計算エンジンとを含み、エネルギ
    分布関数は、Gabor数の大きさを基にした第1の項目
    と、Gabor変換の局所化した離散的ウィンドウ関数に従
    って変化する第2の項目とを含むことを特徴とする請求
    の範囲第1項記載の信号分析装置。
  6. 【請求項6】前記プロセッサは、特定の係数Cm,nに対
    するエネルギ分布関数の第2の項目に用いられる値の予
    め計算したテーブルを格納し、特定の係数の座標m,nの
    近傍の局所化したウィンドウ関数における係数の座標i,
    kに応答してアクセスすることができる記憶手段を含む
    ことを特徴とする請求の範囲第5項記載の信号分析装
    置。
  7. 【請求項7】前記プロセッサは、第1の項目に応答して
    各係数Cm,nに対する更新領域を定め、各係数に対する
    分布は定められた更新領域内で計算されることを特徴と
    する請求の範囲第5項記載の信号分析装置。
  8. 【請求項8】前記供給源から供給された系列を受け取る
    ために接続され、プロセッサに重なったフレームの集合
    として系列を供給する入力バッファを含むことを特徴と
    する請求の範囲第1項記載の信号分析装置。
  9. 【請求項9】前記プロセッサに接続され、スペクトルを
    処理する第2のプロセッサを含むことを特徴とする請求
    の範囲第1項記載の信号分析装置。
  10. 【請求項10】前記第2プロセッサは、前記スペクトル
    に応答して入力信号を区分する手段を含むことを特徴と
    する請求の範囲第9項に記載の信号分析装置。
  11. 【請求項11】前記第2のプロセッサは、前記スペクト
    ルに応答して入力信号の分析に用いるためにスペクトル
    を表す画像を表示する手段を含むことを特徴とする請求
    の範囲第9項に記載の信号分析装置。
  12. 【請求項12】前記プロセッサから時間変化するスペク
    トルを受け取り、時間変化するスペクトルを前記第2の
    プロセッサに供給するために接続された出力バッファ記
    憶システムを含むことを特徴とする請求の範囲第9項に
    記載の信号分析装置。
  13. 【請求項13】前記供給源は、アナログ・デジタル変換
    器を含むことを特徴とする請求の範囲第1項に記載の信
    号分析装置。
  14. 【請求項14】前記供給源は、折返し除去フィルタを含
    むことを特徴とする請求の範囲第13項に記載の信号分析
    装置。
  15. 【請求項15】入力端を有し、該入力端に供給された入
    力信号を表すデジタル信号系列を生成する変換器と、 該変換器で生成された系列を受け取るように接続され、
    重なったフレームの集合を供給する入力記憶手段と、 該入力記憶手段に接続され、局所化したガウス離散ウィ
    ンドウ関数hおよびhと同様な離散補助関数γにより離
    散Gabor変換を行うことにより、前記重なったフレーム
    の集合に応答して直交的離散Gabor変換の係数Cm,nを計
    算する第1の手段と、 該第1の手段に接続され、クロス項を削除した(cross
    −term deleted)Wigner−Ville分布を計算することに
    より、前記係数に応答して入力信号エネルギの時間変化
    するスペクトルを生成する第2の手段と、 該第2の手段に接続され、スペクトルを処理する分析手
    段と を有することを特徴とする信号分析装置。
  16. 【請求項16】hが約0.5(L/12)以下の分散σ
    有する(Lはウィンドウmの長さ)ガウス関数であるこ
    とを特徴とする請求の範囲第15項に記載の信号分析装
    置。
  17. 【請求項17】前記第1の手段で計算された特定のGabo
    r変換の係数の分布の寄与が、Gabor変換の特定の係数の
    大きさの基にした第1の項目およびGabor変換の局所化
    したウィンドウ関数に対して変化する第2の項目の関数
    として計算され、前記第2の手段は、前記第2の項目に
    使用される値の予め計算したテーブルを格納し、特定の
    Gabor変換係数の座標m,nの近傍の局所化したウィンドウ
    関数の座標i,kに応答してアクセスすることができる記
    憶手段を含むことを特徴とする請求の範囲第15項に記載
    の信号分析装置。
  18. 【請求項18】前記第2の手段は、前記第1の項目に応
    答して各係数Cm,nに対する更新領域を定め、定められ
    た更新領域内のみで各係数の分布を計算する手段を含む
    ことを特徴とする請求の範囲第17項に記載の信号分析装
    置。
  19. 【請求項19】前記第1の手段から分布データを受け取
    り、分布データを前記第2の手段に供給するために接続
    された出力バッファ記憶システムを含むことを特徴とす
    る請求の範囲第15項に記載の信号分析装置。
  20. 【請求項20】前記変換器は、折返し除去フィルタを含
    むことを特徴とする請求の範囲第15項に記載の信号分析
    装置。
  21. 【請求項21】入力信号の時間で変化するスペクトルを
    計算する方法において、 入力信号を表すデジタル信号系列を得、 前記デジタル信号の系列をサンプリングして長さLの複
    数のウィンドウを定義し、各複数のウィンドウは複数の
    デジタル信号を含み、各ウィンドウは隣接するウィンド
    ウに関連したデジタル信号でΔMの長さシフトし、 サンプリング間隔ΔNとΔMで、各ウィンドウmに対し
    てn=0からN−1まで、複数のウィンドウm=0から
    M−1まで、ΔNΔMがLより小さい場合において、直
    交的離散Gabor変換の係数Cm,nを計算し、 係数Cm,nに応答して入力信号のエネルギの時間で変化
    するスペクトルを計算し、 時間変化するスペクトルに応答し、入力信号を分析する
    ことを特徴とする入力信号の時間で変化するスペクトル
    を計算する方法。
  22. 【請求項22】前記スペクトルを計算するステップは、
    クロス項を削除した(cross−term deleted)Wigner−V
    ille分布を計算することを含むことを特徴とする請求の
    範囲第21項に記載の方法。
  23. 【請求項23】前記直交的離散Gabor変換の係数Cm,n
    計算するステップは、非周期で局所化したウィンドウ関
    数hおよびhと同様な離散補助関数γに基づいて直交的
    離散Gabor変換を実行することを含むことを特徴とする
    請求の範囲第21項に記載の方法。
  24. 【請求項24】hとγは重ならないことを特徴とする請
    求の範囲第23項に記載の方法。
  25. 【請求項25】hはガウス関数であることを特徴とする
    請求の範囲第23項に記載の方法。
  26. 【請求項26】前記ガウス関数hは、約0.5(L/12)
    以下の分散σを有する(Lはウィンドウmの長さ)こ
    とを特徴とする請求の範囲第25項に記載の方法。
  27. 【請求項27】γは最小自乗誤差の意味でhと同様であ
    ることを特徴とする請求の範囲第23項に記載の方法。
  28. 【請求項28】前記スペクトルを計算するステップは、 時間変化するスペクトルの点mΔM,nΔNの周囲の領域
    内の各点i,kに対応する項を有するルックアップ・テー
    ブルを含み、項は時間変化するスペクトルを計算するの
    に使用する予め計算した値を格納しており、 前記ルックアップ・テーブル内の値を使用して時間変化
    するスペクトルを計算する ことを特徴とする請求の範囲第21項に記載の方法。
  29. 【請求項29】前記時間変化するスペクトルを計算する
    ステップは、 各係数Cm,nに対して、係数の大きさに応答して係数の
    周囲に更新領域を定め、定められた更新領域内で各係数
    に対して時間変化するスペクトルを計算することを特徴
    とする請求の範囲第21項に記載の方法。
  30. 【請求項30】前記時間変化するスペクトルを計算する
    ステップは、 各係数Cm,nに対して、係数の大きさに応答して係数の
    周囲に更新領域を定め、時間変化するスペクトルの点m
    ΔM,nΔNの周囲の領域内の各点i,kに対応する項を有す
    るルックアップ・テーブルを含み、項はクロス項を削除
    した(cross−term deleted)Wigner−Ville分布を計算
    するのに使用する予め計算した値を格納しており、ルッ
    クアップ・テーブルは最大の更新領域を有する係数のク
    ロス項を削除した(cross−term deleted)Wigner−Vil
    le分布を計算するために十分な項目を含んでおり、 ルックアップ・テーブル内の値を用いて時間変化するス
    ペクトルを計算することを特徴とする請求の範囲第21項
    に記載の方法。
  31. 【請求項31】前記分析するステップは、スペクトルの
    画像表現を生成することを含むことを特徴とする請求の
    範囲第21項に記載の方法。
  32. 【請求項32】前記分析するステップは、スペクトルに
    応答して入力信号系列を時間上で区分することを含むこ
    とを特徴とする請求の範囲第21項に記載の方法。
JP5516728A 1992-03-17 1993-03-17 時間により変化するスペクトル分析に関する方法および装置 Expired - Lifetime JP2697957B2 (ja)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
US07/851,725 US5353233A (en) 1992-03-17 1992-03-17 Method and apparatus for time varying spectrum analysis
US851,725 1992-03-17

Publications (2)

Publication Number Publication Date
JPH07505949A JPH07505949A (ja) 1995-06-29
JP2697957B2 true JP2697957B2 (ja) 1998-01-19

Family

ID=25311508

Family Applications (1)

Application Number Title Priority Date Filing Date
JP5516728A Expired - Lifetime JP2697957B2 (ja) 1992-03-17 1993-03-17 時間により変化するスペクトル分析に関する方法および装置

Country Status (6)

Country Link
US (1) US5353233A (ja)
EP (1) EP0632899B1 (ja)
JP (1) JP2697957B2 (ja)
AT (1) ATE193769T1 (ja)
DE (1) DE69328828T2 (ja)
WO (1) WO1993019378A1 (ja)

Families Citing this family (26)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5736631A (en) * 1996-03-11 1998-04-07 Turf Diagnostics & Design, Inc. Turf impact analysis system
JP3552837B2 (ja) * 1996-03-14 2004-08-11 パイオニア株式会社 周波数分析方法及び装置並びにこれを用いた複数ピッチ周波数検出方法及び装置
JP3266819B2 (ja) * 1996-07-30 2002-03-18 株式会社エイ・ティ・アール人間情報通信研究所 周期信号変換方法、音変換方法および信号分析方法
US5852567A (en) * 1996-07-31 1998-12-22 Hughes Electronics Corporation Iterative time-frequency domain transform method for filtering time-varying, nonstationary wide band signals in noise
US5845241A (en) * 1996-09-04 1998-12-01 Hughes Electronics Corporation High-accuracy, low-distortion time-frequency analysis of signals using rotated-window spectrograms
CA2298789C (en) * 1997-08-14 2004-06-22 Hendry Mechanical Works Electric arc monitoring systems
US5952957A (en) * 1998-05-01 1999-09-14 The United States Of America As Represented By The Secretary Of The Navy Wavelet transform of super-resolutions based on radar and infrared sensor fusion
US6772077B1 (en) 1998-08-10 2004-08-03 Hendry Mechanical Works Electric arc monitoring systems
US6301572B1 (en) 1998-12-02 2001-10-09 Lockheed Martin Corporation Neural network based analysis system for vibration analysis and condition monitoring
US6944318B1 (en) 1999-01-15 2005-09-13 Citicorp Development Center, Inc. Fast matching systems and methods for personal identification
US6434515B1 (en) 1999-08-09 2002-08-13 National Instruments Corporation Signal analyzer system and method for computing a fast Gabor spectrogram
US6754620B1 (en) * 2000-03-29 2004-06-22 Agilent Technologies, Inc. System and method for rendering data indicative of the performance of a voice activity detector
US6366862B1 (en) 2000-04-19 2002-04-02 National Instruments Corporation System and method for analyzing signals generated by rotating machines
US6810341B2 (en) 2000-04-19 2004-10-26 National Instruments Corporation Time varying harmonic analysis including determination of order components
US6332116B1 (en) 2000-04-19 2001-12-18 National Instruments Corporation System and method for analyzing signals of rotating machines
US6850552B2 (en) 2000-12-01 2005-02-01 Klaus Bibl Iterative precision spectrum analysis
DE10110208A1 (de) * 2001-03-02 2002-09-19 Siemens Production & Logistics Verfahren zum Spezifizieren, Ausführen und Analysieren von Verfahrensabläufen beim Erkennen
US7158901B2 (en) * 2004-07-12 2007-01-02 Tektronix, Inc. FFT accelerated display mode
US7958488B2 (en) * 2005-08-16 2011-06-07 National Instruments Corporation Virtual testing in a development environment
US8185316B2 (en) * 2007-05-25 2012-05-22 Prime Geoscience Corporation Time-space varying spectra for seismic processing
US9164131B2 (en) * 2010-05-13 2015-10-20 Tektronix, Inc. Signal recognition and triggering using computer vision techniques
WO2011150069A2 (en) * 2010-05-25 2011-12-01 The General Hospital Corporation Apparatus, systems, methods and computer-accessible medium for spectral analysis of optical coherence tomography images
RU2449298C1 (ru) * 2010-12-23 2012-04-27 ОАО "Концерн "Океанприбор" Способ определения энергии помехи
TWI438416B (zh) * 2011-04-14 2014-05-21 Delta Electronics Inc 使用連續遞移轉換之信號分析系統及方法
EP3153872A1 (en) * 2015-10-06 2017-04-12 ABB Schweiz AG Method and system for determining phasor components of a periodic waveform
CN109633271A (zh) * 2019-01-17 2019-04-16 长沙理工大学 基于变分模态分解和维格纳威尔分布的行波时频分析方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4837578A (en) * 1981-10-29 1989-06-06 California Institute Of Technology Apparatus and method for range detection using the analytic signal identified from the received signal
US4894795A (en) * 1987-04-28 1990-01-16 The United States Of America As Represented By The Secretary Of The Navy High-resolution technique for time-frequency signal analysis using modified wigner-ville analysis
US5046504A (en) * 1989-02-01 1991-09-10 Corazonix Corporation Method and apparatus for analyzing and interpreting electrocardiograms using spectro-temporal mapping
US5291560A (en) * 1991-07-15 1994-03-01 Iri Scan Incorporated Biometric personal identification system based on iris analysis

Also Published As

Publication number Publication date
US5353233A (en) 1994-10-04
DE69328828T2 (de) 2001-02-15
ATE193769T1 (de) 2000-06-15
EP0632899A4 (en) 1996-09-11
EP0632899B1 (en) 2000-06-07
EP0632899A1 (en) 1995-01-11
JPH07505949A (ja) 1995-06-29
WO1993019378A1 (en) 1993-09-30
DE69328828D1 (de) 2000-07-13

Similar Documents

Publication Publication Date Title
JP2697957B2 (ja) 時間により変化するスペクトル分析に関する方法および装置
Gokhale et al. Time domain signal analysis using wavelet packet decomposition approach
US5029509A (en) Musical synthesizer combining deterministic and stochastic waveforms
EP0388104B1 (en) Method for speech analysis and synthesis
Wise et al. Maximum likelihood pitch estimation
US6266003B1 (en) Method and apparatus for signal processing for time-scale and/or pitch modification of audio signals
Hinich et al. Detection of non-Gaussian signals in non-Gaussian noise using the bispectrum
US6804640B1 (en) Signal noise reduction using magnitude-domain spectral subtraction
Favero Compound wavelets: wavelets for speech recognition
Lim et al. Classification of underwater transient signals using mfcc feature vector
Frost Power-spectrum estimation
Sluyter et al. A novel method for pitch extraction from speech and a hardware model applicable to vocoder systems
Al Fajar et al. Analysis of dft and fft signal transformation with hamming window in labview
JP2867769B2 (ja) 音響測定方法およびその装置
Curtis et al. An investigation of several frequency-domain processing methods for enhancing the intelligibility of speech in wideband random noise
Ainsleigh et al. A B-wavelet-based noise-reduction algorithm
CN113655534B (zh) 基于多线性奇异值张量分解核磁共振fid信号噪声抑制方法
CN112505413B (zh) 一种时频分析方法和系统
Fahy et al. Fast Fourier Transforms and power spectra in LabVIEW
Antweiler et al. Simulation of time variant room impulse responses
Byrne et al. An approximation theoretic approach to maximum entropy spectral analysis
Alrubei et al. A new accurate estimator of the frequency using the three-point interpolation of DFT samples
Jesus Acoustic pressure and particle motion power spectrum estimation with Matlab®
JP3321841B2 (ja) ピッチ周波数推定方法及びその装置
Marple et al. Travels through the time-frequency zone: advanced Doppler ultrasound processing techniques

Legal Events

Date Code Title Description
R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20070919

Year of fee payment: 10

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20080919

Year of fee payment: 11

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20090919

Year of fee payment: 12

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20100919

Year of fee payment: 13

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20100919

Year of fee payment: 13

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20110919

Year of fee payment: 14

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20110919

Year of fee payment: 14

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20120919

Year of fee payment: 15

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20130919

Year of fee payment: 16