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
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
- G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
- G01R23/16—Spectrum 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)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
- Radar Systems Or Details Thereof (AREA)
- Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
Description
【発明の詳細な説明】 限定された著作権放棄 この特許文書の開示の部分には、著作権保護の申し立
てをした素材が含まれている。著作権者は、特許の文書
や開示の複写物が、米国特許商標庁の出願書類や記録に
表れるように、だれにより行われようと異議を申し立て
ない。しかし、その他の権利はすべて留保している。
てをした素材が含まれている。著作権者は、特許の文書
や開示の複写物が、米国特許商標庁の出願書類や記録に
表れるように、だれにより行われようと異議を申し立て
ない。しかし、その他の権利はすべて留保している。
発明の背景 発明の分野 この発明は、デジタル信号処理に関するものであり、
特に、時間により変化する周波数成分を有する信号の特
性を定めるための方法および装置に関するものである。
特に、時間により変化する周波数成分を有する信号の特
性を定めるための方法および装置に関するものである。
関連する技術の説明 デジタル信号処理のための典型的なシステムは、入力
アナログ信号をデジタル・サンプル信号系列への変換、
デジタル・サンプル信号系列をデジタル信号に対して配
列やブロックレベルで計算を行うデジタル信号処理エン
ジンへの供給および計算結果に対する分析を含んでい
る。この分野は、音声処理、レーダーやソナーの画像認
識、水中音波探査器その他の技術に急速に拡大してい
る。
アナログ信号をデジタル・サンプル信号系列への変換、
デジタル・サンプル信号系列をデジタル信号に対して配
列やブロックレベルで計算を行うデジタル信号処理エン
ジンへの供給および計算結果に対する分析を含んでい
る。この分野は、音声処理、レーダーやソナーの画像認
識、水中音波探査器その他の技術に急速に拡大してい
る。
離散フーリェ変換（Discrete Fourier Transform DF
T）と高速フーリェ変換（Fast Fourier Transform FF
T）によるDFTの実現は、デジタル信号処理システムにお
ける基本的な処理技術である。これらの技術は、処理さ
れている信号に対して、時間により周波数成分がどのよ
うに変化しているのかを特徴づけることができないとい
う制限がある。
T）と高速フーリェ変換（Fast Fourier Transform FF
T）によるDFTの実現は、デジタル信号処理システムにお
ける基本的な処理技術である。これらの技術は、処理さ
れている信号に対して、時間により周波数成分がどのよ
うに変化しているのかを特徴づけることができないとい
う制限がある。
従来の時間により周波数成分が変化する信号を分析す
るシステムにおいては、所謂短時間フーリェ変換（STF
T）が用いられていた。このアルゴリズムは、入力信号
の多数の短いウィンドのフーリエ変換を計算することに
基づいている。各短い時間のウィンドごとに変換して結
合することにより、入力信号の時間により変化するスペ
クトルが求まる。しかし、STFT技術を用いて求めた時間
により変化するスペクトルは、よい分解能を有していな
い。また、STFT技術は、計算量が多い。
るシステムにおいては、所謂短時間フーリェ変換（STF
T）が用いられていた。このアルゴリズムは、入力信号
の多数の短いウィンドのフーリエ変換を計算することに
基づいている。各短い時間のウィンドごとに変換して結
合することにより、入力信号の時間により変化するスペ
クトルが求まる。しかし、STFT技術を用いて求めた時間
により変化するスペクトルは、よい分解能を有していな
い。また、STFT技術は、計算量が多い。
Gabor変換は別のデジタル信号処理の手法で、信号の
合成された時間−周波数表現を生成する（D.Gabor“The
ory of Communication"J.IEE（London）,Vol.93 No.II
I,November,1946,pp.429−457を参照）。この応用は、
しかしながら、主として２重直交補助ウィンド関数γを
選択することが難しいために、制限されていた。最近、
有限サイクリック離散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を参照）。
合成された時間−周波数表現を生成する（D.Gabor“The
ory of Communication"J.IEE（London）,Vol.93 No.II
I,November,1946,pp.429−457を参照）。この応用は、
しかしながら、主として２重直交補助ウィンド関数γを
選択することが難しいために、制限されていた。最近、
有限サイクリック離散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を参照）。
発明の概要 本発明は、時間的に変化する周波数成分で特徴づけら
れた入力信号の時間的に変化するスペクトルを生成する
信号分析装置を提供するものである。この発明は、ま
た、時間で変化するスペクトルを基にして信号を分析す
る方法である。
れた入力信号の時間的に変化するスペクトルを生成する
信号分析装置を提供するものである。この発明は、ま
た、時間で変化するスペクトルを基にして信号を分析す
る方法である。
したがって、第１の特徴として、本発明は、入力信号
を表すデジタル信号系列を生成する変換器を有する信号
分析装置を含む。デジタル信号系列は、直交的離散Gabo
r変換係数Ｃm,nをその信号系列に応答して計算し、この
係数に応答して入力信号エレルギの時間で変化するスペ
クトル（ここではGaborスペクトルと記述している）を
計算する第１のプロセッサに供給される。最後に、第２
のプロセッサはこのスペクトルをその後の分析や表示の
ために処理する。時間で変化するスペクトルが有用な一
つの特別な分析は、入力信号を別々の構成要素に分解す
ることである。別々の構成要素は色々な技術を用いて各
々分析される。
を表すデジタル信号系列を生成する変換器を有する信号
分析装置を含む。デジタル信号系列は、直交的離散Gabo
r変換係数Ｃm,nをその信号系列に応答して計算し、この
係数に応答して入力信号エレルギの時間で変化するスペ
クトル（ここではGaborスペクトルと記述している）を
計算する第１のプロセッサに供給される。最後に、第２
のプロセッサはこのスペクトルをその後の分析や表示の
ために処理する。時間で変化するスペクトルが有用な一
つの特別な分析は、入力信号を別々の構成要素に分解す
ることである。別々の構成要素は色々な技術を用いて各
々分析される。
第１のプロセッサは入力バッファを有しており、バッ
ファ制御装置により、入力されるデジタル信号系列を受
信し、係数の計算のために適当なフォーマットで供給さ
れる。同時に、プロセッサは出力バッファを有してお
り、それには時間で変化するスペクトルを表すマトリッ
クスの値を、プロセッサ・エンジンにより、格納してい
る。
ファ制御装置により、入力されるデジタル信号系列を受
信し、係数の計算のために適当なフォーマットで供給さ
れる。同時に、プロセッサは出力バッファを有してお
り、それには時間で変化するスペクトルを表すマトリッ
クスの値を、プロセッサ・エンジンにより、格納してい
る。
デジタル信号プロセッサは、ソフトウェアで制御され
た汎用コンピュータかハードウェアで実現されて、係数
と時間で変化するスペクトルを計算するようにプログラ
ムされている。そのルーティンは、次のステップを基に
している。
た汎用コンピュータかハードウェアで実現されて、係数
と時間で変化するスペクトルを計算するようにプログラ
ムされている。そのルーティンは、次のステップを基に
している。
入力信号を表すデジタル信号系列を取り入れるステッ
プ、 デジタル信号系列をサンプリングして、長さＬの複数
のウィンドウを定義し、複数のウィンドウの各々は、複
数のデジタル信号を含み、隣接のウィンドウに関連する
ΔＭの長さのデジタル信号でシフトするステップ、 サンプリング間隔ΔＮとΔＭの直交的離散Gabor変換
係数Ｃm,nを、ｎ＝０からＮ−１まで各ウィンドウｍに
対してと、複数のウィンドウｍ＝０からＭ−１まで、Δ
ＮΔＭがＬより小さい値に対して、計算するステップ、 係数Ｃm,nに応答して、入力信号のエネルギの時間で
変化するスペクトルを計算するステップ、 入力信号を時間で変化するスペクトルに応答して分析
するステップ 直交的離散Gabor変換の係数は、非周期的な局所離散
ウィンドウ関数ｈおよびｈと同様でｈとγは重ならない
ような離散補助関数γを用いて計算される。時間で変化
するスペクトルは、クロス項を削除した（cross−term
deleted）Wigner−Ville分布を用いて計算される。
プ、 デジタル信号系列をサンプリングして、長さＬの複数
のウィンドウを定義し、複数のウィンドウの各々は、複
数のデジタル信号を含み、隣接のウィンドウに関連する
ΔＭの長さのデジタル信号でシフトするステップ、 サンプリング間隔ΔＮとΔＭの直交的離散Gabor変換
係数Ｃm,nを、ｎ＝０からＮ−１まで各ウィンドウｍに
対してと、複数のウィンドウｍ＝０からＭ−１まで、Δ
ＮΔＭがＬより小さい値に対して、計算するステップ、 係数Ｃm,nに応答して、入力信号のエネルギの時間で
変化するスペクトルを計算するステップ、 入力信号を時間で変化するスペクトルに応答して分析
するステップ 直交的離散Gabor変換の係数は、非周期的な局所離散
ウィンドウ関数ｈおよびｈと同様でｈとγは重ならない
ような離散補助関数γを用いて計算される。時間で変化
するスペクトルは、クロス項を削除した（cross−term
deleted）Wigner−Ville分布を用いて計算される。
本発明の別の特徴は、離散ウィンドウ関数は、約0.5
（L/12）２以下の分散σ２を有するガウス関数（Ｌは計
算で使用したウィンドウの長さ）である。
（L/12）２以下の分散σ２を有するガウス関数（Ｌは計
算で使用したウィンドウの長さ）である。
また本発明の別の特徴は、スペクトルが、係数を基に
したエネルギ分布関数で使用した予め計算した項目を有
するルックアップ・テーブルを使用して計算することで
ある。それぞれの係数に対して、スペクトル上の係数の
周囲の更新領域は、係数の大きさに対応して定められて
いる。エネルギ分布は、各係数の更新領域内で計算され
る。予め計算された項目のルックアップ・テーブルは、
最大の大きさの係数の更新領域をカバーするのに十分な
サイズを有するように計算される。この様に、定められ
た大きさのルックアップ・テーブルは、与えられたGabo
rスペクトルエンジンに対して使用することができる。
したエネルギ分布関数で使用した予め計算した項目を有
するルックアップ・テーブルを使用して計算することで
ある。それぞれの係数に対して、スペクトル上の係数の
周囲の更新領域は、係数の大きさに対応して定められて
いる。エネルギ分布は、各係数の更新領域内で計算され
る。予め計算された項目のルックアップ・テーブルは、
最大の大きさの係数の更新領域をカバーするのに十分な
サイズを有するように計算される。この様に、定められ
た大きさのルックアップ・テーブルは、与えられたGabo
rスペクトルエンジンに対して使用することができる。
本発明の信号分析装置および方法は、音声処理、音や
電磁波の通信、音や電磁波の監視、レーダー、ソナー、
およびその他の信号処理の環境で遭遇する時間で変化す
る周波数成分で特徴づけられる信号に対して特に有用で
ある。例えば、複雑に信号のグループの中で、ある周波
数スペクトルを有する信号の始まりと終わりの時間を特
定する能力は、即ち他の複数の人が存在する部屋内で特
定の人がしゃべり始めるときに必要であるが、本発明の
Gaborスペクトル分布装置や方法により大幅に改良され
る。
電磁波の通信、音や電磁波の監視、レーダー、ソナー、
およびその他の信号処理の環境で遭遇する時間で変化す
る周波数成分で特徴づけられる信号に対して特に有用で
ある。例えば、複雑に信号のグループの中で、ある周波
数スペクトルを有する信号の始まりと終わりの時間を特
定する能力は、即ち他の複数の人が存在する部屋内で特
定の人がしゃべり始めるときに必要であるが、本発明の
Gaborスペクトル分布装置や方法により大幅に改良され
る。
本発明の他の特徴と効果は、明細書、図面、特許請求
の範囲をよく読むことで理解される。
の範囲をよく読むことで理解される。
図面の簡単な説明 第１図は、本発明の信号分析装置のブロック図であ
る。
る。
第2A図および第2B図は、本発明のGaborスペクトルを
計算する方法をしめすフローチャートである。
計算する方法をしめすフローチャートである。
第３図は、Gaborスペクトルの計算で使用されるデジ
タル信号のフレームの構成を示す図である。
タル信号のフレームの構成を示す図である。
第４図は、Gaborスペクトルの計算で使用されるサン
プルを行う分布を示す説明図である。
プルを行う分布を示す説明図である。
第５図は、本発明によるGaborスペクトルのエネルギ
分布に計算に使用する更新領域の例を示す図である。
分布に計算に使用する更新領域の例を示す図である。
第６図は、本発明のGaborスペクトルを生成する信号
分析装置のより詳細なブロック図である。
分析装置のより詳細なブロック図である。
第7A図および第7B図は、それぞれ周波数ホッピング信
号および雑音のあるホッピング信号で、時間で変化する
周波数成分を有する入力信号を示した図である。
号および雑音のあるホッピング信号で、時間で変化する
周波数成分を有する入力信号を示した図である。
第8A図および第8B図は、それぞれ本発明を第7A図およ
び第7B図に示した入力信号に適用して生成されるGabor
スペクトルを示すグラフの図である。
び第7B図に示した入力信号に適用して生成されるGabor
スペクトルを示すグラフの図である。
第９図は、本発明の信号分析装置の応用に用いられる
局所化したウィンドウ関数ｈと離散補助関数γを示す図
である。
局所化したウィンドウ関数ｈと離散補助関数γを示す図
である。
発明を実施するための最良の形態 本発明の実施例を図面を用いて詳細に説明する。
第１図は、本発明の信号分析装置を説明するブロック
図である。分析装置は、符号10で参照されるが、アナロ
グ／デジタル変換器11、デジタル信号プロセッサ12、デ
ジタル信号プロセッサに接続された制御部13、Gaborス
ペクトル分析装置14、そして表示装置15を有する。入力
信号ｓ（ｔ）が変換器や他の入力源から線16に供給され
る。アナログ／デジタル変換器11は入力信号ｓ（ｔ）を
表すデジタル値の系列ｓ（ｉ）を発生する。系列ｓ
（ｉ）はデジタル信号プロセッサ12に供給される。デジ
タル信号プロセッサは、デジタル信号処理（DSP）機能
を実行するようにプログラムされた汎用コンピュータを
使用することで、またはこの目的に適合した特別のデジ
タル信号処理（DSP）を行うエンジンで実施される。デ
ジタル信号プロセッサ12は離散Gabor変換（DGT）を、入
力信号ｓ（ｉ）のGabor係数を計算し、この係数に応答
しクロス項を削除した（cross−term deleted）Wigner
−Ville分布（PWDV）を基にして時間および周波数のエ
ネルギ分布を計算することにより実行する。エネルギ分
布GS（i,s）はデジタル信号プロセッサの出力から線17
上に供給される。エネルギ分布GS（i,s）はGaborペクト
ル分析装置14で分析され、表示装置15に供給するために
エネルギ分布の画像表示を生成する。また、Gaborスペ
クトル分析装置14は、他の分析機能、例えば、音声処
理、伝送、画像認識等で使用するために、時間で入力信
号の成分を分解すること等を行うことができる。スペク
トル分析装置14は、プログラムされた汎用コンピュータ
やこの目的のために設計された専用装置である。分析装
置14はデジタル信号プロセッサ12と同じプロセッサ上に
実現してもよい。
図である。分析装置は、符号10で参照されるが、アナロ
グ／デジタル変換器11、デジタル信号プロセッサ12、デ
ジタル信号プロセッサに接続された制御部13、Gaborス
ペクトル分析装置14、そして表示装置15を有する。入力
信号ｓ（ｔ）が変換器や他の入力源から線16に供給され
る。アナログ／デジタル変換器11は入力信号ｓ（ｔ）を
表すデジタル値の系列ｓ（ｉ）を発生する。系列ｓ
（ｉ）はデジタル信号プロセッサ12に供給される。デジ
タル信号プロセッサは、デジタル信号処理（DSP）機能
を実行するようにプログラムされた汎用コンピュータを
使用することで、またはこの目的に適合した特別のデジ
タル信号処理（DSP）を行うエンジンで実施される。デ
ジタル信号プロセッサ12は離散Gabor変換（DGT）を、入
力信号ｓ（ｉ）のGabor係数を計算し、この係数に応答
しクロス項を削除した（cross−term deleted）Wigner
−Ville分布（PWDV）を基にして時間および周波数のエ
ネルギ分布を計算することにより実行する。エネルギ分
布GS（i,s）はデジタル信号プロセッサの出力から線17
上に供給される。エネルギ分布GS（i,s）はGaborペクト
ル分析装置14で分析され、表示装置15に供給するために
エネルギ分布の画像表示を生成する。また、Gaborスペ
クトル分析装置14は、他の分析機能、例えば、音声処
理、伝送、画像認識等で使用するために、時間で入力信
号の成分を分解すること等を行うことができる。スペク
トル分析装置14は、プログラムされた汎用コンピュータ
やこの目的のために設計された専用装置である。分析装
置14はデジタル信号プロセッサ12と同じプロセッサ上に
実現してもよい。
制御ブロック13は、１実施例では、付表Ａのようなコ
ンピュータ・プログラムで構成されている。
ンピュータ・プログラムで構成されている。
付表Ａに示されているプログラムは、LABVIEW（Natio
nal Instrumentsの商標）計測ソフトウェアシステムの
ために書かれたソースコード・ルーティンの一例であ
る。明らかに、広範囲の種類の他のソースコードで実装
することができる。
nal Instrumentsの商標）計測ソフトウェアシステムの
ために書かれたソースコード・ルーティンの一例であ
る。明らかに、広範囲の種類の他のソースコードで実装
することができる。
付表Ａは、GaborSpectrum.c,GaborMagnitude.c、そし
てautoWV.cを含む３つのルーティンで構成されている。
GaborSpectrum.cはGabor係数の計算のためのデータを設
定する。ユーザは、スペクトルを計算するためのオーバ
ーサンプリング率，分散，およびウィンドウの長さのよ
うなパラメータを選択してルーティンを初期化する。パ
ラメータに関するより詳細な説明が下記に記載されてい
る。ルーティンGaborSpectrum.cはルーティンGaborMagn
itude.cを呼び出す。そのルーティンは、複素Gabor係数
を計算し、Gaborスペクトルを計算するのに使用する係
数の大きさを出力として供給する。ルーティンautoWV.c
は、係数に応答してエネルギ分布を計算し、Gaborスペ
クトルを出力する。
てautoWV.cを含む３つのルーティンで構成されている。
GaborSpectrum.cはGabor係数の計算のためのデータを設
定する。ユーザは、スペクトルを計算するためのオーバ
ーサンプリング率，分散，およびウィンドウの長さのよ
うなパラメータを選択してルーティンを初期化する。パ
ラメータに関するより詳細な説明が下記に記載されてい
る。ルーティンGaborSpectrum.cはルーティンGaborMagn
itude.cを呼び出す。そのルーティンは、複素Gabor係数
を計算し、Gaborスペクトルを計算するのに使用する係
数の大きさを出力として供給する。ルーティンautoWV.c
は、係数に応答してエネルギ分布を計算し、Gaborスペ
クトルを出力する。
基本的には、Gaborスペクトル・アルゴリズムは、次
の１−８のステップで構成されている。
の１−８のステップで構成されている。
ステップ1: ウィンドウの長さＬを選択する。ある実装では、値Ｌ
は常に自乗である。例えば、Ｌ＝64,128,256,512,…等
である。Ｌを定めた後、周波数分解能が式（１）により
定められる。
は常に自乗である。例えば、Ｌ＝64,128,256,512,…等
である。Ｌを定めた後、周波数分解能が式（１）により
定められる。
ここでfsは変換器のサンプリング周波数である。
Ｌ＝128またはＬ＝256が音声処理に応用する場合に適
当である。
当である。
ステップ2: オーバーサンプリング率を選択する。オーバーサンプ
リング率は、式（２）により定められる。
リング率は、式（２）により定められる。
ここで、ΔＭは時間領域の増分、ΔＮは周波数領域の
増分である。
増分である。
（ｍΔM,nΔＮ）は、Gabor変換で見出されるガウス関
数の座標である。ここで、ｍは０から（SL/ΔＭ−１）
まで、SLはサンプル長、そして、ｎは０から（L/ΔＮ−
１）までである。図４はGabor変換アルゴリズムにおけ
る係数Ｃm,nの位置を示す図である。オーバーサンプリ
ング率として４が選択される。しかし、他の値も可能で
ある。オーバーサンプリング率が大きくなると、計算量
が増加する。オーバーサンプリング率の４が、出力スペ
クトルの分解能と必要な計算量との間の妥協として良い
ことが証明されている。
数の座標である。ここで、ｍは０から（SL/ΔＭ−１）
まで、SLはサンプル長、そして、ｎは０から（L/ΔＮ−
１）までである。図４はGabor変換アルゴリズムにおけ
る係数Ｃm,nの位置を示す図である。オーバーサンプリ
ング率として４が選択される。しかし、他の値も可能で
ある。オーバーサンプリング率が大きくなると、計算量
が増加する。オーバーサンプリング率の４が、出力スペ
クトルの分解能と必要な計算量との間の妥協として良い
ことが証明されている。
ステップ3: 式（３）となるように、ΔM,ΔＮを選択する。
付表Ａの実装では、Ｌ＝128の場合、ΔＭ＝８および
ΔＮ＝４を選択した。しかし、他の選択も同様に適して
いる。ΔＭとΔＮの選択で、Gaborスペクトルの時間軸
上および周波数軸上のそれぞれの分解能を定まる。
ΔＮ＝４を選択した。しかし、他の選択も同様に適して
いる。ΔＭとΔＮの選択で、Gaborスペクトルの時間軸
上および周波数軸上のそれぞれの分解能を定まる。
ステップ4: L,ΔM,ΔＮが与えられると、ウィンドウ関数の分散σ
２は（４）の式で定められる。
２は（４）の式で定められる。
これは、ウィンドウ関数ｈに対する局所化したガウス
関数を確保する。分散を大きくするとよりよい周波数領
域の分解能が得られる。分散を小さくするとよりよい時
間領域の分解能が得られる。
関数を確保する。分散を大きくするとよりよい周波数領
域の分解能が得られる。分散を小さくするとよりよい時
間領域の分解能が得られる。
ステップ5: Ｌとσ２が与えられると、ガウス・ウィンドウ関数ｈ
［Ｌ］が与えられる。そして、離散補助関数γ［Ｌ］
が、付表Ｂに示されているように、γが最小自乗誤差に
関してｈと同様に定れられる。最適なシステムでは、γ
とｈは示された領域の制限内で重ならない。補助関数γ
［Ｌ］とウィンドウ関数ｈ［Ｌ］の例が第９図に示され
ている。
［Ｌ］が与えられる。そして、離散補助関数γ［Ｌ］
が、付表Ｂに示されているように、γが最小自乗誤差に
関してｈと同様に定れられる。最適なシステムでは、γ
とｈは示された領域の制限内で重ならない。補助関数γ
［Ｌ］とウィンドウ関数ｈ［Ｌ］の例が第９図に示され
ている。
ステップ１−５は、準備的なステップである。実際の
スペクトル計算はステップ６以下で行われる。
スペクトル計算はステップ６以下で行われる。
ステップ6: 第３図は、どのようにして連続的な信号が離散化する
のか、そしてどのようにしてＬ点フレームの離散信号ｓ
（ｉ）の系列に適応するのかを示している。あるフレー
ムのＬ点が求められると、Ｌ点に対するGabor係数の大
きさを以下に従って計算する。
のか、そしてどのようにしてＬ点フレームの離散信号ｓ
（ｉ）の系列に適応するのかを示している。あるフレー
ムのＬ点が求められると、Ｌ点に対するGabor係数の大
きさを以下に従って計算する。
＊は共役複素数を示す。
ｍΔＭの値は入力データの使用され方を示す。次の例
は、この過程を説明するためである。
は、この過程を説明するためである。
ｍ＝−L/ΔＭ＋1,s（−Ｌ＋ΔＭ）,s（−Ｌ＋Δ＋
１），…,s（ΔＭ−１）が用いられる。
１），…,s（ΔＭ−１）が用いられる。
ｍ＝0,s（０）,s（１）…,s（Ｌ−１）が用いられる； ｍ＝1,s（ΔＭ）,s（ΔＭ＋２），…,s（ΔＭ＋Ｌ−
１）が用いられる； ｍ＝2,s（２ΔＭ）,s（２ΔＭ＋２），…,s（２ΔＭ＋
Ｌ−１）が用いられる； … 一般的に、 ｓ（ｍΔＭ）,s（smΔＭ＋１），…,s（ｍΔＭ＋Ｌ−
１）が用いられる。
１）が用いられる； ｍ＝2,s（２ΔＭ）,s（２ΔＭ＋２），…,s（２ΔＭ＋
Ｌ−１）が用いられる； … 一般的に、 ｓ（ｍΔＭ）,s（smΔＭ＋１），…,s（ｍΔＭ＋Ｌ−
１）が用いられる。
各々、係数Ｃm,nがｎ＝0,1,・・・Ｎ−１の場合に対
して計算される。
して計算される。
示されているように、この計算は、FFTの順序の形式
を有している。
を有している。
ステップ7: 係数Ｃm,nがｓ（ｉ）に対して計算された後、信号の
スペクトルがクロス項を削除した（cross−term delete
d）Wigner−Ville分布を用いて計算される。Gabor数Ｃ
m,nは疎に位置しているので、スペクトルの完全な再構
築には、係数Ｃm,n間の点の補間が必要となる。Wigner
−Ville分布のクロス項の干渉（interference）を削除
後、信号のエネルギは以下の式で計算される。
スペクトルがクロス項を削除した（cross−term delete
d）Wigner−Ville分布を用いて計算される。Gabor数Ｃ
m,nは疎に位置しているので、スペクトルの完全な再構
築には、係数Ｃm,n間の点の補間が必要となる。Wigner
−Ville分布のクロス項の干渉（interference）を削除
後、信号のエネルギは以下の式で計算される。
ステップ8: （ｍΔM,nΔＮ）における各Gabor係数|Cm,n｜に対し
て、スペクトル値GS8（i,k）がステップ７の式で計算さ
れ、前のGS8（i,k）と累算される。即ち、第５図に示さ
れているように、点（i,k）におけるエネルギは、係数C
22,C21,C31,C32そしてその他の分布の累算された寄与で
ある。計算機言語の表現を用いて、その累算は以下のよ
うに表される。
て、スペクトル値GS8（i,k）がステップ７の式で計算さ
れ、前のGS8（i,k）と累算される。即ち、第５図に示さ
れているように、点（i,k）におけるエネルギは、係数C
22,C21,C31,C32そしてその他の分布の累算された寄与で
ある。計算機言語の表現を用いて、その累算は以下のよ
うに表される。
式（８）における定数４は正規化係数で、特定の応用
により変化するものである。
により変化するものである。
上記の項は、（ｍΔM,nΔＮ）におけるGabor係数から
離れた点（i,k）では、急速に減少するので、このステ
ップにおける反復は、点（ｍΔM,nΔＮ）の周囲の小さ
い領域に限定される。正確な領域の大きさは、応用に合
わせて調整される。また、上記指数の項は毎回計算する
必要はない。この項は予め計算され、ルックアップ・テ
ーブルに格納される。
離れた点（i,k）では、急速に減少するので、このステ
ップにおける反復は、点（ｍΔM,nΔＮ）の周囲の小さ
い領域に限定される。正確な領域の大きさは、応用に合
わせて調整される。また、上記指数の項は毎回計算する
必要はない。この項は予め計算され、ルックアップ・テ
ーブルに格納される。
スペクトル計算上の重要な問題は、いつ計算を停止す
るか即ち更新領域の大きさである。Ｃm,nの隣接する点
へのエネルギ寄与は、指数関数的に減少するので、（ｍ
ΔM,nΔＮ）から遠く離れた点を計算する必要はない。
このため、各Gabor係数Ｃm,nに対して更新領域を定める
必要がある。更新領域は、第５図に示されるように、ス
ペクトルが計算される必要がある楕円として定義され
る。
るか即ち更新領域の大きさである。Ｃm,nの隣接する点
へのエネルギ寄与は、指数関数的に減少するので、（ｍ
ΔM,nΔＮ）から遠く離れた点を計算する必要はない。
このため、各Gabor係数Ｃm,nに対して更新領域を定める
必要がある。更新領域は、第５図に示されるように、ス
ペクトルが計算される必要がある楕円として定義され
る。
ｉ′＝ｉ−ｍΔＭ、ｋ′＝ｋ−ｎΔＮとすると、式
（９）は以下のように表される。
（９）は以下のように表される。
これは、第５図に60のような２次元の楕円である。
任意の係数に対する更新領域は、以下のように定義さ
れる。
れる。
ここで、|Cm,n|maxは、全てのGabor係数の最大の大き
さである。
さである。
したがって、各係数に対する更新領域は、以下ので表
される。
される。
最大の更新領域は、|Cm,n|2＝|Cm,n|max 2の場合であ
る。このとき、以下の式が導き出せる。
る。このとき、以下の式が導き出せる。
このように、ｉ′maxとｊ′maxは、ルックアップ・テ
ーブルのサイズを与える。正のインデックスと値のみを
格納する必要がある。これは、式（B.2）はｉ′とｋ′
に対して偶関数であるからである。
ーブルのサイズを与える。正のインデックスと値のみを
格納する必要がある。これは、式（B.2）はｉ′とｋ′
に対して偶関数であるからである。
一旦、テーブルが作成されると、式（８）計算はすぐ
できる。対応する項目をテーブル内で探し、その値と|C
m,n|2とを掛ける。計算上のコストが高い指数関数の計
算がこの様にして避けることができる。
できる。対応する項目をテーブル内で探し、その値と|C
m,n|2とを掛ける。計算上のコストが高い指数関数の計
算がこの様にして避けることができる。
第2A図と第2B図は、好ましい実施例に従ってGaborス
ペクトルを計算する基本的な手法が示されている。
ペクトルを計算する基本的な手法が示されている。
プロセッサをまず初期化するルーティンが含まれる
（ブロック50）。この初期化は、上記に概略したステッ
プ１−５を実行することが含まれる。
（ブロック50）。この初期化は、上記に概略したステッ
プ１−５を実行することが含まれる。
次に、Gabor係数が第2A図のブロック51,52,53,54、そ
して55で実現されているアルゴリズムに従って計算され
る。
して55で実現されているアルゴリズムに従って計算され
る。
特に、インデックスｉが−Ｌ＋ΔＭに設定され、イン
デックスｍがL/ΔＭ＋１に設定され、Ｍがゼロに設定さ
れる（ブロック51）。
デックスｍがL/ΔＭ＋１に設定され、Ｍがゼロに設定さ
れる（ブロック51）。
次に、ルーティンは次のＬ点を取り入れ、入力関数ｓ
（ｉ）をｉがゼロより小さい値の場合にゼロに等しいと
設定する（ブロック52）。
（ｉ）をｉがゼロより小さい値の場合にゼロに等しいと
設定する（ブロック52）。
次のステップは、Ｌ点直交的離散Gabor係数をインデ
ックスｎ＝０からＮ−１まで計算する。これは、定数に
等しいｍに対する係数の行（column）を計算する（ブロ
ック53）。次のステップで、ｍは増分され、インデック
スｉはｉ＋ΔＭと等しく設定すれる（ブロック54）。
ックスｎ＝０からＮ−１まで計算する。これは、定数に
等しいｍに対する係数の行（column）を計算する（ブロ
ック53）。次のステップで、ｍは増分され、インデック
スｉはｉ＋ΔＭと等しく設定すれる（ブロック54）。
次に、ルーティンは係数が全てのフレームに対して計
算されたかを、ループのインデックスＭとループで処理
されるサンプルされた系列中のフレームの数とが等しい
かを検査することにより判断する。ルーティンが終了し
ていない場合、ルーティンはブロック55Aに分岐し、こ
こでＭを増分してからブロック52に戻る。信号系列が終
了している場合は、アルゴリズムは各係数のエネルギ分
布を第2B図に示されているルーティンに従って計算する
ことに進む。
算されたかを、ループのインデックスＭとループで処理
されるサンプルされた系列中のフレームの数とが等しい
かを検査することにより判断する。ルーティンが終了し
ていない場合、ルーティンはブロック55Aに分岐し、こ
こでＭを増分してからブロック52に戻る。信号系列が終
了している場合は、アルゴリズムは各係数のエネルギ分
布を第2B図に示されているルーティンに従って計算する
ことに進む。
各係数の周囲のエネルギ分布に基づいて、スペクトル
を計算するために、アルゴリズムは第2A図に示されてい
るルーティンにより生成される係数を受け取る（ブロッ
ク56）。次のステップは、インデックスｍを−L/ΔＭ＋
１に設定する（ブロック57）。次は、インデックスｎを
ゼロに設定する（ブロック58）。
を計算するために、アルゴリズムは第2A図に示されてい
るルーティンにより生成される係数を受け取る（ブロッ
ク56）。次のステップは、インデックスｍを−L/ΔＭ＋
１に設定する（ブロック57）。次は、インデックスｎを
ゼロに設定する（ブロック58）。
次に、係数Ｃm,nに対する更新領域を係数の大きさを
基にして見出す（ブロック59）。
基にして見出す（ブロック59）。
係数に対応する更新領域内の各点（i,k）に対して、
エネルギ分布即ちスペクトルを計算する（ブロック6
0）。次に、インデックスｎを増分する（ブロック6
1）。次のステップで、インデックスｎは値Ｎに対して
テストされ、行中の全ての係数の分布が計算されたを決
める。ｎがＮに等しくない場合、アルゴリズムはブロッ
ク59に戻る。ｎがＮに等しい場合、そのときは行は完了
しており、アルゴリズムはブロック63に進み、インデッ
クスｍを増分する。次のステップ（ブロック64）で、ア
ルゴリズムはインデックスｍが最大値Ｍに達したかを決
める。達していない場合、そのときアルゴリズムは次の
係数のフレームを処理するためにブロック58にループす
る。ｍがＭに等しい場合、そのときアルゴリズムは完了
し、エネルギ分布値はGaborスペクトルを出力するため
に使用される（ブロック65）。最後に、ルーティンは停
止する（ブロック66）。
エネルギ分布即ちスペクトルを計算する（ブロック6
0）。次に、インデックスｎを増分する（ブロック6
1）。次のステップで、インデックスｎは値Ｎに対して
テストされ、行中の全ての係数の分布が計算されたを決
める。ｎがＮに等しくない場合、アルゴリズムはブロッ
ク59に戻る。ｎがＮに等しい場合、そのときは行は完了
しており、アルゴリズムはブロック63に進み、インデッ
クスｍを増分する。次のステップ（ブロック64）で、ア
ルゴリズムはインデックスｍが最大値Ｍに達したかを決
める。達していない場合、そのときアルゴリズムは次の
係数のフレームを処理するためにブロック58にループす
る。ｍがＭに等しい場合、そのときアルゴリズムは完了
し、エネルギ分布値はGaborスペクトルを出力するため
に使用される（ブロック65）。最後に、ルーティンは停
止する（ブロック66）。
第６図は、本発明によるGaborスペクトル分析装置を
示している。Gaborスペクトル分析装置は、時間により
変化する周波数スペクトルにより特徴付けられる入力信
号を分析するのに用いられる装置である。第６図に示さ
れている装置は、典型的な音声処理アルゴリズムで遭遇
するような０〜20kHの範囲の信号に対して特に有用であ
る。しかし、装置は様々な周波数領域で動作するように
設計することができる。このため、同じ基本的な装置を
様々な周波数範囲を有する広い範囲の入力信号に対して
適用できる。
示している。Gaborスペクトル分析装置は、時間により
変化する周波数スペクトルにより特徴付けられる入力信
号を分析するのに用いられる装置である。第６図に示さ
れている装置は、典型的な音声処理アルゴリズムで遭遇
するような０〜20kHの範囲の信号に対して特に有用であ
る。しかし、装置は様々な周波数領域で動作するように
設計することができる。このため、同じ基本的な装置を
様々な周波数範囲を有する広い範囲の入力信号に対して
適用できる。
システムは、折返し除去フィルタとA/D変換器を有す
る、変換器／フィルタ100を含んでいる。折返し除去フ
ィルタは様々なデータ収集システムの構成要素である。
折返し除去フィルタの目的は、連続している信号をデジ
タル化する場合に折返しの問題がないように、不要な高
い周波数成分を除去することである。サンプリング理論
によれば、装置の周波数範囲の最高周波数はｆとする
と、折返し除去フィルタは、ｆ以上の周波数を実質的に
取り除くように設計される。完全なローパス・フィルタ
を設計することができないため、折返し除去フィルタ
は、ｆより少し高いカットオフ周波数であるように設計
されることが多い。サンプリング周波数（連続の信号が
デジタル化されるレート）は、少なくても最高周波数の
２倍でなければならない。第６図に示されている分析装
置は、最高速度は20kHであるが、最小サンプリング・レ
ートは40kHである。安全を見越して、48kHのオーディオ
・ボードのサンプリングが使用され、これにより０〜20
kH間には折返しがないことを補償している。
る、変換器／フィルタ100を含んでいる。折返し除去フ
ィルタは様々なデータ収集システムの構成要素である。
折返し除去フィルタの目的は、連続している信号をデジ
タル化する場合に折返しの問題がないように、不要な高
い周波数成分を除去することである。サンプリング理論
によれば、装置の周波数範囲の最高周波数はｆとする
と、折返し除去フィルタは、ｆ以上の周波数を実質的に
取り除くように設計される。完全なローパス・フィルタ
を設計することができないため、折返し除去フィルタ
は、ｆより少し高いカットオフ周波数であるように設計
されることが多い。サンプリング周波数（連続の信号が
デジタル化されるレート）は、少なくても最高周波数の
２倍でなければならない。第６図に示されている分析装
置は、最高速度は20kHであるが、最小サンプリング・レ
ートは40kHである。安全を見越して、48kHのオーディオ
・ボードのサンプリングが使用され、これにより０〜20
kH間には折返しがないことを補償している。
変換器／フィルタ100の出力は、線101に信号ｓ（ｉ）
として供給される。記憶装置102は、入力信号系列を後
で使用するのに記憶するために使用されている。獲得さ
れ、記憶されている信号に対しては、変換器／フィルタ
100はこの応用には必要がない。記憶装置102は、アナロ
グ信号を表しているデジタル信号系列の供給源として作
動する。
として供給される。記憶装置102は、入力信号系列を後
で使用するのに記憶するために使用されている。獲得さ
れ、記憶されている信号に対しては、変換器／フィルタ
100はこの応用には必要がない。記憶装置102は、アナロ
グ信号を表しているデジタル信号系列の供給源として作
動する。
線101上の離散データは入力バッファ103へ供給され
る。入力バッファは入力バッファ制御装置により制御さ
れる。入力バッファ103の目的は、この後の処理に対し
て十分な点（Ｌ）を蓄積することである。各点Ｌのブロ
ックはフレームと呼ばれる。
る。入力バッファは入力バッファ制御装置により制御さ
れる。入力バッファ103の目的は、この後の処理に対し
て十分な点（Ｌ）を蓄積することである。各点Ｌのブロ
ックはフレームと呼ばれる。
入力バッファ制御装置は、入力バッファを制御する。
その主要な機能は入力バッファを保持し、ルーティンの
その後の処理に対して、オーバーラッピング・データ・
アドレス等を供給することである。選択された時間領域
の増分ΔＭに依存して、入力バッファ制御装置は、ΔＭ
点を現在のフレームに追加し最初のΔＭ点を取り除い
て、新しいフレームを形成する。
その主要な機能は入力バッファを保持し、ルーティンの
その後の処理に対して、オーバーラッピング・データ・
アドレス等を供給することである。選択された時間領域
の増分ΔＭに依存して、入力バッファ制御装置は、ΔＭ
点を現在のフレームに追加し最初のΔＭ点を取り除い
て、新しいフレームを形成する。
Ｌ点データフレームはバッファ103の出力をして線105
上に得られる。線105上のフレームは、離散Gabor変換エ
ンジン106に供給される。エンジン106は、第2A図に示さ
れるアルゴリズムに従ってGabor変換を実行する。エン
ジン106は、一次元の入力信号ｓ（ｉ）をインデックス
（i,k）に対する二次元の時間−周波数領域のGabor係数
に展開する。
上に得られる。線105上のフレームは、離散Gabor変換エ
ンジン106に供給される。エンジン106は、第2A図に示さ
れるアルゴリズムに従ってGabor変換を実行する。エン
ジン106は、一次元の入力信号ｓ（ｉ）をインデックス
（i,k）に対する二次元の時間−周波数領域のGabor係数
に展開する。
Gabor変換エンジン106の出力は、フレームｍのインデ
ックスｎに対する周波数増分の信号系列に対する係数Ｃ
m,nの配列である。この係数の配列は、レジスタ配列の
様な係数バッファ107に供給される。係数バッファ107か
ら、係数の配列が擬似Wigner−Ville分布エンジン108に
供給される。エンジン108は、第2B図に示されているル
ーティンに従ってクロス項削除干渉（cross term delet
ed interference）により擬似Wigner−Ville分布を実行
する。出力は各係数Ｃm,n周囲のエネルギのWigner分布
で定義されたスペクトルである。係数の配列に対するス
ペクトルの集まりが、Gaborスペクトル記録GS（i,k）で
ある。スペクトル記録はスペクトル記録バッファ109
へ、Gaborスペクトル記録バッファ制御装置110の制御で
供給される。Gaborスペクトル記録バッファ制御装置は
スペクトル記録バッファ109を保持し、制御する。
ックスｎに対する周波数増分の信号系列に対する係数Ｃ
m,nの配列である。この係数の配列は、レジスタ配列の
様な係数バッファ107に供給される。係数バッファ107か
ら、係数の配列が擬似Wigner−Ville分布エンジン108に
供給される。エンジン108は、第2B図に示されているル
ーティンに従ってクロス項削除干渉（cross term delet
ed interference）により擬似Wigner−Ville分布を実行
する。出力は各係数Ｃm,n周囲のエネルギのWigner分布
で定義されたスペクトルである。係数の配列に対するス
ペクトルの集まりが、Gaborスペクトル記録GS（i,k）で
ある。スペクトル記録はスペクトル記録バッファ109
へ、Gaborスペクトル記録バッファ制御装置110の制御で
供給される。Gaborスペクトル記録バッファ制御装置は
スペクトル記録バッファ109を保持し、制御する。
擬似Wigner−Ville分布エンジンに接続されて、ルッ
クアップ・テーブルがある。ルックアップ・テーブル
は、上記した分布関数で用いられる項目を記憶してい
る。ルックアップ・テーブルを使用して、アルゴリズム
をリアルタイムに近いGaborスペクトル記録の計算に適
しているようにするために、これらの項目の計算がオフ
ラインで計算される。
クアップ・テーブルがある。ルックアップ・テーブル
は、上記した分布関数で用いられる項目を記憶してい
る。ルックアップ・テーブルを使用して、アルゴリズム
をリアルタイムに近いGaborスペクトル記録の計算に適
しているようにするために、これらの項目の計算がオフ
ラインで計算される。
スペクトル記録バッファ制御装置110は、バッファ109
のデータを、バッファが一杯になると、一時に出力す
る。例えば、データが線112上から表示装置113に供給さ
れ、Gaborスペクトル記録の画像表現が生成され、表示
される。その他に、線112上のデータは、記憶装置114に
後で処理するために供給される。
のデータを、バッファが一杯になると、一時に出力す
る。例えば、データが線112上から表示装置113に供給さ
れ、Gaborスペクトル記録の画像表現が生成され、表示
される。その他に、線112上のデータは、記憶装置114に
後で処理するために供給される。
第６図に示される実施例では、線112上のデータは分
析エンジン115に供給されている。線101上の入力信号系
列ｓ（ｉ）もまた分析装置の入力として供給される。分
析エンジンは、線112上のGaborスペクトル記録データに
応答して入力信号ｓ（ｉ）を区分することができる。そ
して区分された部分の入力信号を独立して分析すること
ができる。
析エンジン115に供給されている。線101上の入力信号系
列ｓ（ｉ）もまた分析装置の入力として供給される。分
析エンジンは、線112上のGaborスペクトル記録データに
応答して入力信号ｓ（ｉ）を区分することができる。そ
して区分された部分の入力信号を独立して分析すること
ができる。
表示されたGaborスペクトル記録の応用や、分析エン
ジン115へのスペクトル記録により生成されるデータ
は、スペクトル分析装置，音声認識システム，音声分析
および合成システム，音声検出システム，水中音響検知
装置，レーダ装置，自動車電話通信，そして電話トーン
検知システムが含まれる。
ジン115へのスペクトル記録により生成されるデータ
は、スペクトル分析装置，音声認識システム，音声分析
および合成システム，音声検出システム，水中音響検知
装置，レーダ装置，自動車電話通信，そして電話トーン
検知システムが含まれる。
第7A図と第7B図は、本発明の付表Ａのようなルーティ
ンを使用して分析される入力信号を示している。第7A図
において、周波数ホッピング信号が示され、これは離散
時間ウィンドウにおける４つの離散周波数成分がある。
第7B図には、信号／雑音比0dBにおける雑音のあるホッ
ピング信号が示されている。
ンを使用して分析される入力信号を示している。第7A図
において、周波数ホッピング信号が示され、これは離散
時間ウィンドウにおける４つの離散周波数成分がある。
第7B図には、信号／雑音比0dBにおける雑音のあるホッ
ピング信号が示されている。
第8A図は、第7A図の周波数ホッピング信号に対するGa
borスペクトルを示しており、等高線はスペクトル上の
時間と周波数の位置における増加するエネルギ値を示し
ている。このように、第8A図に示すグラフが三次元で表
されると、等高線を有する各領域は、紙の上にピークと
して表される。
borスペクトルを示しており、等高線はスペクトル上の
時間と周波数の位置における増加するエネルギ値を示し
ている。このように、第8A図に示すグラフが三次元で表
されると、等高線を有する各領域は、紙の上にピークと
して表される。
ホッピング信号の４つの周波数成分の各々は、Gabor
スペクトル内の時間で局所化している。このような明確
な成分間の局所化と分離のため、Gaborスペクトルは時
間で変化する周波数成分を有する信号を分析するのに特
に有用である。
スペクトル内の時間で局所化している。このような明確
な成分間の局所化と分離のため、Gaborスペクトルは時
間で変化する周波数成分を有する信号を分析するのに特
に有用である。
第8B図は、第7B図に示す雑音のある周波数ホッピング
信号の４重のオーバーサンプリング率によるGaborペク
トルを示している。このように、雑音が相当含まれてい
る信号でも、Gaborスペクトルは周波数成分を比較的よ
く局所化する。Gaborスペクトルは、時間で変化する入
力信号のスペクトルにおいて従来技術よりよいと思われ
る。
信号の４重のオーバーサンプリング率によるGaborペク
トルを示している。このように、雑音が相当含まれてい
る信号でも、Gaborスペクトルは周波数成分を比較的よ
く局所化する。Gaborスペクトルは、時間で変化する入
力信号のスペクトルにおいて従来技術よりよいと思われ
る。
上記の発明の好ましい実施例の記述は、説明や記述の
目的のために行われた。これは、包括的とすること、ま
たは発明を開示した詳細な形態に限定することを意図し
たものではない。多数の修正や変形がこの分野の知識を
有するものにとって明白である。実施例は、この発明の
原理と実際の応用を説明するのに最適であるから選択さ
れた。このため、この発明の色々な実施例と特定の意図
した用途に適した改造が、この分野の当業者にとって可
能である。この発明の範囲は特許請求の範囲で定義され
たものおよびその均等物である。
目的のために行われた。これは、包括的とすること、ま
たは発明を開示した詳細な形態に限定することを意図し
たものではない。多数の修正や変形がこの分野の知識を
有するものにとって明白である。実施例は、この発明の
原理と実際の応用を説明するのに最適であるから選択さ
れた。このため、この発明の色々な実施例と特定の意図
した用途に適した改造が、この分野の当業者にとって可
能である。この発明の範囲は特許請求の範囲で定義され
たものおよびその均等物である。
付表Ｂ 補助関数γの計算 Ｈが階数ｐの満たされた列を有すると仮定すると、QR
分解を用いることにより、 が得られる。ここでＱは直交行列があり、Ｒは上三角行
列である。Ｈの最初の列がhT＝［ｈ（０）,h（１），
…,h（Ｌ−１）］である。
分解を用いることにより、 が得られる。ここでＱは直交行列があり、Ｒは上三角行
列である。Ｈの最初の列がhT＝［ｈ（０）,h（１），
…,h（Ｌ−１）］である。
式（A.1）式は |h,…｜＝|r1,1q1,r1,2q1＋ｒ2,2q2,…｜となる。そ
して、 ｈ＝ｒ1,1q1 （A.2） ここでq1はＱの第１行である。式（A.1）を Ｈγ＝μ （A,3） に代入して、 を得る。式（A.4）において、ｘは、 ｘ＝（RT）-1μ （A.5） により定まり、ベクトルｙは任意である。QTQ＝１であ
るので、 いくつかの注を順に述べる。
して、 ｈ＝ｒ1,1q1 （A.2） ここでq1はＱの第１行である。式（A.1）を Ｈγ＝μ （A,3） に代入して、 を得る。式（A.4）において、ｘは、 ｘ＝（RT）-1μ （A.5） により定まり、ベクトルｙは任意である。QTQ＝１であ
るので、 いくつかの注を順に述べる。
（ａ）ｈ＝γ1,1q1であるので、ウィンドウ関数ｈ
（ｉ）は行列Qxの範囲内である。このため Qy Th＝０ （A.7） となる。
（ｉ）は行列Qxの範囲内である。このため Qy Th＝０ （A.7） となる。
（ｂ）複直交補助ウィンドウ関数γは２つの直交ベク
トルQxx＋Qyyの和である。このため、 ‖γ‖２＝‖ｘ‖２＋‖ｙ‖２ （A.8） である。
トルQxx＋Qyyの和である。このため、 ‖γ‖２＝‖ｘ‖２＋‖ｙ‖２ （A.8） である。
（ｃ）Qxx＝Qy（RT）-1μは、複直交の関係で定ま
り、Qyyは、特定の制限に依存している。Qyy＝０のと
き、解はγ＝Qxx＝γminで、最小のエネルギーを有す
る。
り、Qyyは、特定の制限に依存している。Qyy＝０のと
き、解はγ＝Qxx＝γminで、最小のエネルギーを有す
る。
Γを最小にすると、γは と等しくなる。
関係式 γTh＝xTQx Th＋yTQy Th＝（Qxx）Th （A.11） と式（A.8）により、式（A.10）は と書くことができる。
明らかに、ξの最大値は、‖ｙ‖がゼロ・ベクトルの
とき生じる。
とき生じる。
ｙ＝０を式（A.6）に代入すると γ＝Qx x＝Qx（RT）-1μ （A.13） が得られる。
対応する最小自乗誤差は である。式（A.13）は、ｈに最も類似した最適なγを計
算した解である。解は、‖ｙ‖＝０である必要があり、
式（A.8）はγが最小エネルギーを含むことを示唆して
いる。γは最小エネルギー条件 γ＝HT（HHT）-1μ （A.15） を満足する。
算した解である。解は、‖ｙ‖＝０である必要があり、
式（A.8）はγが最小エネルギーを含むことを示唆して
いる。γは最小エネルギー条件 γ＝HT（HHT）-1μ （A.15） を満足する。
Claims (32)
- 【請求項１】入力信号を表すデジタル信号系列の供給源
と、 前記供給源に接続し、前記系列に応答して直交的離散Ga
bor変換の係数Ｃm,nを計算し、前記係数を使用して入力
信号エネルギの時間変化するスペクトルを計算するプロ
セッサと を有し、前記スペクトルは前記入力信号を分析するのに
使用することを特徴とする信号分析装置。 - 【請求項２】前記プロセッサは、係数を生成する第１の
手段と、スペクトルを計算する第２の手段とを有するこ
とを特徴とする請求の範囲第１項記載の信号分析装置。 - 【請求項３】前記プロセッサは、スペクトルの計算に使
用する項目のテーブルを格納する記憶手段を含むことを
特徴とする請求の範囲第２項記載の信号分析装置。 - 【請求項４】前記プロセッサは、クロス項を削除した
（cross−term deleted）Wigner−Ville分布を計算する
ことを特徴とする請求の範囲第１項記載の信号分析装
置。 - 【請求項５】前記プロセッサは、局所化したガウス離散
ウィンドウ関数ｈおよびｈと同様な離散補助関数γを基
に、直交的離散Gabor変換を行うことにより、係数Ｃm,n
計算する計算エンジンと、各係数Ｃm,nに対するエネル
ギ分布関数を計算する計算エンジンとを含み、エネルギ
分布関数は、Gabor数の大きさを基にした第１の項目
と、Gabor変換の局所化した離散的ウィンドウ関数に従
って変化する第２の項目とを含むことを特徴とする請求
の範囲第１項記載の信号分析装置。 - 【請求項６】前記プロセッサは、特定の係数Ｃm,nに対
するエネルギ分布関数の第２の項目に用いられる値の予
め計算したテーブルを格納し、特定の係数の座標m,nの
近傍の局所化したウィンドウ関数における係数の座標i,
kに応答してアクセスすることができる記憶手段を含む
ことを特徴とする請求の範囲第５項記載の信号分析装
置。 - 【請求項７】前記プロセッサは、第１の項目に応答して
各係数Ｃm,nに対する更新領域を定め、各係数に対する
分布は定められた更新領域内で計算されることを特徴と
する請求の範囲第５項記載の信号分析装置。 - 【請求項８】前記供給源から供給された系列を受け取る
ために接続され、プロセッサに重なったフレームの集合
として系列を供給する入力バッファを含むことを特徴と
する請求の範囲第１項記載の信号分析装置。 - 【請求項９】前記プロセッサに接続され、スペクトルを
処理する第２のプロセッサを含むことを特徴とする請求
の範囲第１項記載の信号分析装置。 - 【請求項１０】前記第２プロセッサは、前記スペクトル
に応答して入力信号を区分する手段を含むことを特徴と
する請求の範囲第９項に記載の信号分析装置。 - 【請求項１１】前記第２のプロセッサは、前記スペクト
ルに応答して入力信号の分析に用いるためにスペクトル
を表す画像を表示する手段を含むことを特徴とする請求
の範囲第９項に記載の信号分析装置。 - 【請求項１２】前記プロセッサから時間変化するスペク
トルを受け取り、時間変化するスペクトルを前記第２の
プロセッサに供給するために接続された出力バッファ記
憶システムを含むことを特徴とする請求の範囲第９項に
記載の信号分析装置。 - 【請求項１３】前記供給源は、アナログ・デジタル変換
器を含むことを特徴とする請求の範囲第１項に記載の信
号分析装置。 - 【請求項１４】前記供給源は、折返し除去フィルタを含
むことを特徴とする請求の範囲第13項に記載の信号分析
装置。 - 【請求項１５】入力端を有し、該入力端に供給された入
力信号を表すデジタル信号系列を生成する変換器と、 該変換器で生成された系列を受け取るように接続され、
重なったフレームの集合を供給する入力記憶手段と、 該入力記憶手段に接続され、局所化したガウス離散ウィ
ンドウ関数ｈおよびｈと同様な離散補助関数γにより離
散Gabor変換を行うことにより、前記重なったフレーム
の集合に応答して直交的離散Gabor変換の係数Ｃm,nを計
算する第１の手段と、 該第１の手段に接続され、クロス項を削除した（cross
−term deleted）Wigner−Ville分布を計算することに
より、前記係数に応答して入力信号エネルギの時間変化
するスペクトルを生成する第２の手段と、 該第２の手段に接続され、スペクトルを処理する分析手
段と を有することを特徴とする信号分析装置。 - 【請求項１６】ｈが約0.5（L/12）２以下の分散σ２を
有する（Ｌはウィンドウｍの長さ）ガウス関数であるこ
とを特徴とする請求の範囲第15項に記載の信号分析装
置。 - 【請求項１７】前記第１の手段で計算された特定のGabo
r変換の係数の分布の寄与が、Gabor変換の特定の係数の
大きさの基にした第１の項目およびGabor変換の局所化
したウィンドウ関数に対して変化する第２の項目の関数
として計算され、前記第２の手段は、前記第２の項目に
使用される値の予め計算したテーブルを格納し、特定の
Gabor変換係数の座標m,nの近傍の局所化したウィンドウ
関数の座標i,kに応答してアクセスすることができる記
憶手段を含むことを特徴とする請求の範囲第15項に記載
の信号分析装置。 - 【請求項１８】前記第２の手段は、前記第１の項目に応
答して各係数Ｃm,nに対する更新領域を定め、定められ
た更新領域内のみで各係数の分布を計算する手段を含む
ことを特徴とする請求の範囲第17項に記載の信号分析装
置。 - 【請求項１９】前記第１の手段から分布データを受け取
り、分布データを前記第２の手段に供給するために接続
された出力バッファ記憶システムを含むことを特徴とす
る請求の範囲第15項に記載の信号分析装置。 - 【請求項２０】前記変換器は、折返し除去フィルタを含
むことを特徴とする請求の範囲第15項に記載の信号分析
装置。 - 【請求項２１】入力信号の時間で変化するスペクトルを
計算する方法において、 入力信号を表すデジタル信号系列を得、 前記デジタル信号の系列をサンプリングして長さＬの複
数のウィンドウを定義し、各複数のウィンドウは複数の
デジタル信号を含み、各ウィンドウは隣接するウィンド
ウに関連したデジタル信号でΔＭの長さシフトし、 サンプリング間隔ΔＮとΔＭで、各ウィンドウｍに対し
てｎ＝０からＮ−１まで、複数のウィンドウｍ＝０から
Ｍ−１まで、ΔＮΔＭがＬより小さい場合において、直
交的離散Gabor変換の係数Ｃm,nを計算し、 係数Ｃm,nに応答して入力信号のエネルギの時間で変化
するスペクトルを計算し、 時間変化するスペクトルに応答し、入力信号を分析する
ことを特徴とする入力信号の時間で変化するスペクトル
を計算する方法。 - 【請求項２２】前記スペクトルを計算するステップは、
クロス項を削除した（cross−term deleted）Wigner−V
ille分布を計算することを含むことを特徴とする請求の
範囲第21項に記載の方法。 - 【請求項２３】前記直交的離散Gabor変換の係数Ｃm,nを
計算するステップは、非周期で局所化したウィンドウ関
数ｈおよびｈと同様な離散補助関数γに基づいて直交的
離散Gabor変換を実行することを含むことを特徴とする
請求の範囲第21項に記載の方法。 - 【請求項２４】ｈとγは重ならないことを特徴とする請
求の範囲第23項に記載の方法。 - 【請求項２５】ｈはガウス関数であることを特徴とする
請求の範囲第23項に記載の方法。 - 【請求項２６】前記ガウス関数ｈは、約0.5（L/12）２
以下の分散σ２を有する（Ｌはウィンドウｍの長さ）こ
とを特徴とする請求の範囲第25項に記載の方法。 - 【請求項２７】γは最小自乗誤差の意味でｈと同様であ
ることを特徴とする請求の範囲第23項に記載の方法。 - 【請求項２８】前記スペクトルを計算するステップは、 時間変化するスペクトルの点ｍΔM,nΔＮの周囲の領域
内の各点i,kに対応する項を有するルックアップ・テー
ブルを含み、項は時間変化するスペクトルを計算するの
に使用する予め計算した値を格納しており、 前記ルックアップ・テーブル内の値を使用して時間変化
するスペクトルを計算する ことを特徴とする請求の範囲第21項に記載の方法。 - 【請求項２９】前記時間変化するスペクトルを計算する
ステップは、 各係数Ｃm,nに対して、係数の大きさに応答して係数の
周囲に更新領域を定め、定められた更新領域内で各係数
に対して時間変化するスペクトルを計算することを特徴
とする請求の範囲第21項に記載の方法。 - 【請求項３０】前記時間変化するスペクトルを計算する
ステップは、 各係数Ｃm,nに対して、係数の大きさに応答して係数の
周囲に更新領域を定め、時間変化するスペクトルの点ｍ
ΔM,nΔＮの周囲の領域内の各点i,kに対応する項を有す
るルックアップ・テーブルを含み、項はクロス項を削除
した（cross−term deleted）Wigner−Ville分布を計算
するのに使用する予め計算した値を格納しており、ルッ
クアップ・テーブルは最大の更新領域を有する係数のク
ロス項を削除した（cross−term deleted）Wigner−Vil
le分布を計算するために十分な項目を含んでおり、 ルックアップ・テーブル内の値を用いて時間変化するス
ペクトルを計算することを特徴とする請求の範囲第21項
に記載の方法。 - 【請求項３１】前記分析するステップは、スペクトルの
画像表現を生成することを含むことを特徴とする請求の
範囲第21項に記載の方法。 - 【請求項３２】前記分析するステップは、スペクトルに
応答して入力信号系列を時間上で区分することを含むこ
とを特徴とする請求の範囲第21項に記載の方法。
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| US851,725 | 1992-03-17 | ||
| US07/851,725 US5353233A (en) | 1992-03-17 | 1992-03-17 | Method and apparatus for time varying spectrum analysis |
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)
| 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 |
| AU763046B2 (en) * | 1997-08-14 | 2003-07-10 | 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 |
| US6810341B2 (en) | 2000-04-19 | 2004-10-26 | National Instruments Corporation | Time varying harmonic analysis including determination of order components |
| US6366862B1 (en) | 2000-04-19 | 2002-04-02 | National Instruments Corporation | System and method for analyzing signals generated by rotating machines |
| 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 |
| EP2575598A2 (en) * | 2010-05-25 | 2013-04-10 | 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)
| 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 |
-
1992
- 1992-03-17 US US07/851,725 patent/US5353233A/en not_active Expired - Lifetime
-
1993
- 1993-03-17 EP EP93907584A patent/EP0632899B1/en not_active Expired - Lifetime
- 1993-03-17 WO PCT/US1993/002499 patent/WO1993019378A1/en active IP Right Grant
- 1993-03-17 AT AT93907584T patent/ATE193769T1/de not_active IP Right Cessation
- 1993-03-17 JP JP5516728A patent/JP2697957B2/ja not_active Expired - Lifetime
- 1993-03-17 DE DE69328828T patent/DE69328828T2/de not_active Expired - Lifetime
Also Published As
| Publication number | Publication date |
|---|---|
| DE69328828T2 (de) | 2001-02-15 |
| DE69328828D1 (de) | 2000-07-13 |
| EP0632899A4 (en) | 1996-09-11 |
| JPH07505949A (ja) | 1995-06-29 |
| US5353233A (en) | 1994-10-04 |
| EP0632899B1 (en) | 2000-06-07 |
| WO1993019378A1 (en) | 1993-09-30 |
| EP0632899A1 (en) | 1995-01-11 |
| ATE193769T1 (de) | 2000-06-15 |
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 | |
| CN117935826B (zh) | 音频升采样方法、装置、设备及存储介质 | |
| 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® | |
| 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 |