JP2015036635A - フーリエ解析による周波数測定方法および周波数測定装置 - Google Patents

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

Info

Publication number
JP2015036635A
JP2015036635A JP2013167665A JP2013167665A JP2015036635A JP 2015036635 A JP2015036635 A JP 2015036635A JP 2013167665 A JP2013167665 A JP 2013167665A JP 2013167665 A JP2013167665 A JP 2013167665A JP 2015036635 A JP2015036635 A JP 2015036635A
Authority
JP
Japan
Prior art keywords
frequency
signal under
window function
under measurement
equation
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.)
Granted
Application number
JP2013167665A
Other languages
English (en)
Other versions
JP5712255B2 (ja
Inventor
大橋 善和
Yoshikazu Ohashi
善和 大橋
徐 有恒
Aritsune Jo
有恒 徐
英也 大橋
Hideya Ohashi
英也 大橋
達也 牧村
Tatsuya Makimura
達也 牧村
大浦 好文
Yoshifumi Oura
好文 大浦
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Kinkei System Corp
Original Assignee
Kinkei System Corp
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 Kinkei System Corp filed Critical Kinkei System Corp
Priority to JP2013167665A priority Critical patent/JP5712255B2/ja
Publication of JP2015036635A publication Critical patent/JP2015036635A/ja
Application granted granted Critical
Publication of JP5712255B2 publication Critical patent/JP5712255B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Measurement Of Current Or Voltage (AREA)
  • Measuring Frequencies, Analyzing Spectra (AREA)

Abstract

【課題】交流波形の周波数を高速かつ高精度に測定できるフーリエ解析による周波数測定方法を提供する。【解決手段】フーリエ変換部4bで得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ窓関数の持つ周波数成分の個数と少なくとも同じ個数の隣接する低周波数側および高周波数側複素周波数成分とを求める周波数演算部4cと、周波数演算部4cで得られた周波数FXと減衰発散係数FYの有効性を判定する判定部4dとを有する。【選択図】図1

Description

この発明は、フーリエ解析による周波数測定方法および周波数測定装置に関する。
従来、交流波形の周波数を測定する第一の周波数測定方法としては、波形の値が零となる零クロス点を見つけて零クロス点から次の零クロス点までの時間幅を図形的に求めてその時間幅を周期とみなしてその逆数として周波数を求めるものがある。
上記第一の周波数測定方法では、零クロス点周辺にノイズ成分が含まれている場合、正しく周期測定ができない。特に、インバーター機器などの最近のパワーエレクトロニクス技術の普及の結果、電源波形に含まれるノイズ成分は増大する方向にある。
そこで、このようなノイズ成分の影響を受けない従来の第二の周波数測定方法として、波形データの周波数スペクトルから演算で周波数を求めるものがある(例えば、特許第2505707号(特許文献1)参照)。
上記第二の周波数測定方法では、周波数スペクトルを算出するのにFFT演算を使うが、リアルタイム性に欠けており、そのままの形ではオンラインでの使用はできなかった。また、周波数の測定精度を上げるためには、交流波形の数サイクルに渡る波形に窓関数を掛けてFFT演算をしなければならないが、演算結果を波形の急激な変化に追従させるには無理があった。
一般に時間幅Tのデータをフーリエ解析(FFT(高速フーリエ変換)またはDFT(離散フーリエ変換))して周波数成分分析を行なった場合、もとの波形に含まれる振動成分が商用周波数成分一つだけの場合でもそれがFFTによって得られる各スペクトルの周波数(FFTのBIN周波数という)からずれていると多数のスペクトルに分解される。
例えば、100msecのデータをFFT分析して10Hz毎の周波数スペクトルを算出した場合、商用周波数成分が51Hzとしたときに得られる結果は、50Hz周辺に出現する多数のスペクトルで表現されることになる。
そもそもフーリエ変換の解析結果は、そこに各周波数成分が有ることを保証するものではなく、解析結果の周波数成分の波形の和がそのデータの区間で元の波形にもっとも近い(二乗平均誤差が最小になる)ことを示しているだけである。
上記特許文献1には、一個の周波数成分が複数のスペクトルに分散した場合、元の波形の周波数が幾らだったのかを算出する方法が書かれている。
この特許文献1は、被測定信号の周波数FXを(式5)で算出するもので、減衰発散係数FYも得ることができる。
Figure 2015036635
ただし、 FX :被測定信号の周波数
Y :被測定信号の減衰発散係数
1 〜Vn :被測定信号の複素周波数成分
1 〜Fn :複素周波数成分V1 〜Vn の各周波数
ここで、波形を表す複素関数f(t)は、実用上は実数関数であるが、時間tの関数であり、次の(式6)のように複素フーリエ級数で表現できるものである。
Figure 2015036635
ただし、Vkは第k番目の複素周波数成分、T[sec]は観測データが存する時間幅である。周波数スペクトルの間隔を△f[Hz]とすると、△f=1/Tの関係がある。また、第k番目のスペクトルの周波数Fk[Hz]は、
k=k△f=k/T
であり、角周波数ωk[rad/sec]とは、
ωk=k△ω=2πk△f
の関係がある。
この第二の周波数測定方法によれば、確かに特定時間断面における周波数成分の中心周波数を荒いサンプリング波形データ(例えば32サンプル/交流の1サイクル程度)でも精度よく求めることができる。
特許第2505707号
しかしながら、上記第二の周波数測定方法(特許文献1)では、特定の時間断面での値だけではなく、周波数の変化をリアルタイムで観測したい場合、一定の時間間隔(例えば2msec毎)で何度もこの手法を適用してFFTを繰り返さねばならず、処理時間が掛かるためにリアルタイム性を持たせることには難点があった。
また、上記第二の周波数測定方法を装置に実装するには、FFTプログラムを組み込む必要があり、種類豊富な演算専用チップを使用する場合、チップ専用のマシン語を習得して直接プログラミングすることは困難なので、どうしてもC言語のようなある程度規格化された高級言語を使わざるを得なかったが、高速性が要求されるFFT演算に高級言語を用いることは、高速性が犠牲になるので好ましいことではない。
そこで、この発明の課題は、交流波形の2サイクル〜4サイクル程度で、かつ、交流波形の1サイクルあたり数十サンプルといった比較的低いサンプリング周波数によりサンプリングされたサンプリング波形データからその交流波形の周波数を高速かつ高精度に測定できるフーリエ解析による周波数測定方法および周波数測定装置を提供することにある。
上記課題を解決するため、第1の発明のフーリエ解析による周波数測定方法は、
交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに対して窓関数演算部により窓関数を掛ける第1ステップと、
上記第1ステップで上記窓関数を掛けた上記サンプリング波形データに対してフーリエ変換部によりフーリエ変換を行う第2ステップと、
上記第2ステップで得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを(式1)に代入して周波数演算部により演算することによって、上記被測定信号の周波数FXと減衰発散係数FYを求める第3ステップと、
上記第3ステップで得られた上記被測定信号の周波数FXと減衰発散係数FYの有効性を(式2)に基づいて判定部により判定する第4ステップと
を有することを特徴とする。
ただし、窓関数演算部については後述の方法により、フーリエ変換を先に実施して、その結果得られる複素スペクトル間の簡単な演算で最初に窓関数を掛けたと同じ効果を奏する手法を適用しても良い。
Figure 2015036635
Figure 2015036635
ただし、 Vk :時間幅T[秒]の被測定信号のV(t)をフーリエ変換して得られる
k/T[Hz]の複素周波数成分(ここでk=1,2,…)
m :上記複素周波数成分のうち絶対値が最大のもの
n :窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
X :被測定信号の周波数
Y :被測定信号の減衰発散係数
T :被測定信号のサンプリング波形データの演算対象部分の時間幅
min:(式1)の演算結果の有効性を判定するための閾値
上記構成によれば、交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに窓関数を掛けた後、窓関数を掛けたサンプリング波形データに対してフーリエ変換を行って、そのフーリエ変換により得られた複素周波数スペクトル値の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接複素周波数成分とを(式1)に代入して、上記被測定信号の周波数と減衰発散係数を求める。また、上記被測定信号の周波数FXと減衰発散係数FYの有効性を(式2)に基づいて判定して、(式2)の条件を満たさないなら周波数FXと減衰発散係数FYの演算結果を無効にすることで、量子化誤差等の影響によって不安定になった計算結果を除外できる。これにより、交流波形の2サイクル〜4サイクル程度で、かつ、交流波形の1サイクルあたり数十サンプルといった比較的低いサンプリング周波数によりサンプリングされたサンプリング波形データからその交流波形の周波数を高速かつ高精度に測定できる。特に、窓関数に含まれる周波数成分の個数相当分だけのスペクトルを演算に用いているので、不要な周波数成分を算出する必要がなく、フーリエ変換に際してFFTを用いて全周波数成分(例えば256サンプルデータなら128個の複素スペクトルデータ)を算出する必要がなく、DFTを用いて必要な周波数成分(窓関数内の周波数成分が1個なら中心周波数成分のスペクトル1個と低周波数側および高周波数側にそれぞれ1個〜2個、全部で3個〜5個のスペクトルで十分)のみを算出することで演算速度を遥かに高速化することができたのである。
また、第2の発明の周波数測定装置は、
交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに窓関数を掛ける窓関数演算部と、
上記窓関数演算部で上記窓関数を掛けた上記サンプリング波形データに対してフーリエ変換を行うフーリエ変換部と、
上記フーリエ変換部で得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを(式1)に代入して演算することによって、上記被測定信号の周波数FXと減衰発散係数FYを求める周波数演算部と、
上記周波数演算部で上記(式1)により得られた上記被測定信号の周波数FXと減衰発散係数FYの有効性を(式2)に基づいて判定する判定部と
を備えたことを特徴とする。
Figure 2015036635
Figure 2015036635
ただし、 Vk :時間幅T[秒]の被測定信号のV(t)をフーリエ変換して得られる
k/T[Hz]の複素周波数成分(ここでk=1,2,…)
m :上記複素周波数成分のうち絶対値が最大のもの
n :窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
X :被測定信号の周波数
Y :被測定信号の減衰発散係数
T :被測定信号のサンプリング波形データの演算対象部分の時間幅
min:(式1)の演算結果の有効性を判定するための閾値
上記構成によれば、第1の発明と同様の効果を奏する。
また、第3の発明のフーリエ解析による周波数測定方法は、
交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに対してフーリエ変換部によりフーリエ変換を行う第1ステップと、
上記第1ステップで得られた複素周波数成分のデータに対して窓関数効果演算部により(式4)による演算を行って、上記第1ステップの前に上記被測定信号のサンプリング波形データに(式3)で表された窓関数を掛けた演算結果と同等の演算結果を得る第2ステップと、
上記第2ステップで得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを(式1)に代入して周波数演算部により演算することによって、上記被測定信号の周波数FXと減衰発散係数FYを求める第3ステップと、
上記第3ステップで得られた上記被測定信号の周波数FXと減衰発散係数FYの有効性を(式2)に基づいて判定部により判定する第4ステップと
を有することを特徴とする。
Figure 2015036635
Figure 2015036635
ただし、 Vk :時間幅T[秒]の被測定信号のV(t)をフーリエ変換して得られる
k/T[Hz]の複素周波数成分(ここでk=1,2,…)
m :上記複素周波数成分のうち絶対値が最大のもの
n :窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
X :被測定信号の周波数
Y :被測定信号の減衰発散係数
T :被測定信号のサンプリング波形データの演算対象部分の時間幅
min:(式1)の演算結果の有効性を判定するための閾値
Figure 2015036635
Figure 2015036635
ただし、 w(t):窓関数
α:窓関数の直流成分の係数
β、γ、δ:窓関数の各周波数成分の係数
Ck :窓関数を掛けない状態でフーリエ変換した結果の
k/T[Hz]成分のフーリエ係数
C’k :上記窓関数w(t)を掛けてフーリエ変換した結果の
k/T[Hz]成分のフーリエ係数
上記構成によれば、交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに対してフーリエ変換を行った後、フーリエ変換で得られた複素周波数成分のデータに対して窓関数効果演算部により(式4)による演算を行って、フーリエ変換の前に被測定信号のサンプリング波形データに(式5)で表された窓関数を掛けた演算結果と同等の演算結果を得るので、第1の発明と同様の効果を奏する。また、最初に窓関数を掛ける必要がないので、1サンプル毎に周波数を算出する必要がある場合に再帰的DFTを利用することによって、演算速度を飛躍的に改善できる。
また、第4の発明の周波数測定装置は、
交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに対してフーリエ変換を行うフーリエ変換部と、
上記フーリエ変換部で得られた複素周波数成分のデータに対して(式4)による演算を行って、上記フーリエ変換前に上記被測定信号のサンプリング波形データに(式3)で表された窓関数を掛けた演算結果と同等の演算結果を得る窓関数効果演算部と、
上記窓関数効果演算部で得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを(式1)に代入して演算することによって、上記被測定信号の周波数FXと減衰発散係数FYを求める周波数演算部と、
上記周波数演算部で上記(式1)により得られた上記被測定信号の周波数FXと減衰発散係数FYの有効性を(式2)に基づいて判定する判定部と
を備えたことを特徴とする。
Figure 2015036635
Figure 2015036635
ただし、 Vk :時間幅T[秒]の被測定信号のV(t)をフーリエ変換して得られる
k/T[Hz]の複素周波数成分(ここでk=1,2,…)
m :上記複素周波数成分のうち絶対値が最大のもの
n :窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
X :被測定信号の周波数
Y :被測定信号の減衰発散係数
T :被測定信号のサンプリング波形データの演算対象部分の時間幅
min:(式1)の演算結果の有効性を判定するための閾値
Figure 2015036635
Figure 2015036635
ただし、 w(t):窓関数
α:窓関数の直流成分の係数
β、γ、δ:窓関数の各周波数成分の係数
Ck :窓関数を掛けない状態でフーリエ変換した結果の
k/T[Hz]成分のフーリエ係数
C’k :上記窓関数w(t)を掛けてフーリエ変換した結果の
k/T[Hz]成分のフーリエ係数
上記構成によれば、第3の発明と同様の効果を奏する。
以上より明らかなように、この発明によれば、交流電圧波形または交流電流波形を表す被測定信号の2サイクル〜4サイクル程度で、かつ、交流波形の1サイクルあたり数十サンプルといった比較的低いサンプリング周波数によりサンプリングされたサンプリング波形データからでも、被測定信号の周波数を高速かつ高精度に測定することができるフーリエ解析による周波数測定方法および周波数測定装置を実現することができる。
図1はこの発明の第1実施形態のフーリエ解析による周波数測定方法を用いた周波数測定装置の要部ブロック図である。 図2はこの発明の第2実施形態のフーリエ解析による周波数測定方法を用いた周波数測定装置の要部ブロック図である。 図3は上記周波数測定装置の周波数演算部の構成を示す図である。 図4は上記周波数測定装置により測定された周波数の変化の検出例を示す図である。 図5は図4に示す周波数変化の要部の拡大図である。 図6は上記第1実施形態の周波数測定装置と上記第2実施形態の周波数測定装置の周波数測定結果の比較例を示す図である。
この発明のフーリエ解析による周波数測定方法および周波数測定装置について、まず、特許第2505707号(特許文献1)に記載された数式(1)であって本明細書における(式5)の算出根拠について説明し、その後に上述の課題を解決する手段および実施形態の詳細を説明する。
いま、交流波形を表す複素関数f(t)を、
Figure 2015036635
とし、複素関数f(t)が一個の周波数成分として、
Figure 2015036635
を持っていたとする。このときの未知の角周波数成分ωxが、
ωx=2πFk または Fk=ωx/2π
のように解析周波数Fkのいずれかと一致している場合は、その成分のみが検出され、他の成分は零なので(式5)が成り立つのは明らかである。
そうでない場合、複素関数f(t)の解析結果は、(式5)のように複数の複素周波数成分(スペクトル成分)の和で表されるが、無限個の計算はできないので、何らかの窓関数を掛けることで第k1個目から第k2個目(k2>k1)までのスペクトルを除く他のスペクトルを十分小さくし、第k1個目から第k2個目までのスペクトルで充分な精度でかつ角周波数ωxが複素数の範囲において、次の(式7)が成立するものとする。
Figure 2015036635
この(式7)は時間tが−tでも成立するので、次の(式8)が得られる。
Figure 2015036635
上記(式7)と(式8)の辺々を掛け合わせると、次の(式9)となる。ここで、(式9)の右辺の第二項のΣの記号の下部の式は、「kおよびlが異なる整数値で、かつ、k1以上k2以下となる全ての組み合わせの総和」を示すものとする。
Figure 2015036635
上記(式9)の左辺のAおよび右辺の第一項が時間tに関係しない定数項であるので、右辺の第二項も定数であることが判る。また、△ω=2π/Tのため、区間t=[0、T]で積分すると、
Figure 2015036635
で表される複素数ベクトルは、原点の周りを(k-l)周(ただし、k≠l)するので、そのベクトル和は零になり、区間t=[0、T]での時間積分値が零であることが判る。
すなわち、(式9)の両辺を時刻t、観測区間[0、T]で積分すると、次の(式10)および(式11)で表現できるが、(式11)の右辺の第二項が零になることが判る。
Figure 2015036635
Figure 2015036635
それ故、(式9)の右辺の第二項は、時間tに関係ない定数であって、その時間積分値が零であることから、もともと常に零であることが判る。よって、(式9)は、
Figure 2015036635
となる。
一方、(式7)の両辺をtで微分することにより、次の(式12)が得られる。
Figure 2015036635
この(式12)に(式8)の両辺を掛け合わせると、(式13)が得られる。
Figure 2015036635
この(式13)の右辺の第二項も(式9)と同様の理由で零になる。
また、左辺のAは各周波数成分の二乗和なので、以下の(式14)が得られる。
Figure 2015036635
ただし、波形自体が0値で無入力状態であれば当然のことながら、その周波数を求めることは意味がない。よって、波形の実効値がある値以下なら(式14)での計算結果も無効であると判断する機能も必要である。一般に歪波の実効値は、各周波数成分の実効値の二乗和の平方根で計算されるが、(式14)の分母は丁度その実効値ベクトルの二乗値になっている。それ故、この点については、(式14)の分母の絶対値がA/D変換器の分解能に比して十分に大きな値であることを保証することで足る。
一方、(式14)と同様の考え方で(式8)の両辺を掛け合わせる代わりに(式12)の辺々を二乗して整理することで次の(式15)も得られる。
Figure 2015036635
(式15)でも同様の結果が得られるが、この場合は得られるのがωxの二乗値であり、その平方根は2個出てくるので、別途判別する必要がある。(式14)の方が簡単で使い勝手が良いので、通常は(式14)または角周波数の代わりに単に周波数で表現した(式5)を使う。等号を用いているが、(式5)や(式14)で等号が成り立つのは、第k1番目〜第k2番目(k1<k2)以外の周波数成分が十分小さく無視できる場合に限る。
一般に、交流波形の周波数がBIN周波数からずれていると、そのスペクトラムは多数の周波数成分に散逸するが、窓関数を掛けるとそれを一定の範囲内に抑えることができる。窓関数は、4〜5桁程度の精度が必要な場合は4サイクル分のサンプリング波形データにハニング窓を掛け、6〜8桁の精度が必要な場合は8サイクル分のサンプリング波形データにナットール窓やブラックマン・ハリス窓を掛ける。
以上が(式5)および(式14)の導出の根拠であるが、この計算をPLD(プログラマブルロジックデバイス)などで実現しようとすると、各周波数成分を求めるFFT(高速フーリエ変換)演算などが複雑で実施困難である。DSP(Digital Signal Processor:デジタル・シグナル・プロセッサ)のような演算専用チップを用いると(式5),(式14)の計算も実施可能であるが、高速性を重視したマシン語によるプログラミングは困難である。このような演算専用チップを用いた計算で高速性を出すためには、例えばC言語のような高級言語によって記述されたソフトウェアをコンパイルして使用するのは効率的ではなく、アセンブリ言語やマシン語で記述するのが適切である。しかしながら、様々な演算用チップが市販されているので、それらを統一的に同じソースコードでプログラミングすることは不可能で、個別にソフトウェアを開発し、多大な時間を掛けてデバッグしているのが現状であった。
<フーリエ変換方式について>
商用周波数の交流波形からその周波数を求める場合は、FFT解析結果のスペクトルの周波数の一つが50Hzまたは60Hzに一致するように、サンプリング周波数とFFT演算を行うデータの時間幅(データの個数)とを選定し、商用周波数に近くかつ最も大きな複素周波数成分のスペクトルとその前後に1本か2本のスペクトルを選んでそれらの複素数値を後述する(式1)に代入すれば良いので、敢えてFFT演算を行ってすべての周波数成分を求める必要はなく、DFT(離散フーリエ変換)演算で十分である。
例えば、50Hzの交流波形を64サンプル/交流1サイクルのサンプリングで256個のサンプリング波形データを得て、それにFFT演算を施すと128本のスペクトルが得られるが、この場合のスペクトル間の周波数間隔は50×64/256=12.5Hzである。
ハニング窓関数を用いた場合、交流波形v(t)は以下の(式16)のように表現できるが、これは(式17)のように3つの周波数成分に分解できる。
Figure 2015036635
ただし、fは交流波形の周波数、Tは観測波形データの時間幅である。
したがって、この3本のスペクトルを表す複素数値を(式5)に代入すれば、商用周波の交流波形から周波数が算出できる。万一、それが本来の周波数(例えば50Hz)からずれていれば、この3本のスペクトルの外側にも零でないスペクトルが現れる。それが、(式5)の演算対象から外れると誤差が発生する。その場合、中心周波数の両側に1本ではなく2本または3本のスペクトルを演算対象とすると誤差は小さくなる。
一般的に窓関数の例としては、以下の(式18)の形をしており、定数成分α以外にM個(M=1〜3程度の整数)の周波数成分を持っている。そのため、それが正弦波交流波形の式と掛け合わされると、その個数M×2個の側帯波を生ずることになる。それ故、それら全てを(式5)の演算対象にすれば精度上問題のない結果が得られることが判る。当然のことながら、測定対象の波形は周波数変動するのであるから、その変動幅も考慮する必要がある。
[窓関数の例]
Figure 2015036635
[係数値の具体例]
・ハニング窓の場合
Figure 2015036635
・ブラックマン・ハリス窓の場合
Figure 2015036635
・ナットール窓の場合
Figure 2015036635
これらの式から明らかなように、例えばハニング窓では周波数成分は直流分αを含めて2個であるので、これを用いて交流波形の周波数成分を分析する場合は、中心周波数のスペクトル(中心複素周波数成分)と、その中心周波数のスペクトルによりも低周波数側の1本のスペクトル(隣接する低周波数側の複素周波数成分)と、高周波数側の1本のスペクトル(隣接する高周波数側の複素周波数成分)の計3本のスペクトルが得られれば計算可能である。また、ブラックマン・ハリス窓やナットール窓では計7本である。
しかるに、スペクトルの本数が増えるとDFTの演算量が増えるので、計算速度が遅くなってしまう。また、スペクトル間の周波数間隔を狭めようとするとデータの時間幅を大きくしなければならず、そうなると周波数の急変に追従させることが難しくなる。そのため、本発明者は、追従性が要求される用途では、スペクトルの本数を極力絞って計算速度と追従性能を改善し、また、追従性能を要求されない用途では、窓関数を考慮したDFTの計算部分をあらかじめ実行しておき、実時間ではより簡単な計算式で周波数を算出する方法を考えた。すなわち、高速性・高追従性を有する方法と、高速性・高精度を実現できる方法とを考案し検証した。
次に、スペクトルの本数を極力絞って計算速度と追従性能を改善する方法について説明する。
サンプリング波形データの時間幅を交流波形の本来の周期の2倍にとり、交流波形2サイクル分のサンプリング波形データについて3本のスペクトルを計算する。交流波形の本来の周波数が50Hzの場合は、25Hz、50Hz、75Hzの3本のスペクトルをDFTの計算式(式19)によって計算すると、次の(式20)に示す計算結果となる。
Figure 2015036635
ここで、k=1,2,3である。
Figure 2015036635
しかし、直流成分や2次高調波成分などが含まれていると、それが(式20)での計算結果に誤差として含まれ、それを分離することはできない。
そのような場合には、サンプリング波形データを交流波形の周期の4倍から8倍に渡って収集し、中心周波数の前後の2〜4本、計5〜9本のスペクトルについて(式20)の計算を行えばよい。
ただし、そうすると、スペクトルの本数やデータ数が共に2〜4倍になるので、演算に必要な時間は4〜16倍になってしまい、演算速度の低下が避けられない。
そのような場合には、係数計算を先に済ませておいて、メモリーに保存しておくことで高速演算を可能にする方法が必要となる。
このように、本発明のフーリエ解析による周波数測定方法は、第1ステップで交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに窓関数を掛けた後、第2ステップで窓関数を掛けたサンプリング波形データに対してフーリエ変換を行って、第3ステップでそのフーリエ変換により得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを次の(式1)に代入して、上記被測定信号の周波数と減衰発散係数を求める。さらに、第4ステップで被測定信号の周波数FXと減衰発散係数FYの有効性を次の(式2)によって判定して、(式2)の条件を満たさないなら周波数FXと減衰発散係数FYの演算結果を無効にすることで、量子化誤差等の影響によって不安定になった計算結果を除外できる。これにより、交流波形の2サイクル〜4サイクル程度で、かつ、交流波形の1サイクルあたり数十サンプルといった比較的低いサンプリング周波数によりサンプリングされたサンプリング波形データからその交流波形の周波数を高速かつ高精度に測定することができる。
Figure 2015036635
Figure 2015036635
ただし、 Vk :時間幅T[秒]の被測定信号のV(t)をフーリエ変換して得られる
k/T[Hz]の複素周波数成分(ここでk=1,2,…)
m :上記複素周波数成分のうち絶対値が最大のもの
n :窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
X :被測定信号の周波数
Y :被測定信号の減衰発散係数
T :被測定信号のサンプリング波形データの演算対象部分の時間幅
min:(式1)の演算結果の有効性を判定するための閾値
ところで、窓関数を掛けることは、周波数スペクトルを一定の周波数範囲に収め、本発明における周波数の演算式における演算項数を限定するうえで必要なステップであるが、先に窓関数を掛けてその後にフーリエ変換することに限定する必要はなく、先にフーリエ変換を行ってその結果に対して例えば(式4)による演算を行うことでも同様の効果を得ることができる。
その根拠について以下に説明する。例えば、角周波数ω0が、
Figure 2015036635
で表され、サンプリング波形データの時間幅Tを商用周波の周期の4倍にとって窓関数w(t)を掛けた場合、フーリエ係数C'kは以下の(式21)のようになる。ここで、f(t)は交流波形を表す複素関数、w(tn)は窓関数、Nはデータ個数である。
Figure 2015036635
窓関数w(tn)としてハニング窓を用いた場合、それを代入すると以下の(式22)のように変形できる。
Figure 2015036635
これは、ハニング窓を掛けた場合のFFTの結果のフーリエ係数C’kと、窓関数を掛けない場合のFFTの結果のフーリエ係数Ckとの間に以下の(式23)の関係があることを示している。
Figure 2015036635
同様に、(式3)の窓関数w(t)に対して(式4)が成立する。これらの結果は、必ずしも先に窓関数を掛けてからフーリエ変換する必要がないことを示しており、使用する窓関数をフーリエ変換後に選択することも可能であることや複数の窓関数の効果の違いを調べる場合にその度毎にフーリエ変換し直す必要がないことを示している。また、最初に窓関数を掛ける必要がなければ再帰的DFTが利用でき、1サンプル毎に周波数を算出する必要がある場合は演算速度が飛躍的に改善される。
Figure 2015036635
Figure 2015036635
よって、本発明では、窓関数を掛ける順序(フーリエ変換の前または後)を限定しないこととする。
<フィルタリング方式について>
一方、本発明者は、上記(式1)が観測値f(t)を用いて例えば(式24)のような形になれば、簡単な論理デバイスでも実現可能と考えて理論式を検討し、各係数を算出して図3に示すブロック構成で周波数演算機能を実現した。
Figure 2015036635
図3の構成は、入力部の二乗演算と出力部の割り算以外は、一般のフィルターなどで良く使用される演算回路であるので、それらのファームウエアがそのまま利用できる。この手法も結局は(式1)を変形したものであり、原理は同じであるが、演算の手順が異なるので周波数変化時の過渡特性は異なっている。
上記(式1)を用いるには、サンプリング波形データからDFTやFFTなどの手法で複素周波数成分を求める手順が必要であり、元の波形データには窓関数を掛けることが必要であった。しかし、これでは毎回窓関数を掛けてFFTをする必要が有り、演算時間が掛かってリアルタイム性が悪いので、本発明者は、(式1)を時間関数上に展開したのである。
複素フーリエ級数展開の式から(式1)のVkは、次の(式25)〜(式28)のように表現できる。
Figure 2015036635
例えば、60Hz±10Hz程度の商用周波数の波形データをサンプリング周波数1920Hz(32サンプル/サイクル)で8サイクル分サンプリングし、256個のサンプリング波形データからDFTやFFTによって60/8=7.5Hz毎に得られる周波数成分スペクトルのうちの60Hzを中心とする1本と、その前後のそれぞれ3本のスペクトル(計7本のスペクトル)から(式5)で周波数を求めるものとする。
この場合、データ個数はN=256、サンプリング波形データの時間幅はT=0.1333…[sec]であり、各スペクトルは7.5Hz、15Hz、22.5Hz、30Hz、37.5Hz、45Hz、52.5Hz、60Hz、67.5Hz、75Hz、82.5Hz、90Hz、…が得られる。このうちの37.5Hz〜82.5Hz(k=5〜11)の周波数成分スペクトルを用いて基本波の周波数を算出する。
ところで、次数の異なる高調波の周波数成分の和の実効値の二乗は、各次数の周波数成分の実効値の二乗和に等しいので、(式1)の各周波数成分は以下の(式29)で表現できる。この(式29)の右辺の括弧内は倍角公式を用いて(式30)のように変形できる。
Figure 2015036635
よって、(式1)の分子は(式31)で表現できると共に、
Figure 2015036635
(式1)の分母は(式32)で表現できることが判る。
Figure 2015036635
上記(式31)および(式32)でΣ記号による加算演算の順序を入れ替えている点は重要である。この入れ替えによって、これまで各周波数におけるDFT演算を各波形データに対して行ってその結果を加算していたものが、各周波数におけるDFT演算の係数を先に加算して置いておき、後で波形データに掛け合わせることができるので、毎回時間のかかる演算をする必要がなく、実際のオンライン演算での演算時間を短縮することができる。
ここで、(式1)の分子の式の係数をAnで表し、分母の式の係数をBnで表すと、(式33)および(式34)が得られる。
Figure 2015036635
Figure 2015036635
この(式33)、(式34)の値は前もって計算しておくことのできる値である。実際の計算では、これに更に窓関数W(2πn/N)を掛けて、(式35)や(式36)の値を計算しておく。
Figure 2015036635
Figure 2015036635
ただし、Wは窓関数であり、例えばハニング窓の場合は(式37)のようになるが、これ以外の窓関数を用いてもよい。
Figure 2015036635
求める周波数FXと減衰発散係数FYは、上記(式35),(式36)の係数を用いて、次の(式38)で算出できる。上記(式35),(式36)は、入力される波形データによって変化する値ではないので、最初に一回係数Anと係数Bnを計算するだけで済む。その計算された係数Anと係数Bnをメモリー部に記憶して次の計算で読み出して用いる。なお、上記(式35),(式36)の計算結果を不揮発性メモリーに保存しておいてもよい。
Figure 2015036635
ただし、 FX :被測定信号の周波数
Y :被測定信号の減衰発散係数
f(t1)〜f(tn):被測定信号の観測値(nは2以上の整数)
1 〜tn :被測定信号のサンプリング波形データの時間
また、基本波成分の実効値Vは、次の(式39)で表現できる。
Figure 2015036635
基本波成分が減衰・発散波形であった場合、√2Vは、t=T/2における瞬間の波形の振幅を表す。
ところで、被測定信号の波形の振幅がA/D変換器の最下位ビット付近の値程度にまで零に近くなった場合は、(式38)の分母も零に近くなり、量子化誤差等の影響によって計算結果が不安定になる。
(式38)の演算結果に5桁の精度を求めるなら、分子・分母とも5桁の有効桁数が必要である。そのため、次の(式40)によって分母の値を評価し、一定値以下なら周波数の計算結果を無効と判定することが必要である。
Figure 2015036635
ただし、 fmin:(式38)の演算結果の有効性を判定するための閾値
また、上記(式35),(式36)で考慮されている周波数成分の下限はk1/T、上限はk2/Tであるので、実際に計算可能な周波数の範囲は、
Figure 2015036635
である(ただし、Tは波形の観測時間幅、Mは窓関数の片側の側帯波の個数)。
また、被測定信号の波形の周波数が変化している場合は、減衰・発散波形でもないのに(式38)の計算結果で減衰発散係数FYの値が零にならないことがある。その場合、周波数の計算結果のFX値も十分な精度が得られない。したがって、被測定信号の周波数FXの値にも、その上限の閾値を設けて、計算結果の妥当性をチェックすることが必要である。すなわち、上記(式40)および次の(式41)を満たさない(式38)の計算結果は無効と判定しなければならない。
Figure 2015036635
ただし、 FXMAX:周波数測定限界(上限)閾値
XMIN:周波数測定限界(下限)閾値
YMAX:減衰発散係数の許容限界閾値
以下、この発明のフーリエ解析による周波数測定方法および周波数測定装置を図示の実施の形態により詳細に説明する。窓関数の処理をフーリエ変換の後で実施する手法については、手順が前後するだけなので説明を省略する(この場合、窓関数演算部は窓関数効果演算部となる)。
〔第1実施形態〕
図1はこの発明の第1実施形態のフーリエ解析による周波数測定方法を用いた周波数測定装置の要部ブロック図を示している。
この第1実施形態の周波数測定装置は、図1に示すように、被測定信号の一例としての交流電圧または交流電流が入力される入力変換器1と、入力変換器1により変換された電圧信号をA/D変換するA/D変換器2と、A/D変換器2により変換されたサンプリング波形データを記憶するメモリー部3と、メモリー部3に記憶されたサンプリング波形データに対して被測定信号の周波数と減衰発散係数を求める処理装置4と、処理装置4により求められた被測定信号の周波数と減衰発散係数を表示する表示器5を備えている。
上記処理装置4は、被測定信号のサンプリング波形データに窓関数を掛ける窓関数演算部4aと、窓関数演算部4aで窓関数を掛けたサンプリング波形データに対してフーリエ変換を行うフーリエ変換部4bと、フーリエ変換部4bで得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを上記(式1)に代入して、上記被測定信号の周波数と減衰発散係数を求める周波数演算部4cと、周波数演算部4cで得られた上記被測定信号の周波数と減衰発散係数の有効性を上記(式2)によって判定する判定部4dとを有する。
上記周波数測定装置は、A/D変換器付マイクロコンピュータボードを使用し、そのマイクロコンピュータでソフトウェアを実行することにより演算部を実現している。また、マイクロコンピュータボードにメモリー部3を実装しており、表示器5に用いられたディスプレー端末に演算結果を表示している。
この構成は実施形態の一例であり、これらが一筐体に収められていても良く、また有線や無線による通信回線を経由して遠隔地で演算結果を表示しても良い。
なお、本発明は、主に商用周波電源の周波数を時々刻々観測することを目的としているが、指数関数的な減衰波形または発散波形でその減衰係数や発散係数を求めることもできる。
図4は上記周波数測定装置により測定された周波数の変化の検出例を示している。正弦波の周期が短くなっているところでは、周波数の検出結果が高くなっている。図4において、横軸は時間[sec]を表し、縦軸は周波数[Hz]を表している。また、測定された周波数は、3スペクトルで算出したものを実線でプロットし、7スペクトルで算出したものを破線でプロットしている。
また、図5は図4に示す周波数変化の要部の拡大図である。図5において、横軸は時間[sec]を表し、縦軸は周波数[Hz]を表している。また、測定された周波数は、3スペクトルで算出したものを実線でプロットし、7スペクトルで算出したものを破線でプロットしている。この図5の縦の補助線が20msec間隔であり、周波数が50Hzの交流波形の1サイクルに相当する。図5に示す周波数算出結果のプロットは、3スペクトル方式の場合は各点の前後の1サイクルデータから算出し、7スペクトル方式の場合は各点の前後の3サイクルデータから算出している。
〔第2実施形態〕
図2はこの発明の第2実施形態のフーリエ解析による周波数測定方法を用いた周波数測定装置の要部ブロック図を示している。この第2実施形態の周波数測定装置は、フィルタリング方式の周波数測定を行う処理装置を除いて第1実施形態の周波数測定装置と同一の構成をしており、同一構成部には同一参照番号を付している。
この第2実施形態の周波数測定装置は、図2に示すように、被測定信号の一例としての交流電圧または交流電流が入力される入力変換器1と、入力変換器1により変換された電圧信号をA/D変換するA/D変換器2と、A/D変換器2により変換されたサンプリング波形データを記憶するメモリー部3と、メモリー部3に記憶されたサンプリング波形データに対して被測定信号の周波数と減衰発散係数を求める処理装置14と、処理装置14により求められた被測定信号の周波数と減衰発散係数を表示する表示器5を備えている。
上記処理装置14は、被測定信号のサンプリング波形データを(式38)に代入して、被測定信号の周波数と減衰発散係数を求める周波数演算部14aと、周波数演算部14aで得られた被測定信号の周波数と減衰発散係数の有効性を(式40),(式41)によって判定する判定部14bとを有する。
図3は周波数測定装置の周波数演算部14aの構成を示している。
この周波数演算部14aは、図3に示すように、サンプリング波形データを表す関数f(k)が入力される乗算器30と、1サイクル分遅延させるn−1個の遅延器31と、n個の乗算器32と、n−1個の加算器33と、n個の乗算器34と、n−1個の加算器35と、除算器36とを有する。上記n個の乗算器32は、入力側から順に関数f(n)〜f(1)に係数An,An-1,An-2,An-3,…,A1を掛ける。また、上記n個の乗算器34は、入力側から順に関数f(n)〜f(1)に係数Bn,Bn-1,Bn-2,Bn-3,…,B1を掛ける。
上記周波数演算部14aにおいて、入力されたサンプリング波形データを表す関数f(k)が乗算器30によりf(k)×f(k)が演算してf(k)を出力して、f(k)がn−1個の遅延器31により順次遅延される。これにより、f(k)がk=1〜nまで揃うと、最終的に除算器36の出力に基本波成分の周波数出力F[Hz]が出力される。このとき、n−1個の加算器35のうちの最終段の加算器35から基本波成分の実効値の二乗値に比例した実効値出力が出力される。
以降、次のサンプリング波形データであるf(k)が周波数演算部14aに入力される毎に(式38)の演算が行われることになる。
図6は上記第1実施形態の周波数測定装置と第2実施形態の周波数測定装置の周波数測定結果の比較例を示しており、詳しくは、第1実施形態の(式1)による演算結果(破線のグラフ)と第2実施形態の(式38)による演算結果(実線のグラフ)とを比較したものである。1回の周波数演算には共に商用周波数の8サイクル分にあたる256個のサンプリング波形データを用いている。
一般に周波数分解能を上げるにはサンプリング波形データの時間幅を長くする必要があるが、データの時間幅を長くすれば周波数の変化に対する追従性能が悪くなる。しかしながら、この発明のフーリエ解析による周波数測定方法および周波数測定装置は、データの時間幅を長くしても追従性が低下しないことが判明した。
この発明のフーリエ解析による周波数測定方法によれば、交流波形の1サイクルについて32サンプル程度の荒いサンプリング周波数でも、6桁程度の精度で周波数が算出できる。特許文献1の特許第2505707号においては演算式の証明がされておらず、定性的にその特徴を利用しているだけであったが、本発明においては、(式1)の正当性を完全に証明している。
さらに、(式1)の正当性の証明に合わせて、フィルタリング方式への理論展開を行っている。このフィルタリング方式は、演算部の計算式が簡単で高速演算が可能な上、ソフトウェアの検証も簡単である。
上記第1実施形態のフィルタリング方式の演算は、オフラインで事前に実施可能な係数の演算とリアルタイムで使用する周波数演算式とに分かれており、個別にデバッグできるので、FFT演算プログラムを組み込みソフトウェアに実装する場合のような難解さを回避できる。周波数演算式の部分のみ毎回演算が必要であるが、フィルターの計算などでよく用いられるブロック図に即した演算ができれば良く、また、係数は事前にメモリー部内に置いておくことができる。掛け算の回数で評価するならNポイントのデータに対して4N回の掛け算で済むので、従前の特許文献1の(式5)の方式で毎回FFT演算の後に周波数演算するより、遥かに高速に周波数を算出できる。さらに、(式5)の方式と同程度の精度でかつそれより追従性が良い。
この追従性は、(式35)、(式36)において、k2を幾らにするかによって決まる。k2の値を大きくすれば追従性は良くなる。近接する2次高調波成分が含まれている場合、その影響による演算誤差が増えるが、データの時間幅を長くすることで周波数分解能を向上させられるので、近接する二次高調波成分による影響を回避できる。また、このフィルタリング方式の周波数測定方法では、演算するデータの時間幅を長くしても、演算結果の周波数の変化に対する追従性能は、殆ど低下しないことが判明し、精度、追従性共に優れた演算手法と言える。
Figure 2015036635
上記第1,第2実施形態では、入力変換器1とA/D変換器2とメモリー部3と処理装置4,14と表示器5を備えた周波数測定装置について説明したが、周波数測定装置はこれに限らず、入力変換器とA/D変換器とメモリー部や表示器を有しない周波数測定装置でもよく、その場合、所定のサンプリング波形データが直接入力されて、被測定信号の周波数FXと減衰発散係数FYを求めると共にその有効性を判定する。
この発明の具体的な実施の形態について説明したが、この発明は上記第1,第2実施形態に限定されるものではなく、この発明の範囲内で種々変更して実施することができる。
また、この発明の周波数測定方法は、
交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データを(式38)に代入して周波数演算部により演算することによって、上記被測定信号の周波数と減衰発散係数を求める第1ステップと、
上記第1ステップで得られた上記被測定信号の周波数と減衰発散係数の有効性を(式40),(式41)によって判定部により判定する第2ステップと
を有することを特徴とする。
Figure 2015036635
Figure 2015036635
Figure 2015036635
Figure 2015036635
Figure 2015036635
ただし、 FX :被測定信号の周波数
Y :被測定信号の減衰発散係数
f(t1)〜f(tn):被測定信号の観測値(nは2以上の整数)
1 〜tn :被測定信号のサンプリング波形データの時間
T :被測定信号のサンプリング波形データの時間幅
min:(式38)の演算結果の有効性を判定するための閾値
XMAX:周波数測定限界(上限)
XMIN:周波数測定限界(下限)
YMAX:減衰発散係数の許容限界
W(t):窓関数
N :データ個数
また、この発明の周波数測定装置は、
交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データを(式38)に代入して演算することによって、上記被測定信号の周波数と減衰発散係数を求める周波数演算部と、
上記周波数演算部で得られた上記被測定信号の周波数と減衰発散係数の有効性を(式40),(式41)によって判定する判定部と
を有することを特徴とする。
Figure 2015036635
Figure 2015036635
Figure 2015036635
Figure 2015036635
Figure 2015036635
ただし、 FX :被測定信号の周波数
Y :被測定信号の減衰発散係数
f(t1)〜f(tn):被測定信号の観測値(nは2以上の整数)
1 〜tn :被測定信号のサンプリング波形データの時間
T :被測定信号のサンプリング波形データの時間幅
min:(式38)の演算結果の有効性を判定するための閾値
XMAX:周波数測定限界(上限)
XMIN:周波数測定限界(下限)
YMAX:減衰発散係数の許容限界
W(t):窓関数
N :データ個数
1…入力変換器
2…A/D変換器
3…メモリー部
4…処理装置
5…表示器
4a…窓関数演算部
4b…フーリエ変換部
4c…周波数演算部
4d…判定部
14…処理装置
14a…周波数演算部
14b…判定部
Figure 2015036635
ただし、 Vk :時間幅T[秒]の被測定信号のV(t)をフーリエ変換して得られる
k/T[Hz]の複素周波数成分(ここでk=1,2,…)
m :上記複素周波数成分のうち絶対値が最大のもの
n :窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
X :被測定信号の周波数
Y :被測定信号の減衰発散係数
T :被測定信号のサンプリング波形データの演算対象部分の時間幅
min:(式1)の演算結果の有効性を判定するための閾値
Figure 2015036635
ただし、 Vk :時間幅T[秒]の被測定信号のV(t)をフーリエ変換して得られる
k/T[Hz]の複素周波数成分(ここでk=1,2,…)
m :上記複素周波数成分のうち絶対値が最大のもの
n :窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
X :被測定信号の周波数
Y :被測定信号の減衰発散係数
T :被測定信号のサンプリング波形データの演算対象部分の時間幅
min:(式1)の演算結果の有効性を判定するための閾値
Figure 2015036635
ただし、 Vk :時間幅T[秒]の被測定信号のV(t)をフーリエ変換して得られる
k/T[Hz]の複素周波数成分(ここでk=1,2,…)
m :上記複素周波数成分のうち絶対値が最大のもの
n :窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
X :被測定信号の周波数
Y :被測定信号の減衰発散係数
T :被測定信号のサンプリング波形データの演算対象部分の時間幅
min:(式1)の演算結果の有効性を判定するための閾値
Figure 2015036635
ただし、 Vk :時間幅T[秒]の被測定信号のV(t)をフーリエ変換して得られる
k/T[Hz]の複素周波数成分(ここでk=1,2,…)
m :上記複素周波数成分のうち絶対値が最大のもの
n :窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
X :被測定信号の周波数
Y :被測定信号の減衰発散係数
T :被測定信号のサンプリング波形データの演算対象部分の時間幅
min:(式1)の演算結果の有効性を判定するための閾値
Figure 2015036635
ただし、 Vk :時間幅T[秒]の被測定信号のV(t)をフーリエ変換して得られる
k/T[Hz]の複素周波数成分(ここでk=1,2,…)
m :上記複素周波数成分のうち絶対値が最大のもの
n :窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
または、その周波数成分の個数よりも大きい整数
X :被測定信号の周波数
Y :被測定信号の減衰発散係数
T :被測定信号のサンプリング波形データの演算対象部分の時間幅
min:(式1)の演算結果の有効性を判定するための閾値

Claims (4)

  1. 交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに対して窓関数演算部により窓関数を掛ける第1ステップと、
    上記第1ステップで上記窓関数を掛けた上記サンプリング波形データに対してフーリエ変換部によりフーリエ変換を行う第2ステップと、
    上記第2ステップで得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを(式1)に代入して周波数演算部により演算することによって、上記被測定信号の周波数FXと減衰発散係数FYを求める第3ステップと、
    上記第3ステップで得られた上記被測定信号の周波数FXと減衰発散係数FYの有効性を(式2)に基づいて判定部により判定する第4ステップと
    を有することを特徴とするフーリエ解析による周波数測定方法。
    Figure 2015036635
    Figure 2015036635
    ただし、 Vk :時間幅T[秒]の被測定信号のV(t)をフーリエ変換して得られる
    k/T[Hz]の複素周波数成分(ここでk=1,2,…)
    m :上記複素周波数成分のうち絶対値が最大のもの
    n :窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
    または、その周波数成分の個数よりも大きい整数
    X :被測定信号の周波数
    Y :被測定信号の減衰発散係数
    T :被測定信号のサンプリング波形データの演算対象部分の時間幅
    min:(式1)の演算結果の有効性を判定するための閾値
  2. 交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに窓関数を掛ける窓関数演算部と、
    上記窓関数演算部で上記窓関数を掛けた上記サンプリング波形データに対してフーリエ変換を行うフーリエ変換部と、
    上記フーリエ変換部で得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを(式1)に代入して演算することによって、上記被測定信号の周波数FXと減衰発散係数FYを求める周波数演算部と、
    上記周波数演算部で上記(式1)により得られた上記被測定信号の周波数FXと減衰発散係数FYの有効性を(式2)に基づいて判定する判定部と
    を備えたことを特徴とする周波数測定装置。
    Figure 2015036635
    Figure 2015036635
    ただし、 Vk :時間幅T[秒]の被測定信号のV(t)をフーリエ変換して得られる
    k/T[Hz]の複素周波数成分(ここでk=1,2,…)
    m :上記複素周波数成分のうち絶対値が最大のもの
    n :窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
    または、その周波数成分の個数よりも大きい整数
    X :被測定信号の周波数
    Y :被測定信号の減衰発散係数
    T :被測定信号のサンプリング波形データの演算対象部分の時間幅
    min:(式1)の演算結果の有効性を判定するための閾値
  3. 交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに対してフーリエ変換部によりフーリエ変換を行う第1ステップと、
    上記第1ステップで得られた複素周波数成分のデータに対して窓関数演算部により(式4)による演算を行って、上記第1ステップの前に上記被測定信号のサンプリング波形データに(式3)で表された窓関数を掛けた演算結果と同等の演算結果を得る第2ステップと、
    上記第2ステップで得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを(式1)に代入して周波数演算部により演算することによって、上記被測定信号の周波数FXと減衰発散係数FYを求める第3ステップと、
    上記第3ステップで得られた上記被測定信号の周波数FXと減衰発散係数FYの有効性を(式2)に基づいて判定部により判定する第4ステップと
    を有することを特徴とするフーリエ解析による周波数測定方法。
    Figure 2015036635
    Figure 2015036635
    ただし、 Vk :時間幅T[秒]の被測定信号のV(t)をフーリエ変換して得られる
    k/T[Hz]の複素周波数成分(ここでk=1,2,…)
    m :上記複素周波数成分のうち絶対値が最大のもの
    n :窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
    または、その周波数成分の個数よりも大きい整数
    X :被測定信号の周波数
    Y :被測定信号の減衰発散係数
    T :被測定信号のサンプリング波形データの演算対象部分の時間幅
    min:(式1)の演算結果の有効性を判定するための閾値
    Figure 2015036635
    Figure 2015036635
    ただし、 w(t):窓関数
    α:窓関数の直流成分の係数
    β、γ、δ:窓関数の各周波数成分の係数
    Ck :窓関数を掛けない状態でフーリエ変換した結果の
    k/T[Hz]成分のフーリエ係数
    C’k :上記窓関数w(t)を掛けてフーリエ変換した結果の
    k/T[Hz]成分のフーリエ係数
  4. 交流電圧波形または交流電流波形を表す被測定信号のサンプリング波形データに対してフーリエ変換を行うフーリエ変換部と、
    上記フーリエ変換部で得られた複素周波数成分のデータに対して(式4)による演算を行って、上記フーリエ変換前に上記被測定信号のサンプリング波形データに(式3)で表された窓関数を掛けた演算結果と同等の演算結果を得る窓関数演算部と、
    上記窓関数演算部で得られた複素周波数成分の絶対値が最大のものを中心とする中心複素周波数成分と、その中心複素周波数成分に隣接する低周波数側および高周波数側にあってかつ上記窓関数の持つ周波数成分の個数と少なくとも同じ個数かまたはそれ以上の個数の隣接する低周波数側および高周波数側の複素周波数成分とを(式1)に代入して演算することによって、上記被測定信号の周波数FXと減衰発散係数FYを求める周波数演算部と、
    上記周波数演算部で上記(式1)により得られた上記被測定信号の周波数FXと減衰発散係数FYの有効性を(式2)に基づいて判定する判定部と
    を備えたことを特徴とする周波数測定装置。
    Figure 2015036635
    Figure 2015036635
    ただし、 Vk :時間幅T[秒]の被測定信号のV(t)をフーリエ変換して得られる
    k/T[Hz]の複素周波数成分(ここでk=1,2,…)
    m :上記複素周波数成分のうち絶対値が最大のもの
    n :窓関数に含まれる直流成分を除く周波数成分の個数と同じか、
    または、その周波数成分の個数よりも大きい整数
    X :被測定信号の周波数
    Y :被測定信号の減衰発散係数
    T :被測定信号のサンプリング波形データの演算対象部分の時間幅
    min:(式1)の演算結果の有効性を判定するための閾値
    Figure 2015036635
    Figure 2015036635
    ただし、 w(t):窓関数
    α:窓関数の直流成分の係数
    β、γ、δ:窓関数の各周波数成分の係数
    Ck :窓関数を掛けない状態でフーリエ変換した結果の
    k/T[Hz]成分のフーリエ係数
    C’k :上記窓関数w(t)を掛けてフーリエ変換した結果の
    k/T[Hz]成分のフーリエ係数
JP2013167665A 2013-08-12 2013-08-12 フーリエ解析による周波数測定方法および周波数測定装置 Active JP5712255B2 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2013167665A JP5712255B2 (ja) 2013-08-12 2013-08-12 フーリエ解析による周波数測定方法および周波数測定装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2013167665A JP5712255B2 (ja) 2013-08-12 2013-08-12 フーリエ解析による周波数測定方法および周波数測定装置

Publications (2)

Publication Number Publication Date
JP2015036635A true JP2015036635A (ja) 2015-02-23
JP5712255B2 JP5712255B2 (ja) 2015-05-07

Family

ID=52687192

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2013167665A Active JP5712255B2 (ja) 2013-08-12 2013-08-12 フーリエ解析による周波数測定方法および周波数測定装置

Country Status (1)

Country Link
JP (1) JP5712255B2 (ja)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106872773A (zh) * 2017-04-25 2017-06-20 中国电子科技集团公司第二十九研究所 一种单载频脉冲信号的多脉冲频率精确测量方法及装置
CN109472051A (zh) * 2018-10-11 2019-03-15 天津大学 基于短时傅里叶变换的硬件木马检测方法
CN113341219A (zh) * 2021-05-19 2021-09-03 北京航空航天大学 恒频交流供电系统频率调制幅度测量方法及装置
CN114034927A (zh) * 2021-11-01 2022-02-11 南京国电南自电网自动化有限公司 一种跟频插值采样的信号测量方法及系统
KR20220041273A (ko) * 2020-09-24 2022-04-01 한국전력공사 상업시설 영업시간 예측모델 생성 장치 및 방법

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2505707B2 (ja) * 1993-11-12 1996-06-12 株式会社近計システム フ―リエ解析による周波数測定方法
JPH10504647A (ja) * 1994-08-25 1998-05-06 シーメンス アクチエンゲゼルシヤフト 電力供給系統の監視方法および装置

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2505707B2 (ja) * 1993-11-12 1996-06-12 株式会社近計システム フ―リエ解析による周波数測定方法
JPH10504647A (ja) * 1994-08-25 1998-05-06 シーメンス アクチエンゲゼルシヤフト 電力供給系統の監視方法および装置

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106872773A (zh) * 2017-04-25 2017-06-20 中国电子科技集团公司第二十九研究所 一种单载频脉冲信号的多脉冲频率精确测量方法及装置
CN109472051A (zh) * 2018-10-11 2019-03-15 天津大学 基于短时傅里叶变换的硬件木马检测方法
CN109472051B (zh) * 2018-10-11 2023-07-04 天津大学 基于短时傅里叶变换的硬件木马检测方法
KR20220041273A (ko) * 2020-09-24 2022-04-01 한국전력공사 상업시설 영업시간 예측모델 생성 장치 및 방법
KR102411743B1 (ko) * 2020-09-24 2022-06-24 한국전력공사 상업시설 영업시간 예측모델 생성 장치 및 방법
CN113341219A (zh) * 2021-05-19 2021-09-03 北京航空航天大学 恒频交流供电系统频率调制幅度测量方法及装置
CN114034927A (zh) * 2021-11-01 2022-02-11 南京国电南自电网自动化有限公司 一种跟频插值采样的信号测量方法及系统

Also Published As

Publication number Publication date
JP5712255B2 (ja) 2015-05-07

Similar Documents

Publication Publication Date Title
JP5712255B2 (ja) フーリエ解析による周波数測定方法および周波数測定装置
Chen et al. Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions
Wen et al. Simple interpolated FFT algorithm based on minimize sidelobe windows for power-harmonic analysis
Xia et al. Single-phase phase angle measurements in electric power systems
Huang et al. Enhancement of digital equivalent voltage flicker measurement via continuous wavelet transform
Bilski et al. Virtual spectrum analyzer based on data acquisition card
US8407268B2 (en) Method for determining an optimum sampling frequency, and a power analyzer performing the method
Radonjic et al. Stochastic measurement of power grid frequency using a two-bit A/D converter
CN109374966A (zh) 一种电网频率估计方法
KR100911685B1 (ko) 코히어런트하지 않게 샘플링된 데이타의 파워 스펙트럼을측정하기 위한 저누설 방법
Xue et al. Subspace-least mean square method for accurate harmonic and interharmonic measurement in power systems
CN106405254B (zh) 一种低频段电磁环境分析方法及装置
Zhang et al. A new regularized adaptive windowed Lomb periodogram for time–frequency analysis of nonstationary signals with impulsive components
JP6090000B2 (ja) 周波数解析装置
Nuccio et al. Assessment of virtual instruments measurement uncertainty
JP5971425B1 (ja) 交流信号解析装置、交流信号解析方法及びプログラム
CN109584902B (zh) 一种音乐节奏确定方法、装置、设备及存储介质
JP2505707B2 (ja) フ―リエ解析による周波数測定方法
CN106501582B (zh) 信号基线处理方法及信号基线处理设备
Wetula A Hilbert transform based algorithm for detection of a complex envelope of a power grid signals-an implementation
JP2010185682A (ja) 一般調和解析装置および周波数分析装置
Nunzi et al. A procedure for highly reproducible measurements of ADC spectral parameters
JP2020118496A (ja) 振動分析装置、振動分析方法及びプログラム
Bravo et al. Development of a smart wavelet-based power quality monitoring system
Srinivasan et al. Multi-sine EIS-drift, non linearity and solution resistance effects

Legal Events

Date Code Title Description
A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20141120

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20150303

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20150309

R150 Certificate of patent or registration of utility model

Ref document number: 5712255

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

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

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

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250