特許 7603304 周波数推定装置、周波数推定方法、及び周波数推定プログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2024-12-12

(45)【発行日】2024-12-20

(54)【発明の名称】周波数推定装置、周波数推定方法、及び周波数推定プログラム

(51)【国際特許分類】

   G01R 23/167 20060101AFI20241213BHJP

   G01R 23/16 20060101ALI20241213BHJP

   G01R 23/173 20060101ALI20241213BHJP

   H04B 1/10 20060101ALI20241213BHJP

   H03H 17/06 20060101ALI20241213BHJP

【ＦＩ】

G01R23/167

G01R23/16 D

G01R23/173 D

H04B1/10 N

H03H17/06 635B

【請求項の数】 4

(21)【出願番号】P 2021002224

(22)【出願日】2021-01-08

(65)【公開番号】P2022107342

(43)【公開日】2022-07-21

【審査請求日】2023-12-11

【新規性喪失の例外の表示】特許法第３０条第２項適用 令和２年５月２１日 ｈｔｔｐｓ：／／ｗｗｗ．ｉｅｉｃｅ．ｏｒｇ／ｋｅｎ／ｉｎｄｅｘ／ｉｅｉｃｅ－ｔｅｃｈｒｅｐ－１２０－３８．ｈｔｍｌ にて公開 令和２年５月２９日 映像情報メディア学会ＩＳＴ／ＭＥ、電子情報通信学会ＭＩ，ＩＥ，ＳＩＰ，ＢｉｏＸ合同研究会 にて公開

(73)【特許権者】

【識別番号】000125370

【氏名又は名称】学校法人東京理科大学

(74)【代理人】

【識別番号】110001519

【氏名又は名称】弁理士法人太陽国際特許事務所

(72)【発明者】

【氏名】相川 直幸

【審査官】島▲崎▼ 純一

(56)【参考文献】

【文献】特開２０００－１８１４７２（ＪＰ，Ａ）

【文献】国際公開第２０１８／１５４７４７（ＷＯ，Ａ１）

【文献】特開平１０－１９７５７５（ＪＰ，Ａ）

【文献】特開平０５－１３６７１０（ＪＰ，Ａ）

【文献】国際公開第２００４／００４２６４（ＷＯ，Ａ１）

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ０１Ｒ ２３／１６７

Ｇ０１Ｒ ２３／１６

Ｇ０１Ｒ ２３／１７３

Ｈ０４Ｂ １／１０

Ｈ０３Ｈ １７／０６

(57)【特許請求の範囲】

【請求項1】

入力信号の周波数を出力するフィルタ部と、
前記入力信号をヒルベルト変換するヒルベルト変換部と、
前記入力信号及び前記ヒルベルト変換部から出力された変換信号から解析信号を生成する解析信号生成部と、
前記解析信号に基づいて、瞬時周波数を算出する瞬時周波数算出部と、
前記入力信号の周波数の偶数倍の周波数成分が除去されるように、前記瞬時周波数に基づいて前記フィルタ部のフィルタ係数を設定するフィルタ係数設定部と、
を備え、
前記変換信号は、前記ヒルベルト変換の次数に基づく誤差を含む信号である
周波数推定装置。

【請求項2】

前記瞬時周波数が、前記入力信号をサンプリングするサンプリング周波数に対応するナイキスト周波数よりも大きい周波数の場合、前記フィルタ係数設定部は、前記瞬時周波数及び前記入力信号の周波数の偶数倍の周波数成分の折り返しノイズが除去されるように前記フィルタ係数を設定する
請求項１記載の周波数推定装置。

【請求項3】

コンピュータが、
入力信号をヒルベルト変換し、
前記入力信号及び前記ヒルベルト変換された変換信号から解析信号を生成し、
前記解析信号に基づいて、瞬時周波数を算出し、
前記入力信号の周波数の偶数倍の周波数成分が除去されるように、前記瞬時周波数に基づいて前記入力信号の周波数を出力するフィルタ部のフィルタ係数を設定し、
前記フィルタ係数に基づいて前記入力信号の周波数を出力する
処理を実行し、
前記変換信号は、前記ヒルベルト変換の次数に基づく誤差を含む信号である
周波数推定方法。

