特許 5712255 フーリエ解析による周波数測定方法および周波数測定装置

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

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

(11)【特許番号】5712255

(24)【登録日】2015年3月13日

(45)【発行日】2015年5月7日

(54)【発明の名称】フーリエ解析による周波数測定方法および周波数測定装置

(51)【国際特許分類】

   G01R 23/16 20060101AFI20150416BHJP

【ＦＩ】

   G01R23/16 D

【請求項の数】4

【全頁数】29

(21)【出願番号】特願2013-167665(P2013-167665)

(22)【出願日】2013年8月12日

(65)【公開番号】特開2015-36635(P2015-36635A)

(43)【公開日】2015年2月23日

【審査請求日】2013年9月6日

(73)【特許権者】

【識別番号】592061599

【氏名又は名称】株式会社近計システム

(74)【代理人】

【識別番号】100084146

【弁理士】

【氏名又は名称】山崎 宏

(74)【代理人】

【識別番号】100081422

【弁理士】

【氏名又は名称】田中 光雄

(74)【代理人】

【識別番号】100118625

【弁理士】

【氏名又は名称】大畠 康

(74)【代理人】

【識別番号】100183232

【弁理士】

【氏名又は名称】山崎 敏行

(72)【発明者】

【氏名】大橋 善和

(72)【発明者】

【氏名】徐 有恒

(72)【発明者】

【氏名】大橋 英也

(72)【発明者】

【氏名】牧村 達也

(72)【発明者】

【氏名】大浦 好文

【審査官】 菅藤 政明

(56)【参考文献】

【文献】 特許第２５０５７０７（ＪＰ，Ｂ２）

【文献】 特表平１０−５０４６４７（ＪＰ，Ａ）

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

Ｇ０１Ｒ ２３／１６

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

【請求項1】

交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに対して窓関数演算部により窓関数を掛ける第１ステップと、
上記第１ステップで上記窓関数を掛けた上記サンプリング波形データに対してフーリエ変換部によりフーリエ変換を行う第２ステップと、
上記第２ステップで得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを(式１)に代入して周波数演算部により演算することによって、上記被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yを求める第３ステップと、
上記第３ステップで得られた上記被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yの有効性を(式２)に基づいて判定部により判定する第４ステップと
を有することを特徴とするフーリエ解析による周波数測定方法。

【数1】

【数2】

