PRIMO信号処理研究所 / Synchro PRIMO Lab.

周波数測定、位相差測定に関する新しい数理。
Please contact to snishie@mac.com.

Synchro PRIMOによる系統周波数(50/60Hz)の観察

2016-04-28 18:52:46 | 信号処理

 当社の方法を用いると、コンセントに来ている交流の周波数(=50/60Hz)の挙動を精密に観測できます。

 普段は意識することはないのですが、電力網には多数の発電機がつながり、また電力消費量は常に変動しています。電力会社では「系統周波数」を正しく50/60Hzに維持するため「需要と同じ量の発電」を常に行っています。

 しかし、秒単位でみると、消費量と発電量は厳密には同じではありません。分単位・1時間単位でつじつまが合うように実際は制御されています。

 では消費量と発電量が合わないとどうなるか?

・ 消費量 > 発電量 だと周波数は低下していきます。

・ 消費量 < 発電量 だと周波数は上昇していきます。

 系統周波数の変化をみることで、瞬間的なロードバランス(消費量と発電量の差)を観察することができます。

 電力系統の周波数の、振幅・周波数・位相角を実時間計測する装置を「同期位相計測装置」(Phasor Measurement Unit = PMU)と呼び、現在さかんに研究開発が行われています。


 下記は当社PMUによる系統周波数のモニター実験です。

(1) 東京電力(50Hz)のモニター

 画面はほぼ実時間で表示されています。

 グラフ:右上が波形、左上が位相角(周波数が50/60Hzより高いと右回りに回る)、左下が系統周波数の"カスケード表示"です。1目盛りが2秒間,1/10秒おきの計測です。

 系統周波数グラフですが、49.9Hz ~ 50.1 Hzのレンジです。グラフの1ピクセルは1mHz(1/1000Hz)を表しており、目盛り1マスが10mHz (1/100 Hz)の変化となっています。東京電力では、めったなことで±100mHzを超えることはないようです。

https://www.youtube.com/watch?v=yR3R7KRuTn4

 

(2) 北陸電力(60Hz)のモニター

 今度は60Hz地区より北陸電力のデータです。収集は富山県高岡市内で夜間。

 周波数変動(グラフ右下)をごらんください。妙に波打っています。変動も東京電力とくらべて乱高下ぎみです。60Hz地域は九州電力・中国電力・四国電力・関西電力・中部電力・北陸電力がずべて「系統連携」されており、周波数は非常に正確に同期されています。 ・・・が、実際には電力会社間の電力の流れ(潮流)の変化の関係で、「系統動揺」と呼ばれる現象が発生しています。正確には動揺ではなく何百キロも離れた「同期発電機」間の同期過渡現象というべきでしょう。しかしこの動揺が大きくなりすぎると様々な問題を引き起こします。

 60Hz地域で実験した経験では、変動が激しいときは±100mHzを飛び越しているシーンを結構みかけます。また比較的急激な変化(といっても、グラフで見ていて、おっとっと・・・と言う感じの10mHz程度の急降下です)が夜間から早朝にかけてよく見られます。 下記のグラフでもそんなシーンが何カ所かあります。

https://www.youtube.com/watch?v=n3DQLr6NHn0


下記は英国における運用の様子を紹介するビデオのようです。

How the National Grid responds to demand 「どのように電力系統は需要に対して対応するか」

https://www.youtube.com/watch?v=UTM2Ck6XWHg

 動画2:10あたりから英国内の系統周波数の実測データが紹介されていますが、ナレーションでいう"フランスからの供給が止まったとき"の周波数低下はなんと -0.3Hz 。一応安全圏内ではありますが、こんな周波数低下は私の日本での実験では出会ったことがありません。

 当社のPMUはこの動画に紹介されているものよりも、はるかに短時間(動画では1秒、当社のものは1/10秒)で同程度の分解能(1mHz)を実現しています。緊急時の周波数低下の検出にはおそらく威力を発揮するでしょう。


 昨年9月、欧州の各地で電力波形を収集してきてまいりました。(フランス南部・イタリア・スイス・ドイツ南部・北部、オーストリア南部・Wien郊外、ベルギー)

 欧州はすべてが50Hzで同一の系統で運用されており、広大なエリアでコンセントから見える周波数はすべて同一なのです。エリアによっては、近隣の発電所が火力中心であったり水力中心であったり、また系統のど真ん中なのか、枝線の末端なのかによって、波形・周波数挙動は全く異なっておりました。 機会あれが紹介したいと思います。


波形データから求める簡易的な位相差計算方法 (2)

2016-04-28 18:49:09 | 信号処理

周波数が既知の波形なら、今回説明する方法で簡単に位相差が計算可能です。

過去の記事 波形データから求める簡易的な位相差計算方法 では、周波数・振幅が既知の場合を紹介しました。今回は周波数のみが既知の場合です。

・ 周波数が正確に分かっていれば半周期のデータでも十分です。

・ ノイズの少ない波形で1サイクルが10サンプル程度で表現されている信号であれば、最低3サンプルで計算可能です。(実用上は1サイクル程度の平均化が好ましい)


手順を紹介します。

・与える信号を x[n], y[n] とします。サンプリング周波数はfsとします

・x,y には位相差φがあり、周波数はf とします。式で書くと;

 x[k] = Ax sin(2πf/fs・k)

 y[k] = Ay sin(2πf/fs・k -φ)   という信号です。