【請求項4】

コンピュータに、
入力信号をヒルベルト変換し、
前記入力信号及び前記ヒルベルト変換された変換信号から解析信号を生成し、
前記解析信号に基づいて、瞬時周波数を算出し、
前記入力信号の周波数の偶数倍の周波数成分が除去されるように、前記瞬時周波数に基づいて前記入力信号の周波数を出力するフィルタ部のフィルタ係数を設定し、
前記フィルタ係数に基づいて前記入力信号の周波数を出力する
処理を実行させ、
前記変換信号は、前記ヒルベルト変換の次数に基づく誤差を含む信号である
周波数推定プログラム。

【発明の詳細な説明】

【技術分野】

【0001】

開示の技術は、周波数推定装置、周波数推定方法、及び周波数推定プログラムに関する。

【背景技術】

【0002】

正弦波は最も基本的な波形の１つであり、その周波数を推定することはレダーやソナー、通信、医学などの分野で研究されてきた問題である。例えば、電子式はかりの分野において、単一周波数を信号を用いる音叉式センサは分解能、安定性、及び再現性が優れており、近年注目を集めている。

【0003】

音叉式センサは、音叉２個をつないだ形状の振動子に荷重を加えると振動数が変化する現象を利用した荷重センサである。金属製の音叉と、これに荷重を伝達するリング構造のみで構成され、荷重によってセンサの振動数が変化する。この周波数を測定し、得られた周波数に対してディジタル処理を行うことでノイズを除去し、重量値に変換する。

【0004】

このような単一周波数測定の主な手法としては、周波数カウンタを用いる手法と瞬時周波数による手法が挙げられる。前者の手法で現在主流となっているのが、被測定信号の１周期を基準クロックを用いて測定し、その逆数から周波数を求めるレシプロカル方式の周波数カウンタである。この方式では、被測定信号の周波数が変化すると計測値である周波数の出力間隔は不等間隔になってしまう。そのため、ディジタル処理を行うためには補間処理等の工夫が必要となる。

【0005】

後者の瞬時周波数は一般的にヒルベルト変換器を用いて求められる。ヒルベルト変換器を用いて瞬時周波数を推定する手法では、サンプルごとに周波数の推定値を出力することができ、周波数の時間変化が計測可能となる。このため、瞬時周波数推定は、回転機の故障診断や非線形及び非定常過程の理解にも用いられている。ヒルベルト変換器の入力を実部、出力を虚部とした複素信号は解析信号と呼ばれ、解析信号における或る瞬間の位相が瞬時位相である。この瞬時位相を時間微分し、２πで割ったものが瞬時周波数となる。

【0006】

ここで、ディジタル信号の場合、位相の時間微分は隣接するサンプル間の差分から計算できるため、周期を求めることなく周波数を推定することが可能となる。従って、得られる周波数はサンプル毎に出力される。また、瞬時周波数は非線形及び非定常過程の詳細なメカニズムを理解するためにも活用されている。

【0007】

実際に用いる解析信号を求めるために使用するヒルベルト変換器は、有限次数で設計する必要がある。そのため、理想のヒルベルト変換器の振幅特性に対して、近似誤差としてリプルを有する。その結果、ヒルベルト変換器を用いて推定された瞬時周波数には、本来の周波数成分以外にリプルに由来する成分が含まれるため、高精度な周波数推定をする際に問題となる。この誤差はヒルベルト変換器の次数を増やすことで抑制することが可能であるが、次数の増加は遅延時間の増大、回路規模の増加を伴うため問題となる。

【0008】