ただし、 Ｖ_k ：時間幅Ｔ[秒]の被測定信号のＶ(t)をフーリエ変換して得られる
ｋ／Ｔ[Ｈｚ］の複素周波数成分(ここでｋ＝１,２,…)
Ｖ_m ：上記複素周波数成分のうち絶対値が最大のもの
ｎ ：窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
Ｆ_X ：被測定信号の周波数
Ｆ_Y ：被測定信号の減衰発散係数
Ｔ ：被測定信号のサンプリング波形データの演算対象部分の時間幅
Ｖ_min：(式１)の演算結果の有効性を判定するための閾値

【請求項2】

交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに窓関数を掛ける窓関数演算部と、
上記窓関数演算部で上記窓関数を掛けた上記サンプリング波形データに対してフーリエ変換を行うフーリエ変換部と、
上記フーリエ変換部で得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを(式１)に代入して演算することによって、上記被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yを求める周波数演算部と、
上記周波数演算部で上記(式１)により得られた上記被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yの有効性を(式２)に基づいて判定する判定部と
を備えたことを特徴とする周波数測定装置。

【数3】

【数4】

ただし、 Ｖ_k ：時間幅Ｔ[秒]の被測定信号のＶ(t)をフーリエ変換して得られる
ｋ／Ｔ[Ｈｚ］の複素周波数成分(ここでｋ＝１,２,…)
Ｖ_m ：上記複素周波数成分のうち絶対値が最大のもの
ｎ ：窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
Ｆ_X ：被測定信号の周波数
Ｆ_Y ：被測定信号の減衰発散係数
Ｔ ：被測定信号のサンプリング波形データの演算対象部分の時間幅
Ｖ_min：(式１)の演算結果の有効性を判定するための閾値

【請求項3】

交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに対してフーリエ変換部によりフーリエ変換を行う第１ステップと、
上記第１ステップで得られた複素周波数成分のデータに対して窓関数演算部により(式４)による演算を行って、上記第１ステップの前に上記被測定信号のサンプリング波形データに(式３)で表された窓関数を掛けた演算結果と同等の演算結果を得る第２ステップと、
上記第２ステップで得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを(式１)に代入して周波数演算部により演算することによって、上記被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yを求める第３ステップと、
上記第３ステップで得られた上記被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yの有効性を(式２)に基づいて判定部により判定する第４ステップと
を有することを特徴とするフーリエ解析による周波数測定方法。

【数5】

【数6】

ただし、 Ｖ_k ：時間幅Ｔ[秒]の被測定信号のＶ(t)をフーリエ変換して得られる
ｋ／Ｔ[Ｈｚ］の複素周波数成分(ここでｋ＝１,２,…)
Ｖ_m ：上記複素周波数成分のうち絶対値が最大のもの
ｎ ：窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
Ｆ_X ：被測定信号の周波数
Ｆ_Y ：被測定信号の減衰発散係数
Ｔ ：被測定信号のサンプリング波形データの演算対象部分の時間幅
Ｖ_min：(式１)の演算結果の有効性を判定するための閾値

【数7】

【数8】

ただし、 w(t)：窓関数
α：窓関数の直流成分の係数
β、γ、δ：窓関数の各周波数成分の係数
C_k ：窓関数を掛けない状態でフーリエ変換した結果の
ｋ／Ｔ[Ｈz]成分のフーリエ係数
C’_k ：上記窓関数w(t)を掛けてフーリエ変換した結果の
ｋ／Ｔ[Ｈz]成分のフーリエ係数

【請求項4】

交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに対してフーリエ変換を行うフーリエ変換部と、
上記フーリエ変換部で得られた複素周波数成分のデータに対して(式４)による演算を行って、上記フーリエ変換前に上記被測定信号のサンプリング波形データに(式３)で表された窓関数を掛けた演算結果と同等の演算結果を得る窓関数演算部と、
上記窓関数演算部で得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを(式１)に代入して演算することによって、上記被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yを求める周波数演算部と、
上記周波数演算部で上記(式１)により得られた上記被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yの有効性を(式２)に基づいて判定する判定部と
を備えたことを特徴とする周波数測定装置。

【数9】

【数10】

ただし、 Ｖ_k ：時間幅Ｔ[秒]の被測定信号のＶ(t)をフーリエ変換して得られる
ｋ／Ｔ[Ｈｚ］の複素周波数成分(ここでｋ＝１,２,…)
Ｖ_m ：上記複素周波数成分のうち絶対値が最大のもの
ｎ ：窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
Ｆ_X ：被測定信号の周波数
Ｆ_Y ：被測定信号の減衰発散係数
Ｔ ：被測定信号のサンプリング波形データの演算対象部分の時間幅
Ｖ_min：(式１)の演算結果の有効性を判定するための閾値

【数11】

【数12】

ただし、 w(t)：窓関数
α：窓関数の直流成分の係数
β、γ、δ：窓関数の各周波数成分の係数
C_k ：窓関数を掛けない状態でフーリエ変換した結果の
ｋ／Ｔ[Ｈz]成分のフーリエ係数
C’_k ：上記窓関数w(t)を掛けてフーリエ変換した結果の
ｋ／Ｔ[Ｈz]成分のフーリエ係数

【発明の詳細な説明】

【技術分野】

【0001】

この発明は、フーリエ解析による周波数測定方法および周波数測定装置に関する。

【背景技術】

【0002】

従来、交流波形の周波数を測定する第一の周波数測定方法としては、波形の値が零となる零クロス点を見つけて零クロス点から次の零クロス点までの時間幅を図形的に求めてその時間幅を周期とみなしてその逆数として周波数を求めるものがある。

【0003】

上記第一の周波数測定方法では、零クロス点周辺にノイズ成分が含まれている場合、正しく周期測定ができない。特に、インバーター機器などの最近のパワーエレクトロニクス技術の普及の結果、電源波形に含まれるノイズ成分は増大する方向にある。

【0004】

そこで、このようなノイズ成分の影響を受けない従来の第二の周波数測定方法として、波形データの周波数スペクトルから演算で周波数を求めるものがある(例えば、特許第２５０５７０７号(特許文献１)参照)。

【0005】

上記第二の周波数測定方法では、周波数スペクトルを算出するのにＦＦＴ演算を使うが、リアルタイム性に欠けており、そのままの形ではオンラインでの使用はできなかった。また、周波数の測定精度を上げるためには、交流波形の数サイクルに渡る波形に窓関数を掛けてＦＦＴ演算をしなければならないが、演算結果を波形の急激な変化に追従させるには無理があった。

【0006】

一般に時間幅Ｔのデータをフーリエ解析(ＦＦＴ(高速フーリエ変換)またはＤＦＴ(離散フーリエ変換))して周波数成分分析を行なった場合、もとの波形に含まれる振動成分が商用周波数成分一つだけの場合でもそれがＦＦＴによって得られる各スペクトルの周波数(ＦＦＴのＢＩＮ周波数という)からずれていると多数のスペクトルに分解される。

【0007】

例えば、１００ｍｓｅｃのデータをＦＦＴ分析して１０Ｈz毎の周波数スペクトルを算出した場合、商用周波数成分が５１Ｈzとしたときに得られる結果は、５０Ｈz周辺に出現する多数のスペクトルで表現されることになる。

【0008】

そもそもフーリエ変換の解析結果は、そこに各周波数成分が有ることを保証するものではなく、解析結果の周波数成分の波形の和がそのデータの区間で元の波形にもっとも近い(二乗平均誤差が最小になる)ことを示しているだけである。

【0009】

上記特許文献１には、一個の周波数成分が複数のスペクトルに分散した場合、元の波形の周波数が幾らだったのかを算出する方法が書かれている。

【0010】

この特許文献１は、被測定信号の周波数Ｆ_Xを(式５)で算出するもので、減衰発散係数Ｆ_Yも得ることができる。

【0011】

【数1】

ただし、 Ｆ_X ：被測定信号の周波数
Ｆ_Y ：被測定信号の減衰発散係数
Ｖ₁ 〜Ｖ_n ：被測定信号の複素周波数成分
Ｆ₁ 〜Ｆ_n ：複素周波数成分Ｖ₁ 〜Ｖ_n の各周波数

【0012】

ここで、波形を表す複素関数ｆ(t)は、実用上は実数関数であるが、時間ｔの関数であり、次の(式６)のように複素フーリエ級数で表現できるものである。

【0013】

【数2】

【0014】

ただし、Ｖ_kは第ｋ番目の複素周波数成分、Ｔ[ｓｅｃ]は観測データが存する時間幅である。周波数スペクトルの間隔を△ｆ[Ｈz]とすると、△ｆ＝１／Ｔの関係がある。また、第ｋ番目のスペクトルの周波数Ｆ_k[Ｈz]は、
Ｆ_k＝ｋ△ｆ＝ｋ／Ｔ
であり、角周波数ω_k[ｒａｄ／ｓｅｃ]とは、
ω_k＝ｋ△ω＝２πｋ△ｆ
の関係がある。

【0015】

この第二の周波数測定方法によれば、確かに特定時間断面における周波数成分の中心周波数を荒いサンプリング波形データ(例えば３２サンプル／交流の１サイクル程度)でも精度よく求めることができる。

【先行技術文献】

【特許文献】

【0016】

【特許文献1】特許第２５０５７０７号

【発明の概要】

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

【0017】

しかしながら、上記第二の周波数測定方法(特許文献１)では、特定の時間断面での値だけではなく、周波数の変化をリアルタイムで観測したい場合、一定の時間間隔(例えば２ｍｓｅｃ毎)で何度もこの手法を適用してＦＦＴを繰り返さねばならず、処理時間が掛かるためにリアルタイム性を持たせることには難点があった。

【0018】

また、上記第二の周波数測定方法を装置に実装するには、ＦＦＴプログラムを組み込む必要があり、種類豊富な演算専用チップを使用する場合、チップ専用のマシン語を習得して直接プログラミングすることは困難なので、どうしてもＣ言語のようなある程度規格化された高級言語を使わざるを得なかったが、高速性が要求されるＦＦＴ演算に高級言語を用いることは、高速性が犠牲になるので好ましいことではない。

【0019】

そこで、この発明の課題は、交流波形の２サイクル〜４サイクル程度で、かつ、交流波形の１サイクルあたり数十サンプルといった比較的低いサンプリング周波数によりサンプリングされたサンプリング波形データからその交流波形の周波数を高速かつ高精度に測定できるフーリエ解析による周波数測定方法および周波数測定装置を提供することにある。

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

【0020】

上記課題を解決するため、第１の発明のフーリエ解析による周波数測定方法は、
交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに対して窓関数演算部により窓関数を掛ける第１ステップと、
上記第１ステップで上記窓関数を掛けた上記サンプリング波形データに対してフーリエ変換部によりフーリエ変換を行う第２ステップと、
上記第２ステップで得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを(式１)に代入して周波数演算部により演算することによって、上記被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yを求める第３ステップと、
上記第３ステップで得られた上記被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yの有効性を(式２)に基づいて判定部により判定する第４ステップと
を有することを特徴とする。

【0021】

ただし、窓関数演算部については後述の方法により、フーリエ変換を先に実施して、その結果得られる複素スペクトル間の簡単な演算で最初に窓関数を掛けたと同じ効果を奏する手法を適用しても良い。

【0022】

【数3】

【0023】

【数4】

ただし、 Ｖ_k ：時間幅Ｔ[秒]の被測定信号のＶ(t)をフーリエ変換して得られる
ｋ／Ｔ[Ｈｚ］の複素周波数成分(ここでｋ＝１,２,…)
Ｖ_m ：上記複素周波数成分のうち絶対値が最大のもの
ｎ ：窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
Ｆ_X ：被測定信号の周波数
Ｆ_Y ：被測定信号の減衰発散係数
Ｔ ：被測定信号のサンプリング波形データの演算対象部分の時間幅
Ｖ_min：(式１)の演算結果の有効性を判定するための閾値

【0024】

上記構成によれば、交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに窓関数を掛けた後、窓関数を掛けたサンプリング波形データに対してフーリエ変換を行って、そのフーリエ変換により得られた複素周波数スペクトル値の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接複素周波数成分とを(式１)に代入して、上記被測定信号の周波数と減衰発散係数を求める。また、上記被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yの有効性を(式２)に基づいて判定して、(式２)の条件を満たさないなら周波数Ｆ_Xと減衰発散係数Ｆ_Yの演算結果を無効にすることで、量子化誤差等の影響によって不安定になった計算結果を除外できる。これにより、交流波形の２サイクル〜４サイクル程度で、かつ、交流波形の１サイクルあたり数十サンプルといった比較的低いサンプリング周波数によりサンプリングされたサンプリング波形データからその交流波形の周波数を高速かつ高精度に測定できる。特に、窓関数に含まれる周波数成分の個数相当分だけのスペクトルを演算に用いているので、不要な周波数成分を算出する必要がなく、フーリエ変換に際してＦＦＴを用いて全周波数成分(例えば２５６サンプルデータなら１２８個の複素スペクトルデータ）を算出する必要がなく、ＤＦＴを用いて必要な周波数成分(窓関数内の周波数成分が１個なら中心周波数成分のスペクトル１個と低周波数側および高周波数側にそれぞれ１個〜２個、全部で３個〜５個のスペクトルで十分)のみを算出することで演算速度を遥かに高速化することができたのである。

【0025】

また、第２の発明の周波数測定装置は、
交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに窓関数を掛ける窓関数演算部と、
上記窓関数演算部で上記窓関数を掛けた上記サンプリング波形データに対してフーリエ変換を行うフーリエ変換部と、
上記フーリエ変換部で得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを(式１)に代入して演算することによって、上記被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yを求める周波数演算部と、
上記周波数演算部で上記(式１)により得られた上記被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yの有効性を(式２)に基づいて判定する判定部と
を備えたことを特徴とする。

【0026】

【数5】

【0027】

【数6】

ただし、 Ｖ_k ：時間幅Ｔ[秒]の被測定信号のＶ(t)をフーリエ変換して得られる
ｋ／Ｔ[Ｈｚ］の複素周波数成分(ここでｋ＝１,２,…)
Ｖ_m ：上記複素周波数成分のうち絶対値が最大のもの
ｎ ：窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
Ｆ_X ：被測定信号の周波数
Ｆ_Y ：被測定信号の減衰発散係数
Ｔ ：被測定信号のサンプリング波形データの演算対象部分の時間幅
Ｖ_min：(式１)の演算結果の有効性を判定するための閾値

【0028】

上記構成によれば、第１の発明と同様の効果を奏する。

【0029】

また、第３の発明のフーリエ解析による周波数測定方法は、
交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに対してフーリエ変換部によりフーリエ変換を行う第１ステップと、
上記第１ステップで得られた複素周波数成分のデータに対して窓関数効果演算部により(式４)による演算を行って、上記第１ステップの前に上記被測定信号のサンプリング波形データに(式３)で表された窓関数を掛けた演算結果と同等の演算結果を得る第２ステップと、
上記第２ステップで得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを(式１)に代入して周波数演算部により演算することによって、上記被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yを求める第３ステップと、
上記第３ステップで得られた上記被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yの有効性を(式２)に基づいて判定部により判定する第４ステップと
を有することを特徴とする。

【0030】

【数7】

【0031】

【数8】

ただし、 Ｖ_k ：時間幅Ｔ[秒]の被測定信号のＶ(t)をフーリエ変換して得られる
ｋ／Ｔ[Ｈｚ］の複素周波数成分(ここでｋ＝１,２,…)
Ｖ_m ：上記複素周波数成分のうち絶対値が最大のもの
ｎ ：窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
Ｆ_X ：被測定信号の周波数
Ｆ_Y ：被測定信号の減衰発散係数
Ｔ ：被測定信号のサンプリング波形データの演算対象部分の時間幅
Ｖ_min：(式１)の演算結果の有効性を判定するための閾値

【0032】

【数9】

【0033】

【数10】

ただし、 w(t)：窓関数
α：窓関数の直流成分の係数
β、γ、δ：窓関数の各周波数成分の係数
C_k ：窓関数を掛けない状態でフーリエ変換した結果の
ｋ／Ｔ[Ｈz]成分のフーリエ係数
C’_k ：上記窓関数w(t)を掛けてフーリエ変換した結果の
ｋ／Ｔ[Ｈz]成分のフーリエ係数

【0034】

上記構成によれば、交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに対してフーリエ変換を行った後、フーリエ変換で得られた複素周波数成分のデータに対して窓関数効果演算部により(式４)による演算を行って、フーリエ変換の前に被測定信号のサンプリング波形データに(式５)で表された窓関数を掛けた演算結果と同等の演算結果を得るので、第１の発明と同様の効果を奏する。また、最初に窓関数を掛ける必要がないので、１サンプル毎に周波数を算出する必要がある場合に再帰的ＤＦＴを利用することによって、演算速度を飛躍的に改善できる。

【0035】

また、第４の発明の周波数測定装置は、
交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに対してフーリエ変換を行うフーリエ変換部と、
上記フーリエ変換部で得られた複素周波数成分のデータに対して(式４)による演算を行って、上記フーリエ変換前に上記被測定信号のサンプリング波形データに(式３)で表された窓関数を掛けた演算結果と同等の演算結果を得る窓関数効果演算部と、
上記窓関数効果演算部で得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを(式１)に代入して演算することによって、上記被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yを求める周波数演算部と、
上記周波数演算部で上記(式１)により得られた上記被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yの有効性を(式２)に基づいて判定する判定部と
を備えたことを特徴とする。

【0036】

【数11】

【0037】

【数12】

ただし、 Ｖ_k ：時間幅Ｔ[秒]の被測定信号のＶ(t)をフーリエ変換して得られる
ｋ／Ｔ[Ｈｚ］の複素周波数成分(ここでｋ＝１,２,…)
Ｖ_m ：上記複素周波数成分のうち絶対値が最大のもの
ｎ ：窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
Ｆ_X ：被測定信号の周波数
Ｆ_Y ：被測定信号の減衰発散係数
Ｔ ：被測定信号のサンプリング波形データの演算対象部分の時間幅
Ｖ_min：(式１)の演算結果の有効性を判定するための閾値

【0038】

【数13】

【0039】

【数14】

ただし、 w(t)：窓関数
α：窓関数の直流成分の係数
β、γ、δ：窓関数の各周波数成分の係数
C_k ：窓関数を掛けない状態でフーリエ変換した結果の
ｋ／Ｔ[Ｈz]成分のフーリエ係数
C’_k ：上記窓関数w(t)を掛けてフーリエ変換した結果の
ｋ／Ｔ[Ｈz]成分のフーリエ係数

【0040】

上記構成によれば、第３の発明と同様の効果を奏する。

【発明の効果】

【0041】

以上より明らかなように、この発明によれば、交流電圧波形または交流電流波形を表す被測定信号の２サイクル〜４サイクル程度で、かつ、交流波形の１サイクルあたり数十サンプルといった比較的低いサンプリング周波数によりサンプリングされたサンプリング波形データからでも、被測定信号の周波数を高速かつ高精度に測定することができるフーリエ解析による周波数測定方法および周波数測定装置を実現することができる。

【図面の簡単な説明】

【0042】

【図1】図１はこの発明の第１実施形態のフーリエ解析による周波数測定方法を用いた周波数測定装置の要部ブロック図である。

【図2】図２はこの発明の第２実施形態のフーリエ解析による周波数測定方法を用いた周波数測定装置の要部ブロック図である。

【図3】図３は上記周波数測定装置の周波数演算部の構成を示す図である。

【図4】図４は上記周波数測定装置により測定された周波数の変化の検出例を示す図である。

【図5】図５は図４に示す周波数変化の要部の拡大図である。

【図6】図６は上記第１実施形態の周波数測定装置と上記第２実施形態の周波数測定装置の周波数測定結果の比較例を示す図である。

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

【0043】

この発明のフーリエ解析による周波数測定方法および周波数測定装置について、まず、特許第２５０５７０７号(特許文献１)に記載された数式(１)であって本明細書における(式５)の算出根拠について説明し、その後に上述の課題を解決する手段および実施形態の詳細を説明する。

【0044】

いま、交流波形を表す複素関数ｆ(t)を、

とし、複素関数ｆ(t)が一個の周波数成分として、

を持っていたとする。このときの未知の角周波数成分ω_xが、
ω_x＝２πＦ_k または Ｆ_k＝ω_x／２π
のように解析周波数Ｆ_kのいずれかと一致している場合は、その成分のみが検出され、他の成分は零なので(式５)が成り立つのは明らかである。

【0045】

そうでない場合、複素関数ｆ(t)の解析結果は、(式５)のように複数の複素周波数成分(スペクトル成分)の和で表されるが、無限個の計算はできないので、何らかの窓関数を掛けることで第k1個目から第k2個目(k2＞k1)までのスペクトルを除く他のスペクトルを十分小さくし、第k1個目から第k2個目までのスペクトルで充分な精度でかつ角周波数ω_xが複素数の範囲において、次の(式７)が成立するものとする。

【0046】

【数15】

【0047】

この(式７)は時間ｔが−ｔでも成立するので、次の(式８)が得られる。

【0048】

【数16】

【0049】

上記(式７)と(式８)の辺々を掛け合わせると、次の(式９)となる。ここで、(式９)の右辺の第二項のΣの記号の下部の式は、「ｋおよびｌが異なる整数値で、かつ、k1以上k2以下となる全ての組み合わせの総和」を示すものとする。

【0050】

【数17】

【0051】

上記(式９)の左辺のＡ^２および右辺の第一項が時間ｔに関係しない定数項であるので、右辺の第二項も定数であることが判る。また、△ω＝２π／Ｔのため、区間ｔ＝[０、Ｔ]で積分すると、

で表される複素数ベクトルは、原点の周りを(k-l)周(ただし、ｋ≠ｌ)するので、そのベクトル和は零になり、区間ｔ＝[０、Ｔ]での時間積分値が零であることが判る。

【0052】

すなわち、(式９)の両辺を時刻ｔ、観測区間[０、Ｔ]で積分すると、次の(式１０)および(式１１)で表現できるが、(式１１)の右辺の第二項が零になることが判る。

【0053】

【数18】

【0054】

【数19】

【0055】

それ故、(式９)の右辺の第二項は、時間ｔに関係ない定数であって、その時間積分値が零であることから、もともと常に零であることが判る。よって、(式９)は、

【0056】

【数20】

となる。

【0057】

一方、(式７)の両辺をｔで微分することにより、次の(式１２)が得られる。

【0058】

【数21】

【0059】

この(式１２)に(式８)の両辺を掛け合わせると、(式１３)が得られる。

【0060】

【数22】

【0061】

この(式１３)の右辺の第二項も(式９)と同様の理由で零になる。

【0062】

また、左辺のA^２は各周波数成分の二乗和なので、以下の(式１４)が得られる。

【0063】

【数23】

【0064】

ただし、波形自体が０値で無入力状態であれば当然のことながら、その周波数を求めることは意味がない。よって、波形の実効値がある値以下なら(式１４)での計算結果も無効であると判断する機能も必要である。一般に歪波の実効値は、各周波数成分の実効値の二乗和の平方根で計算されるが、(式１４)の分母は丁度その実効値ベクトルの二乗値になっている。それ故、この点については、(式１４)の分母の絶対値がＡ／Ｄ変換器の分解能に比して十分に大きな値であることを保証することで足る。

【0065】

一方、(式１４)と同様の考え方で(式８)の両辺を掛け合わせる代わりに(式１２)の辺々を二乗して整理することで次の(式１５)も得られる。

【0066】

【数24】

【0067】

(式１５)でも同様の結果が得られるが、この場合は得られるのがω_xの二乗値であり、その平方根は２個出てくるので、別途判別する必要がある。(式１４)の方が簡単で使い勝手が良いので、通常は(式１４)または角周波数の代わりに単に周波数で表現した(式５)を使う。等号を用いているが、(式５)や(式１４)で等号が成り立つのは、第k1番目〜第k2番目(k1＜k2)以外の周波数成分が十分小さく無視できる場合に限る。

【0068】

一般に、交流波形の周波数がＢＩＮ周波数からずれていると、そのスペクトラムは多数の周波数成分に散逸するが、窓関数を掛けるとそれを一定の範囲内に抑えることができる。窓関数は、４〜５桁程度の精度が必要な場合は４サイクル分のサンプリング波形データにハニング窓を掛け、６〜８桁の精度が必要な場合は８サイクル分のサンプリング波形データにナットール窓やブラックマン・ハリス窓を掛ける。

【0069】

以上が(式５)および(式１４)の導出の根拠であるが、この計算をＰＬＤ(プログラマブルロジックデバイス)などで実現しようとすると、各周波数成分を求めるＦＦＴ(高速フーリエ変換)演算などが複雑で実施困難である。ＤＳＰ(Digital Signal Processor：デジタル・シグナル・プロセッサ)のような演算専用チップを用いると(式５),(式１４)の計算も実施可能であるが、高速性を重視したマシン語によるプログラミングは困難である。このような演算専用チップを用いた計算で高速性を出すためには、例えばＣ言語のような高級言語によって記述されたソフトウェアをコンパイルして使用するのは効率的ではなく、アセンブリ言語やマシン語で記述するのが適切である。しかしながら、様々な演算用チップが市販されているので、それらを統一的に同じソースコードでプログラミングすることは不可能で、個別にソフトウェアを開発し、多大な時間を掛けてデバッグしているのが現状であった。

【0070】

＜フーリエ変換方式について＞
商用周波数の交流波形からその周波数を求める場合は、ＦＦＴ解析結果のスペクトルの周波数の一つが５０Ｈzまたは６０Ｈzに一致するように、サンプリング周波数とＦＦＴ演算を行うデータの時間幅(データの個数)とを選定し、商用周波数に近くかつ最も大きな複素周波数成分のスペクトルとその前後に１本か２本のスペクトルを選んでそれらの複素数値を後述する(式１)に代入すれば良いので、敢えてＦＦＴ演算を行ってすべての周波数成分を求める必要はなく、ＤＦＴ(離散フーリエ変換)演算で十分である。

【0071】

例えば、５０Ｈzの交流波形を６４サンプル／交流１サイクルのサンプリングで２５６個のサンプリング波形データを得て、それにＦＦＴ演算を施すと１２８本のスペクトルが得られるが、この場合のスペクトル間の周波数間隔は５０×６４／２５６＝１２．５Ｈzである。

【0072】

ハニング窓関数を用いた場合、交流波形ｖ(t)は以下の(式１６)のように表現できるが、これは(式１７)のように３つの周波数成分に分解できる。

【0073】

【数25】

ただし、ｆは交流波形の周波数、Ｔは観測波形データの時間幅である。

【0074】

したがって、この３本のスペクトルを表す複素数値を(式５)に代入すれば、商用周波の交流波形から周波数が算出できる。万一、それが本来の周波数(例えば５０Ｈz)からずれていれば、この３本のスペクトルの外側にも零でないスペクトルが現れる。それが、(式５)の演算対象から外れると誤差が発生する。その場合、中心周波数の両側に１本ではなく２本または３本のスペクトルを演算対象とすると誤差は小さくなる。

【0075】

一般的に窓関数の例としては、以下の(式１８)の形をしており、定数成分α以外にＭ個(Ｍ＝１〜３程度の整数)の周波数成分を持っている。そのため、それが正弦波交流波形の式と掛け合わされると、その個数Ｍ×２個の側帯波を生ずることになる。それ故、それら全てを(式５)の演算対象にすれば精度上問題のない結果が得られることが判る。当然のことながら、測定対象の波形は周波数変動するのであるから、その変動幅も考慮する必要がある。

【0076】

[窓関数の例]

【0077】

【数26】

【0078】

[係数値の具体例]
・ハニング窓の場合

【0079】

【数27】

【0080】

・ブラックマン・ハリス窓の場合

【0081】

【数28】

【0082】

・ナットール窓の場合

【0083】

【数29】

【0084】

これらの式から明らかなように、例えばハニング窓では周波数成分は直流分αを含めて２個であるので、これを用いて交流波形の周波数成分を分析する場合は、中心周波数のスペクトル(中心複素周波数成分)と、その中心周波数のスペクトルによりも低周波数側の１本のスペクトル(隣接する低周波数側の複素周波数成分)と、高周波数側の１本のスペクトル(隣接する高周波数側の複素周波数成分)の計３本のスペクトルが得られれば計算可能である。また、ブラックマン・ハリス窓やナットール窓では計７本である。

【0085】

しかるに、スペクトルの本数が増えるとＤＦＴの演算量が増えるので、計算速度が遅くなってしまう。また、スペクトル間の周波数間隔を狭めようとするとデータの時間幅を大きくしなければならず、そうなると周波数の急変に追従させることが難しくなる。そのため、本発明者は、追従性が要求される用途では、スペクトルの本数を極力絞って計算速度と追従性能を改善し、また、追従性能を要求されない用途では、窓関数を考慮したＤＦＴの計算部分をあらかじめ実行しておき、実時間ではより簡単な計算式で周波数を算出する方法を考えた。すなわち、高速性・高追従性を有する方法と、高速性・高精度を実現できる方法とを考案し検証した。

【0086】

次に、スペクトルの本数を極力絞って計算速度と追従性能を改善する方法について説明する。

【0087】

サンプリング波形データの時間幅を交流波形の本来の周期の２倍にとり、交流波形２サイクル分のサンプリング波形データについて３本のスペクトルを計算する。交流波形の本来の周波数が５０Ｈzの場合は、２５Ｈz、５０Ｈz、７５Ｈzの３本のスペクトルをＤＦＴの計算式(式１９)によって計算すると、次の(式２０)に示す計算結果となる。

【0088】

【数30】

ここで、ｋ＝１,２,３である。

【0089】

【数31】

【0090】

しかし、直流成分や２次高調波成分などが含まれていると、それが(式２０)での計算結果に誤差として含まれ、それを分離することはできない。

【0091】

そのような場合には、サンプリング波形データを交流波形の周期の４倍から８倍に渡って収集し、中心周波数の前後の２〜４本、計５〜９本のスペクトルについて(式２０)の計算を行えばよい。

【0092】

ただし、そうすると、スペクトルの本数やデータ数が共に２〜４倍になるので、演算に必要な時間は４〜１６倍になってしまい、演算速度の低下が避けられない。

【0093】

そのような場合には、係数計算を先に済ませておいて、メモリーに保存しておくことで高速演算を可能にする方法が必要となる。

【0094】

このように、本発明のフーリエ解析による周波数測定方法は、第１ステップで交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに窓関数を掛けた後、第２ステップで窓関数を掛けたサンプリング波形データに対してフーリエ変換を行って、第３ステップでそのフーリエ変換により得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを次の(式１)に代入して、上記被測定信号の周波数と減衰発散係数を求める。さらに、第４ステップで被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yの有効性を次の(式２)によって判定して、(式２)の条件を満たさないなら周波数Ｆ_Xと減衰発散係数Ｆ_Yの演算結果を無効にすることで、量子化誤差等の影響によって不安定になった計算結果を除外できる。これにより、交流波形の２サイクル〜４サイクル程度で、かつ、交流波形の１サイクルあたり数十サンプルといった比較的低いサンプリング周波数によりサンプリングされたサンプリング波形データからその交流波形の周波数を高速かつ高精度に測定することができる。

【0095】

【数32】

【0096】

【数33】

ただし、 Ｖ_k ：時間幅Ｔ[秒]の被測定信号のＶ(t)をフーリエ変換して得られる
ｋ／Ｔ[Ｈｚ］の複素周波数成分(ここでｋ＝１,２,…)
Ｖ_m ：上記複素周波数成分のうち絶対値が最大のもの
ｎ ：窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
Ｆ_X ：被測定信号の周波数
Ｆ_Y ：被測定信号の減衰発散係数
Ｔ ：被測定信号のサンプリング波形データの演算対象部分の時間幅
Ｖ_min：(式１)の演算結果の有効性を判定するための閾値

【0097】

ところで、窓関数を掛けることは、周波数スペクトルを一定の周波数範囲に収め、本発明における周波数の演算式における演算項数を限定するうえで必要なステップであるが、先に窓関数を掛けてその後にフーリエ変換することに限定する必要はなく、先にフーリエ変換を行ってその結果に対して例えば(式４）による演算を行うことでも同様の効果を得ることができる。

【0098】

その根拠について以下に説明する。例えば、角周波数ω₀が、

で表され、サンプリング波形データの時間幅Ｔを商用周波の周期の４倍にとって窓関数ｗ(t)を掛けた場合、フーリエ係数Ｃ'_kは以下の(式２１)のようになる。ここで、ｆ(t)は交流波形を表す複素関数、ｗ(t_n)は窓関数、Ｎはデータ個数である。

【0099】

【数34】

【0100】

窓関数ｗ(t_n)としてハニング窓を用いた場合、それを代入すると以下の(式２２)のように変形できる。

【0101】

【数35】

【0102】

これは、ハニング窓を掛けた場合のＦＦＴの結果のフーリエ係数C’_kと、窓関数を掛けない場合のＦＦＴの結果のフーリエ係数C_kとの間に以下の(式２３)の関係があることを示している。

【0103】

【数36】

【0104】

同様に、(式３)の窓関数ｗ(t)に対して(式４)が成立する。これらの結果は、必ずしも先に窓関数を掛けてからフーリエ変換する必要がないことを示しており、使用する窓関数をフーリエ変換後に選択することも可能であることや複数の窓関数の効果の違いを調べる場合にその度毎にフーリエ変換し直す必要がないことを示している。また、最初に窓関数を掛ける必要がなければ再帰的ＤＦＴが利用でき、１サンプル毎に周波数を算出する必要がある場合は演算速度が飛躍的に改善される。

【0105】

【数37】

【0106】

【数38】

【0107】

よって、本発明では、窓関数を掛ける順序(フーリエ変換の前または後)を限定しないこととする。

【0108】

＜フィルタリング方式について＞
一方、本発明者は、上記(式１)が観測値ｆ(t)を用いて例えば(式２４)のような形になれば、簡単な論理デバイスでも実現可能と考えて理論式を検討し、各係数を算出して図３に示すブロック構成で周波数演算機能を実現した。

【0109】

【数39】

【0110】

図３の構成は、入力部の二乗演算と出力部の割り算以外は、一般のフィルターなどで良く使用される演算回路であるので、それらのファームウエアがそのまま利用できる。この手法も結局は(式１)を変形したものであり、原理は同じであるが、演算の手順が異なるので周波数変化時の過渡特性は異なっている。

【0111】

上記(式１)を用いるには、サンプリング波形データからＤＦＴやＦＦＴなどの手法で複素周波数成分を求める手順が必要であり、元の波形データには窓関数を掛けることが必要であった。しかし、これでは毎回窓関数を掛けてＦＦＴをする必要が有り、演算時間が掛かってリアルタイム性が悪いので、本発明者は、(式１)を時間関数上に展開したのである。

【0112】

複素フーリエ級数展開の式から(式１)のＶ_kは、次の(式２５)〜(式２８)のように表現できる。

【0113】

【数40】

【0114】

例えば、６０Ｈz±１０Ｈz程度の商用周波数の波形データをサンプリング周波数１９２０Ｈz(３２サンプル／サイクル)で８サイクル分サンプリングし、２５６個のサンプリング波形データからＤＦＴやＦＦＴによって６０／８＝７.５Ｈz毎に得られる周波数成分スペクトルのうちの６０Ｈzを中心とする１本と、その前後のそれぞれ３本のスペクトル(計７本のスペクトル)から(式５)で周波数を求めるものとする。

【0115】

この場合、データ個数はＮ＝２５６、サンプリング波形データの時間幅はＴ＝０．１３３３…[sec]であり、各スペクトルは７.５Ｈz、１５Ｈz、２２.５Ｈz、３０Ｈz、３７.５Ｈz、４５Ｈz、５２.５Ｈz、６０Ｈz、６７.５Ｈz、７５Ｈz、８２.５Ｈz、９０Ｈz、…が得られる。このうちの３７.５Ｈz〜８２.５Ｈz(ｋ＝５〜１１)の周波数成分スペクトルを用いて基本波の周波数を算出する。

【0116】

ところで、次数の異なる高調波の周波数成分の和の実効値の二乗は、各次数の周波数成分の実効値の二乗和に等しいので、(式１)の各周波数成分は以下の(式２９)で表現できる。この(式２９)の右辺の括弧内は倍角公式を用いて(式３０)のように変形できる。

【0117】

【数41】

【0118】

よって、(式１)の分子は(式３１)で表現できると共に、

【0119】

【数42】

【0120】

(式１)の分母は(式３２)で表現できることが判る。

【0121】

【数43】

【0122】

上記(式３１)および(式３２)でΣ記号による加算演算の順序を入れ替えている点は重要である。この入れ替えによって、これまで各周波数におけるＤＦＴ演算を各波形データに対して行ってその結果を加算していたものが、各周波数におけるＤＦＴ演算の係数を先に加算して置いておき、後で波形データに掛け合わせることができるので、毎回時間のかかる演算をする必要がなく、実際のオンライン演算での演算時間を短縮することができる。

【0123】

ここで、(式１)の分子の式の係数をＡ_nで表し、分母の式の係数をＢ_nで表すと、(式３３)および(式３４)が得られる。

【0124】

【数44】

【0125】

【数45】

【0126】

この(式３３)、(式３４)の値は前もって計算しておくことのできる値である。実際の計算では、これに更に窓関数Ｗ(２πｎ／Ｎ)を掛けて、(式３５)や(式３６)の値を計算しておく。

【0127】

【数46】

【0128】

【数47】

ただし、Ｗは窓関数であり、例えばハニング窓の場合は(式３７)のようになるが、これ以外の窓関数を用いてもよい。

【0129】

【数48】

【0130】

求める周波数Ｆ_Xと減衰発散係数Ｆ_Yは、上記(式３５),(式３６)の係数を用いて、次の(式３８)で算出できる。上記(式３５),(式３６)は、入力される波形データによって変化する値ではないので、最初に一回係数Ａ_nと係数Ｂ_nを計算するだけで済む。その計算された係数Ａ_nと係数Ｂ_nをメモリー部に記憶して次の計算で読み出して用いる。なお、上記(式３５),(式３６)の計算結果を不揮発性メモリーに保存しておいてもよい。

【0131】

【数49】

ただし、 Ｆ_X ：被測定信号の周波数
Ｆ_Y ：被測定信号の減衰発散係数
ｆ(ｔ₁)〜ｆ(ｔ_n)：被測定信号の観測値(ｎは２以上の整数)
ｔ₁ 〜ｔ_n ：被測定信号のサンプリング波形データの時間

【0132】

また、基本波成分の実効値Ｖ_Ｘは、次の(式３９)で表現できる。

【0133】

【数50】

【0134】

基本波成分が減衰・発散波形であった場合、√２Ｖ_Ｘは、ｔ＝Ｔ／２における瞬間の波形の振幅を表す。

【0135】

ところで、被測定信号の波形の振幅がＡ／Ｄ変換器の最下位ビット付近の値程度にまで零に近くなった場合は、(式３８)の分母も零に近くなり、量子化誤差等の影響によって計算結果が不安定になる。

【0136】

(式３８)の演算結果に５桁の精度を求めるなら、分子・分母とも５桁の有効桁数が必要である。そのため、次の(式４０)によって分母の値を評価し、一定値以下なら周波数の計算結果を無効と判定することが必要である。

【0137】

【数51】

ただし、 ｆ_min：(式３８)の演算結果の有効性を判定するための閾値

【0138】

また、上記(式３５),(式３６)で考慮されている周波数成分の下限はk1／Ｔ、上限はk2／Ｔであるので、実際に計算可能な周波数の範囲は、

である(ただし、Ｔは波形の観測時間幅、Ｍは窓関数の片側の側帯波の個数)。

【0139】

また、被測定信号の波形の周波数が変化している場合は、減衰・発散波形でもないのに(式３８)の計算結果で減衰発散係数Ｆ_Yの値が零にならないことがある。その場合、周波数の計算結果のＦ_X値も十分な精度が得られない。したがって、被測定信号の周波数Ｆ_Xの値にも、その上限の閾値を設けて、計算結果の妥当性をチェックすることが必要である。すなわち、上記(式４０)および次の(式４１)を満たさない(式３８)の計算結果は無効と判定しなければならない。

【0140】

【数52】

ただし、 Ｆ_XMAX：周波数測定限界(上限)閾値
Ｆ_XMIN：周波数測定限界(下限)閾値
Ｆ_YMAX：減衰発散係数の許容限界閾値

【0141】

以下、この発明のフーリエ解析による周波数測定方法および周波数測定装置を図示の実施の形態により詳細に説明する。窓関数の処理をフーリエ変換の後で実施する手法については、手順が前後するだけなので説明を省略する(この場合、窓関数演算部は窓関数効果演算部となる)。

【0142】

〔第１実施形態〕
図１はこの発明の第１実施形態のフーリエ解析による周波数測定方法を用いた周波数測定装置の要部ブロック図を示している。

【0143】

この第１実施形態の周波数測定装置は、図１に示すように、被測定信号の一例としての交流電圧または交流電流が入力される入力変換器１と、入力変換器１により変換された電圧信号をＡ／Ｄ変換するＡ／Ｄ変換器２と、Ａ／Ｄ変換器２により変換されたサンプリング波形データを記憶するメモリー部３と、メモリー部３に記憶されたサンプリング波形データに対して被測定信号の周波数と減衰発散係数を求める処理装置４と、処理装置４により求められた被測定信号の周波数と減衰発散係数を表示する表示器５を備えている。

【0144】

上記処理装置４は、被測定信号のサンプリング波形データに窓関数を掛ける窓関数演算部４aと、窓関数演算部４aで窓関数を掛けたサンプリング波形データに対してフーリエ変換を行うフーリエ変換部４bと、フーリエ変換部４bで得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを上記(式１)に代入して、上記被測定信号の周波数と減衰発散係数を求める周波数演算部４cと、周波数演算部４cで得られた上記被測定信号の周波数と減衰発散係数の有効性を上記(式２)によって判定する判定部４dとを有する。

【0145】

上記周波数測定装置は、Ａ／Ｄ変換器付マイクロコンピュータボードを使用し、そのマイクロコンピュータでソフトウェアを実行することにより演算部を実現している。また、マイクロコンピュータボードにメモリー部３を実装しており、表示器５に用いられたディスプレー端末に演算結果を表示している。

【0146】

この構成は実施形態の一例であり、これらが一筐体に収められていても良く、また有線や無線による通信回線を経由して遠隔地で演算結果を表示しても良い。

【0147】

なお、本発明は、主に商用周波電源の周波数を時々刻々観測することを目的としているが、指数関数的な減衰波形または発散波形でその減衰係数や発散係数を求めることもできる。

【0148】

図４は上記周波数測定装置により測定された周波数の変化の検出例を示している。正弦波の周期が短くなっているところでは、周波数の検出結果が高くなっている。図４において、横軸は時間[ｓｅｃ]を表し、縦軸は周波数[Ｈz]を表している。また、測定された周波数は、３スペクトルで算出したものを実線でプロットし、７スペクトルで算出したものを破線でプロットしている。

【0149】

また、図５は図４に示す周波数変化の要部の拡大図である。図５において、横軸は時間[ｓｅｃ]を表し、縦軸は周波数[Ｈz]を表している。また、測定された周波数は、３スペクトルで算出したものを実線でプロットし、７スペクトルで算出したものを破線でプロットしている。この図５の縦の補助線が２０ｍｓｅｃ間隔であり、周波数が５０Ｈzの交流波形の1サイクルに相当する。図５に示す周波数算出結果のプロットは、３スペクトル方式の場合は各点の前後の１サイクルデータから算出し、７スペクトル方式の場合は各点の前後の３サイクルデータから算出している。

【0150】

〔第２実施形態〕
図２はこの発明の第２実施形態のフーリエ解析による周波数測定方法を用いた周波数測定装置の要部ブロック図を示している。この第２実施形態の周波数測定装置は、フィルタリング方式の周波数測定を行う処理装置を除いて第１実施形態の周波数測定装置と同一の構成をしており、同一構成部には同一参照番号を付している。

【0151】

この第２実施形態の周波数測定装置は、図２に示すように、被測定信号の一例としての交流電圧または交流電流が入力される入力変換器１と、入力変換器１により変換された電圧信号をＡ／Ｄ変換するＡ／Ｄ変換器２と、Ａ／Ｄ変換器２により変換されたサンプリング波形データを記憶するメモリー部３と、メモリー部３に記憶されたサンプリング波形データに対して被測定信号の周波数と減衰発散係数を求める処理装置１４と、処理装置１４により求められた被測定信号の周波数と減衰発散係数を表示する表示器５を備えている。

【0152】

上記処理装置１４は、被測定信号のサンプリング波形データを(式３８)に代入して、被測定信号の周波数と減衰発散係数を求める周波数演算部１４aと、周波数演算部１４aで得られた被測定信号の周波数と減衰発散係数の有効性を(式４０),(式４１)によって判定する判定部１４bとを有する。

【0153】

図３は周波数測定装置の周波数演算部１４aの構成を示している。

【0154】

この周波数演算部１４aは、図３に示すように、サンプリング波形データを表す関数ｆ(k)が入力される乗算器３０と、１サイクル分遅延させるｎ−１個の遅延器３１と、ｎ個の乗算器３２と、ｎ−１個の加算器３３と、ｎ個の乗算器３４と、ｎ−１個の加算器３５と、除算器３６とを有する。上記ｎ個の乗算器３２は、入力側から順に関数ｆ(n)〜ｆ(1)に係数Ａ_n,Ａ_n-1,Ａ_n-2,Ａ_n-3,…,Ａ₁を掛ける。また、上記ｎ個の乗算器３４は、入力側から順に関数ｆ(n)〜ｆ(1)に係数Ｂ_n,Ｂ_n-1,Ｂ_n-2,Ｂ_n-3,…,Ｂ₁を掛ける。

【0155】

上記周波数演算部１４aにおいて、入力されたサンプリング波形データを表す関数ｆ(k)が乗算器３０によりｆ(k)×ｆ(k)が演算してｆ(k)^２を出力して、ｆ(k)^２がｎ−１個の遅延器３１により順次遅延される。これにより、ｆ(k)^２がｋ＝１〜ｎまで揃うと、最終的に除算器３６の出力に基本波成分の周波数出力Ｆ[Ｈz]が出力される。このとき、ｎ−１個の加算器３５のうちの最終段の加算器３５から基本波成分の実効値の二乗値に比例した実効値^２出力が出力される。

【0156】

以降、次のサンプリング波形データであるｆ(k)が周波数演算部１４aに入力される毎に(式３８)の演算が行われることになる。

【0157】

図６は上記第１実施形態の周波数測定装置と第２実施形態の周波数測定装置の周波数測定結果の比較例を示しており、詳しくは、第１実施形態の(式１)による演算結果(破線のグラフ)と第２実施形態の(式３８)による演算結果(実線のグラフ)とを比較したものである。１回の周波数演算には共に商用周波数の８サイクル分にあたる２５６個のサンプリング波形データを用いている。

【0158】

一般に周波数分解能を上げるにはサンプリング波形データの時間幅を長くする必要があるが、データの時間幅を長くすれば周波数の変化に対する追従性能が悪くなる。しかしながら、この発明のフーリエ解析による周波数測定方法および周波数測定装置は、データの時間幅を長くしても追従性が低下しないことが判明した。

【0159】

この発明のフーリエ解析による周波数測定方法によれば、交流波形の１サイクルについて３２サンプル程度の荒いサンプリング周波数でも、６桁程度の精度で周波数が算出できる。特許文献１の特許第２５０５７０７号においては演算式の証明がされておらず、定性的にその特徴を利用しているだけであったが、本発明においては、(式１)の正当性を完全に証明している。

【0160】

さらに、(式１)の正当性の証明に合わせて、フィルタリング方式への理論展開を行っている。このフィルタリング方式は、演算部の計算式が簡単で高速演算が可能な上、ソフトウェアの検証も簡単である。

【0161】

上記第１実施形態のフィルタリング方式の演算は、オフラインで事前に実施可能な係数の演算とリアルタイムで使用する周波数演算式とに分かれており、個別にデバッグできるので、ＦＦＴ演算プログラムを組み込みソフトウェアに実装する場合のような難解さを回避できる。周波数演算式の部分のみ毎回演算が必要であるが、フィルターの計算などでよく用いられるブロック図に即した演算ができれば良く、また、係数は事前にメモリー部内に置いておくことができる。掛け算の回数で評価するならＮポイントのデータに対して４Ｎ回の掛け算で済むので、従前の特許文献１の(式５)の方式で毎回ＦＦＴ演算の後に周波数演算するより、遥かに高速に周波数を算出できる。さらに、(式５)の方式と同程度の精度でかつそれより追従性が良い。

【0162】

この追従性は、(式３５)、(式３６)において、k2を幾らにするかによって決まる。k2の値を大きくすれば追従性は良くなる。近接する２次高調波成分が含まれている場合、その影響による演算誤差が増えるが、データの時間幅を長くすることで周波数分解能を向上させられるので、近接する二次高調波成分による影響を回避できる。また、このフィルタリング方式の周波数測定方法では、演算するデータの時間幅を長くしても、演算結果の周波数の変化に対する追従性能は、殆ど低下しないことが判明し、精度、追従性共に優れた演算手法と言える。

【0163】

【表1】

【0164】

上記第１,第２実施形態では、入力変換器１とＡ／Ｄ変換器２とメモリー部３と処理装置４,１４と表示器５を備えた周波数測定装置について説明したが、周波数測定装置はこれに限らず、入力変換器とＡ／Ｄ変換器とメモリー部や表示器を有しない周波数測定装置でもよく、その場合、所定のサンプリング波形データが直接入力されて、被測定信号の周波数Ｆ_Xと減衰発散係数Ｆ_Yを求めると共にその有効性を判定する。

【0165】

この発明の具体的な実施の形態について説明したが、この発明は上記第１,第２実施形態に限定されるものではなく、この発明の範囲内で種々変更して実施することができる。

【0166】

また、この発明の周波数測定方法は、
交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データを(式３８)に代入して周波数演算部により演算することによって、上記被測定信号の周波数と減衰発散係数を求める第１ステップと、
上記第１ステップで得られた上記被測定信号の周波数と減衰発散係数の有効性を(式４０),(式４１)によって判定部により判定する第２ステップと
を有することを特徴とする。

【0167】

【数53】

【0168】

【数54】

【0169】

【数55】

【0170】

【数56】

【0171】

【数57】

ただし、 Ｆ_X ：被測定信号の周波数
Ｆ_Y ：被測定信号の減衰発散係数
ｆ(ｔ₁)〜ｆ(ｔ_n)：被測定信号の観測値(ｎは２以上の整数)
ｔ₁ 〜ｔ_n ：被測定信号のサンプリング波形データの時間
Ｔ ：被測定信号のサンプリング波形データの時間幅
ｆ_min：(式３８)の演算結果の有効性を判定するための閾値
Ｆ_XMAX：周波数測定限界(上限)
Ｆ_XMIN：周波数測定限界(下限)
Ｆ_YMAX：減衰発散係数の許容限界
Ｗ(t)：窓関数
Ｎ ：データ個数

【0172】

また、この発明の周波数測定装置は、
交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データを(式３８)に代入して演算することによって、上記被測定信号の周波数と減衰発散係数を求める周波数演算部と、
上記周波数演算部で得られた上記被測定信号の周波数と減衰発散係数の有効性を(式４０),(式４１)によって判定する判定部と
を有することを特徴とする。

【0173】

【数58】

【0174】

【数59】

【0175】

【数60】

【0176】

【数61】

【0177】

【数62】

ただし、 Ｆ_X ：被測定信号の周波数
Ｆ_Y ：被測定信号の減衰発散係数
ｆ(ｔ₁)〜ｆ(ｔ_n)：被測定信号の観測値(ｎは２以上の整数)
ｔ₁ 〜ｔ_n ：被測定信号のサンプリング波形データの時間
Ｔ ：被測定信号のサンプリング波形データの時間幅
ｆ_min：(式３８)の演算結果の有効性を判定するための閾値
Ｆ_XMAX：周波数測定限界(上限)
Ｆ_XMIN：周波数測定限界(下限)
Ｆ_YMAX：減衰発散係数の許容限界
Ｗ(t)：窓関数
Ｎ ：データ個数

【符号の説明】

【0178】

１…入力変換器
２…Ａ／Ｄ変換器
３…メモリー部
４…処理装置
５…表示器
４a…窓関数演算部
４b…フーリエ変換部
４c…周波数演算部
４d…判定部
１４…処理装置
１４a…周波数演算部
１４b…判定部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】