パラメータ A, f, φ が分かっていれば波形 x[n], y[n]は簡単に計算できますが、 その逆の x[n], y[n]からパラメータA,f, φを求めるのは極めて困難です。

 今回の説明では、 「f のみが既知」「周波数が同一の2つの信号」を条件に位相差φと振幅 Ax,Ayの求め方を解説します。


それでは、 x[0],x[1],x[2] .... , y[0],y[1],y[2].....から  振幅 Ax, Ay, 位相差φを求めてみましょう。


STEP1:  x[n],y[n] をそれぞれ「直交化」します。

 私が考案したトリッキーな方法が便利です。

(1) まず秘密の定数 A を以下のように計算します。周波数fから下記のようにθを求め、1/cos θの形の係数Aを求めておきます。

  θ = π/2 - (2π * f/fs)  

  A = 1/(2*cos(θ))   のように"定数A"求めます。

(2) 直交信号の作成:Qを合成します。先ほどのAを使います。

x[n] から (I1[n], Q1[n])を、y[n] から (I2[n],Q2[n])を計算します。

  I1[k] = x[k]

  Q1[k] = (x[k-1] - x[k+1]) *A

・ 信号yも同じように計算し、

  I2[k] = y[k]

  Q2[k] = (y[k-1] - y[k+1]) *A

となります。Aは共通です。

x[n]は  X=(I1[n],Q1[n]) 、 y[n]は  Y[n]=(I2[n],Q2[n])のようにベクトルで表現できるようになりました。


STEP 2: 絶対値の計算

 IQ信号にすると、その絶対値は瞬時振幅となります。自乗和の平方根で振幅が求められます。

   Ax[k] = √(I1[k]^2 + Q1[k]^2) 

   Ay[k] = √(I2[k]^2 + Q2[k]^2) 


STEP 3: 直交化したX,Yを"ベクトル"とみたて、内積D・外積Cを計算します。 

  内積: D[k]= I1[k]・I2[k] + Q1[k]・Q2[k]

  外積: C[k]= I1[k]・Q2[k] - I2[k]・Q1[k]

 位相差が一定の信号であれば, D[k], C[k]は定数になります。実際の計算では多少ばらつくので適当な区間で平均化するといいでしょう。 ノイジーな波形の場合は、平均化して、D, C をもとめておきます。

内積は"ドット積”、外積は"クロス積"とも呼ばれているそうです。


STEP4:  位相角を求めるステップです。

 cos , sin を併用する「Phasor形式」で計算してみます。幾何学的に、"2つのベクトルの作る三角形の面積"は、二辺の長さの積 * cos θ = 1/2 *内積×cosθ、または 1/2*外積×sinθで表現できることを利用しています。言われてみれば当たり前なのですが、この対称性が実はポイントです。

 cos φ = D/ (Ax・Ay)

 sin φ = C/ (Ax・Ay)

cos φ, sin φからφを求めるには用途に応じた逆三角関数の計算を行います。

 arcsinを使うと -π/2 ~ +π/2, arccos を使うと 0 ~ ±π の範囲が使いやすいです。計算が安定しているのは cos φ = 0 の周辺、 sin φ =0 の周辺です。 cos, sin = ±1近傍では 逆三角関数の演算はエラー含みとなります。


 

以上でおわりです。

 今回は「周波数のみ既知」としましたが、別の記事で紹介している「瞬時周波数計算」を使用すれば、周波数が未知の信号間の位相差も計算可能です。また、微小に周波数の異なる2つの波形の位相差が刻々と変化する様子も計算できます。 

今回は概要だけを紹介しましたが、後日、さまざまな計算例を説明したいと思います。


ご質問・ご要望は  snishie [at] synchro-primo.jp までお寄せください。 説明用のExcel Template 他を用意してあります。

Synchro PRIMO 法のコンセプトは当社ホームページ http://synchro-primo.jp をご覧ください。


サンプリング周波数 44.1 kHz の謎

2016-04-05 20:01:56 | 信号処理

 瞬時周波数が理論的に分かる一種の校正用の信号を考えていて一つ発見したことがあります。

 オーディオ用のサンプリング周波数, 44.1 kHz は素因数分解すると、( 2 × 3 × 5 × 7 )^2 となっています。大学院時代に「くし形フィルタ」を扱っており、ノッチ周波数の設定のために適当な遅延段数を設定していたのですが、なにげに割り切れる数が多いので不思議に思っていました。まじめに素因数分解してみると、こんな結果。

 元々は、CDが登場した時の時代背景があり、TVの水平偏向周波数 =15.7 kHz が関与しているのは知っていましたが、こんな形の素数の積になっていたのは驚きです。ネットで検索するとすでに周知の事実でした。

 現在考えているものは、周波数比が、「互いに素」となる周波数成分を混合すると、非常に長い周期の、しかも瞬時周波数が理論値で計算することのできる信号が作成可能なこと。 44.1 kHz のサンプリング周波数において、2Hz, 3Hz, 5Hz, 7Hz の正弦波を混合すると、1秒おきに一定のパターンで繰り返される瞬時周波数の挙動が得られるはずです。またサンプリング周波数となんらかの整数分の1の比となる挙動かもしれません。ただ・・・どのような振幅比でミックスすればよいかは分かりません。

 このような「校正用」の信号定義ができれば、いわゆる「F0推定」の正しさの議論は大きく進歩すると考えています。「F0"推定"であり"測定"はできない、なぜなら真の値が不明だから」といわれ続けていますが、理論値のあきらかな標準信号は定義できれば方法も変わってくるでしょう。

 参考;時間-周波数解析, L.Cohen, 吉川昭訳, 朝倉書店.