本発明者らは、推定された瞬時周波数に含まれる振動成分の周波数が、入力信号の周波数の偶数倍となることを示し、これを伝送零点を有するＦＩＲ（Ｆｉｎｉｔｅ Ｉｍｐｕｌｓｅ Ｒｅｓｐｏｎｓｅ）フィルタを用いて除去することで、推定精度の向上が達成できることを提案した（非特許文献１参照）。

【先行技術文献】

【非特許文献】

【0009】

【文献】高尾圭祐，名取隆廣，宮田統馬，相川直幸，“有限次数ヒルベルト変換器を用いた高精度瞬時周波数の一推定法，” 第34 回信号処理シンポジウム，vol.B2-2，pp.237-241，Nov. 2019．

【発明の概要】

【発明が解決しようとする課題】

【0010】

しかしながら、非特許文献１記載の方法は、入力信号の周波数が変化するたびにフィルタを再設計する必要がある、という問題があった。

【0011】

開示の技術は、上記事実を考慮して成されたものであり、入力信号の周波数が変化するたびにフィルタを再設計することなく周波数を推定することができる周波数推定装置、周波数推定方法、及び周波数推定プログラムを得ることを目的とする。

【課題を解決するための手段】

【0012】

上記目的を達成するため、第１態様に係る周波数推定装置は、入力信号の周波数を出力するフィルタ部と、前記入力信号をヒルベルト変換するヒルベルト変換部と、前記入力信号及び前記ヒルベルト変換部から出力された変換信号から解析信号を生成する解析信号生成部と、前記解析信号に基づいて、瞬時周波数を算出する瞬時周波数算出部と、前記入力信号の周波数の偶数倍の周波数成分が除去されるように、前記瞬時周波数に基づいて前記フィルタ部のフィルタ係数を設定するフィルタ係数設定部と、を備える。

【0013】

第１態様に係る周波数推定装置において、前記瞬時周波数が、前記入力信号をサンプリングするサンプリング周波数に対応するナイキスト周波数よりも大きい周波数の場合、前記フィルタ係数設定部は、前記瞬時周波数及び前記入力信号の周波数の偶数倍の周波数成分の折り返しノイズが除去されるように前記フィルタ係数を設定するようにしてもよい。

【0014】

第１態様に係る周波数推定装置において、前記変換信号は、前記ヒルベルト変換の次数に基づく誤差を含む信号としてもよい。

【0015】

第２態様に係る周波数推定方法は、コンピュータが、入力信号をヒルベルト変換し、前記入力信号及び前記ヒルベルト変換された変換信号から解析信号を生成し、前記解析信号に基づいて、瞬時周波数を算出し、前記入力信号の周波数の偶数倍の周波数成分が除去されるように、前記瞬時周波数に基づいて前記入力信号の周波数を出力するフィルタ部のフィルタ係数を設定し、前記フィルタ係数に基づいて前記入力信号の周波数を出力する処理を実行する。

【0016】

第３態様に係る周波数推定プログラムは、コンピュータに、入力信号をヒルベルト変換し、前記入力信号及び前記ヒルベルト変換された変換信号から解析信号を生成し、前記解析信号に基づいて、瞬時周波数を算出し、前記入力信号の周波数の偶数倍の周波数成分が除去されるように、前記瞬時周波数に基づいて前記入力信号の周波数を出力するフィルタ部のフィルタ係数を設定し、前記フィルタ係数に基づいて前記入力信号の周波数を出力する処理を実行させる。

【発明の効果】

【0017】

開示の技術によれば、入力信号の周波数が変化するたびにフィルタを再設計することなく周波数を推定することができる、という効果を有する。

【図面の簡単な説明】

【0018】

【図1】ヒルベルト変換器の振幅特性を示す図である。

【図2】図１の特性を有するヒルベルト変換器を用いて推定される瞬時周波数を示す図である。

【図3】４２次のヒルベルト変換器を用いて推定された瞬時周波数のフーリエスペクトルを示す図である。

【図4】伝送零点について説明するための図である。

【図5】周波数推定装置の機能ブロック図である。

【図6】周波数推定処理のフローチャート図である。

【図7】可変フィルタの振幅特性を示す図である。

【図8】可変フィルタの振幅特性を示す図である。

【図9】瞬時周波数の推定結果を示す図である。

【図10】図９の一部を拡大した図である。

【図11】図９の一部を拡大した図である。

【図12】チャープ信号の瞬時周波数の推定結果を示す図である。

【図13】図１２の一部を拡大した図である。

【図14】折り返しノイズについて説明するための図である。

【発明を実施するための形態】

【0019】

以下、本発明の実施形態について図面を参照しながら詳細に説明する。

【0020】

＜開示の技術の概要＞

【0021】

本実施形態では、前述した非特許文献１記載の方法では、入力信号の周波数が変化するたびにフィルタを再設計する必要があるという問題に対して、伝送零点が可変な低域通過ＦＩＲフィルタを用いて振動成分を除去する方法について説明する。

【0022】

以下で説明する伝送零点が可変な低域通過ＦＩＲフィルタは、阻止域特性を変化させることで、入力信号の周波数が変化してもその偶数倍に伝送零点が変化するフィルタである。

【0023】

Ｎ次の直線位相ＦＩＲフィルタを用いたヒルベルト変換器の零位相振幅特性Ｈ_０（ｅ^ｊω）は次式で表される。

【0024】

・・・（１）

【0025】

ここで、ｈ（ｎ）はフィルタ係数である。

【0026】

下記参考文献１に記載されたＲｅｍｅｚアルゴリズムを用いて設計された４２次及び６０次のヒルベルト変換器の振幅特性を図１に示す。ただし、サンプリング周波数は５０ｋＨｚ、低域通過域端周波数、高域通過域端周波数は１．２５ｋＨｚ、２７．７５ｋＨｚとする。図１から明らかなように、有限次数のヒルベルト変換器のリプルは次数が大きいほど小さくなる。

【0027】

(参考文献１)

【0028】

J. McClellan, T. Parks, and L. Rabiner, “A computer program for designing optimum FIR linear phase digital filters,”IEEE Transactions on Audio and Electroacoustics, vol.21, no.6, pp.506-526, Dec. 1973. Conference Name:IEEE Transactions on Audio and Electroacoustics.

【0029】

正弦波の入力信号ｘ（ｔ）は次式で表される。

【0030】

・・・（２）

【0031】

ここで、Ａは振幅、ω_１は角周波数を示す。また、図１に示すように、ある角周波数ωにおける理想の振幅特性と実際の振幅特性との差をδとする。そのとき上記（２）式の入力信号ｘ（ｔ）に対して、有限次数のヒルベルト変換器を用いて得られる解析信号

は、次式で表される。なお、以下では解析信号をｘ＾（ｔ）と表す場合がある。

【0032】

・・・（３）

【0033】

ただし、ｘ_ｈ（ｔ）はヒルベルト変換器の出力である。この解析信号から瞬時位相は次式で表される。

【0034】

・・・（４）

【0035】

そして、瞬時周波数（瞬時角周波数）

は次式で表される。なお、以下ではω~（ｔ）と表す場合がある。

【0036】

・・・（５）

【0037】

上記（５）式より推定される瞬時周波数はヒルベルト変換器のリプルがδ＝０となる角周波数ω_１以外ではｃｏｓ項が残り、真の角周波数ω_１とはならないことが判る。従って、有限次数のヒルベルト変換器を用いて得られる瞬時周波数は、真の角周波数と別に振動成分を有することが判る。上記非特許文献１に記載されているように、真の角周波数ω_１とそれ以外の角周波数との関係は、次式で表される。

【0038】

・・・（６）

【0039】

ただし、ａ_ｎは次式で表される。

【0040】

・・・（７）

【0041】

また、ａ_ｎは、以下の漸化式を用いて求められる。

【0042】

・・・（８）
ただし、

【0043】

上記（６）～（８）式より、推定された瞬時角周波数ω~（ｔ）は、直流成分として入力信号の角周波数ω_１を持ち、δに依存する角周波数がω_１の偶数倍となる余弦成分の和として表現できることがわかる。従って、有限次数ヒルベルト変換器を用いて推定された瞬時角周波数には、δ≠０の場合、入力信号角周波数の偶数倍の角周波数成分が含まれる。よって、有限次数のヒルベルト変換器を用いて高精度な角周波数推定を行うには、推定された瞬時角周波数の偶数倍の角周波数成分を除去すれば良い。

【0044】

例として１８００Ｈｚの入力信号ｘ（ｔ）を入力した時、図１の特性を有するヒルベルト変換器を用いて推定される瞬時周波数を図２に示す。これより、実際に推定された瞬時周波数が真の周波数である１８００Ｈｚ付近で振動していることが判る。さらに、４２次のヒルベルト変換器を用いて推定された瞬時周波数のフーリエスペクトルを図３に示す。このグラフからも推定された瞬時周波数は入力信号の周波数の２倍である３６００Ｈｚ、４倍である７２００Ｈｚの成分を持つことが判る。上記非特許文献１記載の方法では、この振動成分を除去することによって，推定精度の向上を可能としている。しかしながら、このフィルタでは、周波数が変化するたびにフィルタを再設計する必要がある。

【0045】

以下、瞬時角周波数に含まれる偶数倍の振動成分を除去するためのフィルタについて述べる。ヒルベルト変換器に入力される単一正弦波の角周波数は未知であり、回転機の故障診断を行うような場合など、角周波数が変化することも考慮する必要がある。また、単純に伝送零点を有するフィルタをカスケード接続した場合、ある周波数帯で信号が増幅される可能性があり、これはノイズを含む信号を扱う際に問題となる。そこで、本実施形態では、推定値に含まれる振動成分を効率よく除去するためのフィルタとして、図４に示すような伝送零点位置が可変な低域通過ＦＩＲフィルタを用いる。図４に示すようなフィルタの振幅特性を次式で与える。

【0046】

・・・（９）

【0047】

ただし、Ｂは伝送零点の数、Ｍはフィルタ次数の半分の次数である。ここで、ｃ（ｋ、α）はフィルタ係数であり、可変パラメータが変化するとフィルタ係数ｃ（ｋ、α）が変化するため、以下のように可変パラメータを用いたＬ次の多項式として近似する。

【0048】

・・・（１０）

【0049】

ここで、ｇ（ｋ、ｌ）は多項式係数である。また、αは、想定される入力信号の正規化角周波数における下限値であるω'_ｍｉｎから上限値であるω’_ｍａｘの範囲内で連続的に変化し、その偶数倍の位置に伝送零点を置く。上記（９）、（１０）式を行列の形で表すと次式のように表される。

【0050】

・・・（１１）

【0051】

ただし、Ｔは転置を表し、Ｃ、Ａ、Ｐ、Ｇは、それぞれ下記（１２）～（１６）式で表される。

【0052】

・・・（１２）

・・・（１３）

・・・（１４）

・・・（１５）

・・・（１６）

【0053】

ここで、フィルタの理想特性は、図４より次式で定義される。

【0054】

・・・（１７）

【0055】

ただし、ω_ｐ、ω_ｓはそれぞれ通過域端正規化角周波数、阻止域端正規化角周波数である。

【0056】

この時、誤差関数ｅ（ω、α）は次式で表される。

【0057】

ｅ（ω、α）＝Ｄ（ω）－Ａ（ω，α） ・・・（１８）

【0058】

従って、周波数領域におけるフィルタの設計問題は以下で表される。

【0059】

・・・（１９）

【0060】

ただし、Ωは近似帯域の周波数範囲であり、遷移域は評価対象外となるため、Ωは次式で表される。

【0061】

・・・（２０）

【0062】

ここで、計算機を用いて上述の設計問題を解くためにαをＱ点に離散化した場合、上記（１８）式は、次式で表される。

【0063】

ｅ（ω、α_ｑ）＝Ｄ（ω）－Ａ（ω，α_ｑ） ・・・（２１）

【0064】

ただし、１≦ｑ≦Ｑである。ここで、λ（ω，α_ｑ）を正の最大許容誤差とすると、上記（１９）の設計問題は、下記式のように変形することができる。

【0065】

minimaize λ（ω，α_ｑ） ・・・（２２）

【0066】

subject to ｅ（ω、α_ｑ）^２≦λ（ω，α_ｑ） ・・・（２３）

【0067】

さらに、上記（２３）式は次式のように変形できる。

【0068】

λ（ω、α_ｑ）－ｅ（ω，α_ｑ）^２≧０ ・・・（２４）

【0069】

上記（２２）、（２３）式を行列の形で表すと次式で表される。

【0070】

・・・（２５）

【0071】

上記（２５）式より、Γ（ω，α_ｑ）は多項式係数Ｇと正の最大許容誤差λ（ω，α_ｑ）に関して線形であり、半正定値である、さらに、上記（２５）式の近似帯域ＱをＭ分割に、αをＱ分割に離散化した半正定値条件は次式で表される、

【0072】

・・・（２６）

【0073】

ここで、ｍは１≦ｍ≦Ｍの範囲内の値である。従って、上記（２２）、（２３）の最適化問題をω、αについて離散化すると次式のように定式化される。

【0074】

minimaize ｄｘ^Ｔ ・・・（２７）

【0075】

subject to Ｕ（ｘ）≧０ ・・・（２８）

【0076】

ただし、ｄ、ｘは次式で表される。

【0077】

・・・（２９）

・・・（３０）

【0078】

行列Ｕ（ｘ）はｘに関して線形であるため、上記（２７）～（３０）式は半正定値問題である。

【0079】

＜周波数推定装置の構成及び作用＞

【0080】

図５は、周波数推定装置１０の機能構成を示す図である。図５に示すように、周波数推定装置１０は、ヒルベルト変換部１２、解析信号生成部１４、瞬時周波数算出部１６、フィルタ係数設定部１８、及びフィルタ部２０を備える。

【0081】

ヒルベルト変換部１２は、入力信号ｘ（ｔ）をヒルベルト変換する。入力信号ｘ（ｔ）は、上記（２）式で表される。ヒルベルト変換部１２は、前述した有限次数のヒルベルト変換器で構成される。変換信号ｘ_ｈ（ｔ）は、ヒルベルト変換の次数に基づく誤差、すなわち前述した誤差δを含む信号である。

【0082】

解析信号生成部１４は、入力信号ｘ（ｔ）及びヒルベルト変換部１２から出力された変換信号ｘ_ｈ（ｔ）から上記（３）式により解析信号ｘ＾（ｔ）を生成する。

【0083】

瞬時周波数算出部１６は、解析信号生成部１４により生成された解析信号ｘ＾（ｔ）に基づいて、瞬時周波数ω~（ｔ）を上記（５）式により算出する。

【0084】

フィルタ係数設定部１８は、入力信号（ｔ）の周波数の偶数倍の周波数成分が除去されるように、瞬時周波数ω~（ｔ）に基づいて上記（１０）式によりフィルタ部２０のフィルタ係数ｃ（ｋ，α）を設定する。

【0085】

フィルタ部２０は、フィルタ係数設定部１８により設定されたフィルタ係数ｃ（ｋ，α）に基づいて、入力信号ｘ（ｔ）の周波数を出力する。フィルタ部２０は、前述した伝送零点位置が可変な低域通過ＦＩＲフィルタにより構成される。具体的には、振幅特性が、上記（９）式で表される振幅特性Ａ（ω、α）であるＦＩＲフィルタで構成される。

【0086】

周波数推定装置１０に含まれる各機能部は、例えばＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ）等で構成されるが、これに限られない。例えば周波数推定装置１０をＣＰＵ（Central Processing Unit）、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、及び入出力インターフェース（Ｉ／Ｏ）等を備えたコンピュータで構成してもよい。この場合、周波数推定プログラムをＲＯＭに予め記憶しておき、ＣＰＵがＲＯＭから周波数推定プログラムを読み込んで周波数推定処理を実行することにより、各機能部が実現される構成としてもよい。

【0087】

次に、本実施形態の作用として、周波数推定装置１０で実行される周波数推定処理について図６に示すフローチャートを参照して説明する。なお、図６の周波数推定処理は、入力信号ｘ（ｔ）が入力されている間は繰り返し実行される。

【0088】

ステップＳ１００では、入力信号ｘ（ｔ）をヒルベルト変換する。

【0089】

ステップＳ１０２では、入力信号ｘ（ｔ）及びステップＳ１００でヒルベルト変換された変換信号ｘ_ｈ（ｔ）から上記（３）式により解析信号ｘ＾（ｔ）を生成する。

【0090】

ステップＳ１０４では、ステップＳ１０２で生成された解析信号ｘ＾（ｔ）に基づいて、瞬時周波数ω~（ｔ）を上記（５）式により算出する。

【0091】

ステップＳ１０６では、入力信号（ｔ）の周波数の偶数倍の周波数成分が除去されるように、瞬時周波数ω~（ｔ）に基づいて上記（１０）式によりフィルタ処理のフィルタ係数ｃ（ｋ，α）を設定する。

【0092】

ステップＳ１０８では、ステップＳ１０６で設定されたフィルタ係数ｃ（ｋ，α）に基づいて、フィルタ処理を実行し、入力信号ｘ（ｔ）の周波数を出力する。

【0093】

このように、本実施形態では、伝送零点が可変な低域通過ＦＩＲフィルタとして、阻止域特性を変化させることで入力信号の周波数が変化してもその偶数倍に伝送零点が変化するフィルタを用いて入力信号ｘ（ｔ）の周波数を推定する。これにより、入力信号ｘ（ｔ）の周波数が変化するたびにフィルタを再設計することなく周波数を推定することができる。

【0094】

＜実施例＞

【0095】

以下、実施例について説明する。

【0096】

ここでは、前述した伝送零点位置が可変な低域通過ＦＩＲフィルタ（以下、可変フィルタと称する）を用いることで、瞬時周波数を高精度に推定できることをシミュレーションによって示す。本シミュレーションでは入力信号として振幅が１で周波数が１６００Ｈｚから２０００Ｈｚの範囲で不連続に変化する信号を用いた。サンプリング周波数を５０ｋＨｚとし、ヒルベルト変換器の次数は４２次、低域通過域端周波数、高域通過域端周波数はそれぞれ１．２５ｋＨｚ、２７．７５ｋＨｚとした。また、可変フィルタの通過域端周波数、阻止域端周波数はそれぞれ２．５Ｈｚ、３０００Ｈｚとし、各種パラメータとして、次数２Ｎ＝１２、多項式次数Ｌ＝３、伝送零点の数Ｂ＝２、３とし、可変パラメータαの範囲は入力信号の周波数が１６００～２０００Ｈｚとなる範囲とした。このような設計仕様に対し、評価点を通過域で１０等分、阻止域で９９０等分し、可変パラメータの変化に関して１０等分に離散化し設計を行った。この時、ノッチの数Ｂ＝３として設計された可変フィルタの振幅特性を図７、８に示す。これらはそれぞれ、正弦波の入力信号の周波数が１６８０Ｈｚ、１９００Ｈｚである場合の振幅特性である。そのため、伝送零点の位置はこれらの入力信号の周波数の２、４、６倍の位置となり、これを破線で示す。これより得られた可変フィルタは準等リプル特性が得られていることがわかる。

【0097】

また、可変フィルタのパラメータは瞬時周波数の推定値を用いて変更する必要がある。しかし、フィルタ処理前の瞬時周波数を用いると振動成分の影響で伝送零点位置も変動してしまう。そのため、初期値に上記（５）式により算出されたωを用いて、その後可変フィルタから出力された周波数をパラメータとして用いた。

【0098】

入力信号に対して瞬時周波数を推定した結果を図９に示す。ここではヒルベルト変換器を用いて推定された瞬時周波数、可変フィルタを用いて振動成分を除去した場合の結果を重ねて示している。また、推定結果を拡大したグラフを図１０、１１に示す。これらより、入力信号の周波数によってヒルベルト変換器のリプル、すなわち誤差δが異なるため、推定精度も異なることがわかる。例えば、図１０に示した１７００Ｈｚは図１のω_１であり、そのときδ＝－１２．１８×１０^－３である。また、図１１に示した１９００Ｈｚは図１のω_２であり、そのときδ＝－２．３５×１０^－３である。図１０、１１からも判るように、伝送零点の数は、誤差δの大きさや精度によって決める必要がある。

【0099】

図１２には、本開示の技術を、時間と共に周波数が増加するチャープ信号に適用して周波数を推定した結果を示した。図１３には、図１２の一部を拡大した図を示した。図１２、１３に示すように、ヒルベルト変換のみを行った場合（ＨＴ ｏｎｌｙ）と比較して、Ｂ＝２、Ｂ＝３の場合の何れも理想（ｉｄｅａｌ）の周波数に追従できていることが判った。

【0100】

以上より、本開示の技術によって瞬時周波数に含まれる振動成分を効果的に除去できていることが確認できた。

【0101】

＜変形例＞

【0102】

前述したように、瞬時周波数に生じる振動成分の周波数は、入力信号の周波数の偶数倍となり、この偶数倍の瞬時周波数の振動成分を除去する必要がある。ここで、瞬時周波数が、入力信号をサンプリングするサンプリング周波数に対応するナイキスト周波数よりも大きい周波数の場合、所謂折り返しノイズが発生する。

【0103】

例えば図１４に示すように、入力信号ｘ（ｔ）の周波数ｆ_０＝４ｋＨｚ、サンプリング周波数ｆ_Ｓ＝１０ｋＨｚ、ナイキスト周波数ｆ_Ｎ＝５ｋＨｚであったとする。この場合、周波数ｆ_０の２倍である８ｋＨｚ、４倍である１６ｋＨｚ・・・、の瞬時周波数の振動成分を除去する必要がある。

【0104】

ここで、例えば８ｋＨｚの瞬時周波数については、折り返し周波数としてのナイキスト周波数ｆ_Ｎを対称軸として２ｋＨｚの周波数に折り返しノイズが発生するため、２ｋＨｚの周波数の折り返しノイズも除去する必要がある。

【0105】

そこで、瞬時周波数が、入力信号をサンプリングするサンプリング周波数に対応するナイキスト周波数よりも大きい周波数の場合、フィルタ係数設定部１８は、瞬時周波数及び入力信号の周波数の偶数倍の周波数成分の折り返しノイズが除去されるようにフィルタ係数を設定するようにしてもよい。これにより、精度良く周波数を推定することができる。

【0106】

なお、本実施形態は、請求項に係る発明を限定するものではなく、また実施の形態の中で説明されている特徴の組み合わせの全てが発明の解決手段に必須であるとは限らない。前述した実施の形態には種々の段階の発明が含まれており、開示される複数の構成要件の組み合わせにより種々の発明が抽出される。実施の形態に示される全構成要件から幾つかの構成要件が削除されても、効果が得られる限りにおいて、この幾つかの構成要件が削除された構成が発明として抽出され得る。

【符号の説明】

【0107】

１０ 周波数推定装置
１２ ヒルベルト変換部
１４ 解析信号生成部
１６ 瞬時周波数算出部
１８ フィルタ係数設定部
２０ フィルタ部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】