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

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

簡易式 位相差計算方法 (Excel) Vo.2 解説

2016-12-17 23:06:41 | 信号処理

 

 今回説明する計算は、「2信号の周波数が異なる」ものについてです。信号2の周波数を高めにしてあります。信号1・信号2の間の「位相差」は少しずつ大きくなります。下図のシミュレーションでは位相差の変化がキリのいいような周波数差をわざとつけています。

手順を順番に解説します。

STEP1 2つの信号を「直交化」します。

 ヒルベルト変換とか・・・必要? とか聞かれそうですが。単一成分で、かつ「周波数既知」なら「差分を用いた巧みな方法」があります。Xn から I,Qをつくるには

I = Xn,  Q = (X[n-1]-X[n+1]) / (2*sinΩ)

とすることで可能です。Excelでやると下図のようになります。 信号1、信号2 それぞれ同じようにIQ化します。 振幅を確かめると確かに合っています。

相対角周波数 Ω は、 Ω = 2πf/fs のように求めます。(fs:サンプリング周波数)

STEP 2-1  作成したIQ信号から「内積」をつくります。

s1 =(I_1, Q_1), s2 = (I_2, Q_2) の複素信号の内積 (Dot)は 

Dot = I_1 * I_2 + Q_1 * Q_2 と計算できます。

STEP 2-2  作成したIQ信号から「外積」をつくります。

s1 =(I_1, Q_1), s2 = (I_2, Q_2) の複素信号の外積 (Cross) は 

Cross = I_1 * Q_2 - Q_1 * I_2 と計算できます。

 

STEP 3 位相差(角)を求めます。

 内積、外積のどちらを使っても位相差 (φ) を求めることができます。下図の例では「外積」をつかったものです。 STEP 2-2 で求めた外積を、 信号1,信号2の振幅の積で割ります。

sin φ = Cross / ( Amp1 * Amp2 )

逆sin を計算して角度が求められます。設定した角度とピタリ一致します。

 

 


 

 オシロスコープにリサージュ図を描かせる方法ですと、異なる周波数では、図形が安定しません。「異なる周波数の間で位相差を測って意味あるのか」とお考えのかたもいらしゃるでしょう。しかし周波数がずれているから位相角も変化していると考えるのが自然です。

 逆に、「位相差が変化しないから、周波数は一致している」ともいえます。今まで私たちは位相差をダイレクトに計算する方法を習っておらず、周波数と時間からスタートしていました。位相角を直接測る方法はリサージュ以降、進歩していなかったのかもしれません。

 代替え的には・・・信号がゼロクロスする時刻をみればおおざっぱな位相差はわかるでしょう。しかし、本計算例で示したとおり、信号が数値で与えられているならば、瞬時の位相差はこんなに簡単に計算できるのです。その「位相差の時間変化率」が「周波数差」であり、 2π・Δf = dφ/dt となるのです。

本計算によると、瞬時の位相差が3サンプルで計算できます。

 今回の例を、連続系で考えてみると非常に単純です。直交化するには微分して振幅を補正すればいいだけ。連続系では微分するとωが係数として飛び出しますが、本例の離散系では sinΩとなっているだけの話です。内積・外積を「幾何学的(つまり図形)」にみれば、2つのベクトルのなす角度をcos でみるかsinでみるかの違いだけです。

 角度を360度シームレスに扱うためには、本例のようにcos, sin を同時にあつかえば2πをまたぐ瞬間、arctanを使った場合はπ/2 での符号の反転、など考えなくてすみます。1回転をこえる場合の回転数の検知は簡単です。 内積・外積を実用的に使いこなすには、±90 度で sin を使う計算、 0-180 度の範囲ならcos を使うものが便利です。また感度の良いところは sin =0, cos =0 の周辺ですので、計算条件を工夫することで精密な角度計算ができるのです。

 今回の解説では、2つの信号の「周波数が既知」としましたが、この周波数がわからなければ・・・なにもできません。

 そんなわけで、Synchro PRIMO法は、なによりも、「瞬時周波数」を計算することを重要視しています。

 「周波数がわかれば/直交化でき/振幅もわかり/ベクトル積を計算すれば角度もわかる」・・・というのが原理なのです。

 


簡易式 位相差計算方法 (Excel) Vo.2 (予告)

2016-12-16 08:01:35 | 信号処理

(予告)・・・記事のアップは数日間お待ちくださいませ。

引き続き、簡易式 位相差計算法をご紹介します。

今後は制約条件がちょっと楽になります。

  • 2信号の周波数は、「それぞれ既知」 ・・・ つまり周波数の異なる波形での位相差計算を試みます。周波数がΔf異なれば、位相角はT秒間に 2πΔf T ずつズレていきます。
  • 今回は、振幅は「未知」でかまいません。

今回、計算する方法はオーソドックスな「ベクトル演算」をベースにしたもので、リサージュ外積は登場しません(たぶん)

 


簡易式 位相差計算方法 (Excel)

2016-12-15 01:09:09 | 信号処理

(謹告)近日中に、本Excelの内容をPythonに移植したものを公開します。 20-Oct-2022.

簡易式位相差計算をExcelでやったものを紹介します。

「簡易式」の制限は以下の通りです

  • 周波数が既知
  • 2つの信号の周波数は同じ、位相差のみある。
  • 2つの信号の振幅が既知
  • 位相差は -90 度〜 90 度まで。 (89度と91度の区別はできません)
こんな感じ。図の左上に設定したパラメータがあります。 位相差は 23.45 deg に設定。
時系列データはだらだらと続きます。

 

 こうすれば、簡単に計算できます!!

STEP1の外積演算 の説明

Cross[n] = X[n-1] * Y[n] - X[n] * Y[n-1] と計算します。 Excel 上ではたすき掛けののように見えます

STEP2 振幅で割る

Q[n] = Cross[n] / (Ax * Ay)  です。本例では Ax =1.0,  Ay=1.3 としています。

STEP3 相対角周波数で割る。 

P[n] = Q[n] / sin(Ω)

= ここの周波数は、既知としてある周波数の「相対角周波数」表記です。 Ω= 2π f /fs と計算します。 

STEP4 上で求めた P[n] が位相差 φ としたときの sinφとなっています。

  φ[n] = arcsin ( P[n] ) / PI * 180   ・・・と計算すると、deg. による角度が求められます。

 


 

今回の説明では、制限が沢山ありましたが、最終的には、任意の2つの波形について、

1)それぞれの周波数、2) それぞれの瞬時振幅 3)瞬時の位相差 

のExcelによる計算例を紹介したいと思います。

本例は、制限は多いですが、従来のような「ゼロクロス点」を推定する必要はまったくありません。

時系列波形から、(本例のような)「外積演算」を用いて導出できる全く新しい方法です。


SynchroPRIMO Ver.3.1 公開

2016-12-14 01:22:58 | 信号処理

 改良とバージョンアップを積み重ねてきたSynchro PRIMOですが、ソフトウェア構造を一新し(別名リファクタリング・・・中がGdgdになってきたので書き換えた・・・)、Ver 3 として公開しました.

 ちょっとした電気工作をし、100V(50/60Hz)を 1Vppくらいに降圧、 信号をオーディオデバイスに突っ込むと、"激安・高性能PMU" のできあがりです。(注:市販品はめちゃ高価です。)計算量は少ないのでラズパイでも動くでしょう。計算方法は、一応解析解を計算するものです!DFTの改良とか、カルマンフィルタなどで誤魔化しているのではありません。(欧米のPMUメーカーさん大丈夫ですか? ww) 

 ただし、AD変換器のサンプリング周波数が測定不確かさ・安定度の決め手となります。時刻あわせはNTPでは無理。GPSの1pps.信号を同時記録する準備中です。

(YouTube)   https://youtu.be/26-dTldBA-c 

(プログラム構造) 

 一応Win32で作りました。(MFCはよく分からない) グラフィックス・GUIは、”ベタ”でコーディングしてあります。

(改良点)

  • グラフィックスのちらつきを解消 (Callback処理でミスっていた・・・
  • オーディオ読み込み/前処理/主処理/グラフィックスCallbackをフル・マルチスレッド化(デバッグには時間かかったです)
  • グラフィックスパネルを統合化 (1:波形とフェーザー、2:周波数専用、3: VFチャートと命名した周波数・電圧の散布図、あと2面を予備に)
  • あまり使っていなかったバッファを廃止。省メモリ。
  • オーディオデバイスからの読み込みを改良。バッファリング最小容量を1/10秒程度に。起動時の遅延を減らす(できればASIOにしたかった)
  • セーブしたデータの解析モードを装備(これで、後でゆっくりと解析できるように)
  • 前処理フィルタのラインナップを整備. Intel Core i5 (2 Core) 程度に最適化(種類ふやしたが、効き目はかわらんかった) 
  • 前処理フィルタは線形でないとあかんかった・・・(Ramp状の周波数変化時、まっすぐなはずの測定値が波うってしまう)
  • 周波数変化率(ROCOF)の計算式を変更し軽量化(一応、最小自乗法によるものです)
  • 波形表示パネルに、小細工。全8CH分のON/OFFをパネルにて操作可能に。
  • モニタした波形の全データセーブ(WAV形式で容量が大きいため、最大24時間)
  • 計測した周波数(Hz)と電圧値(p.u)をまとめてWAVファイルに記録(長時間モニタで威力を発揮)
  • グラフィックスで使用する色は日本の伝統色をベースに使用。(黒だけでも4種類使っています)

(性能)

  • 電力系統周波数測定で、200測定/秒(50Hz), 240測定/秒(60Hz)に改善。おおまかには1/4サイクルで測定。ただし、前処理フィルタの平均時間窓長の関係で、まあこんなもの
  • Synchro PRIMO計算コアモジュール(これが商売ネタ)を改訂。APIを分かりやすく・・・もうじき、Class化できる?

(検証) 

  • 24時間連続運転は安定。(細かいバグがまだあるかもしれません)
  • 24 bit ADCでのモニタをすると・・・VLF帯(およそ8kHz-19kHz)の謎の信号(超スローのFSKにみえる)が多数ひっかかった・・・電力計インバータ由来なのか、「言及してはいけないもの」なのかは不明。
  • パルス状に変化する謎の周波数変動を深夜から早朝にかけ多数検知。位相が瞬間的に1度飛んで(進み方向)いた・・・まれに「2段飛び」もあった。(※系統内で発電所の電力が加わる瞬間?)
  • 市内の数キロはなれた場所で(6.6kVは別変電所)同時モニタしてみたが・・・位相角の変化の様子が合わない。線路インピーダンスの時間変化で、見かけの位相角が影響をうけていただけだった(1度よりはるかに小さい値です) 室内の別コンセントでも似た状況がわずかに発生します。
  • ADCの2CHに同一の信号を入れても、位相差=0とならない。ケーブル長が1m違うと下の桁がまったく合わない。

 

 


5点の値から周波数を求めるExcel計算式

2016-12-13 20:58:34 | 信号処理

 

時系列波形が「こんな具合に」Excel 上にあったとします。せいぜい1/4周期分しかない・・・

「謎の計算式」(クリックして見てください)をはりつけると、周波数がピタリ計算できます。ただし、雑音があると計算不能になります。

やっていることは、リサージュ外積 Ls3, Ls1 の比をとってあとはごにょごにょ・・・・


「リサージュ外積」なるもの

2016-12-13 20:15:15 | 信号処理

 

 リサージュ図形はみなさまご存じと思いますが、「リサージュ外積」とは何でしょうか?

これは、私が勝手に命名した演算です。。。。もともとは、 k段の遅延器の両端の位相差を測ろうとして「リサージュ図形」をとってみたら・・・という発想から生まれたものです。

 k段の遅延を仮定するリサージュ外積 Ls _k(Xn) は下記のように定義されます。Ls_k という一般形式は面倒なので、Ls1, Ls3 を示しますと、 

  Ls1(Xn) = X[n]*X[n] - X[n-1]*X[n+1] 

  Ls3(Xn) = X[n-1]*X[n+1] - X[n-2]*X[n+2]   と、きれいな形をしています。 時系列データXnの自乗に相当する次元をもっており、感覚的にはエネルギーだと思ってください。

 

Xn = A* sin (Ωn - φ) 「ただし、A:振幅, Ω:相対角周波数= ω/fs, φ:位相角」とすると、

Ls1 = A^2 * sin(Ω) * sin(Ω)

Ls3 = A^2 * sin(Ω) *sin(3Ω)    という奇妙な形に変換できます。

 n(つまり時間)によらず、「周波数」と「振幅」のみの関数になるのです。ここがポイントで。 r = Ls3 / Ls1 を計算すると、 r は周波数Ωのみの関数となるわけで、この原理をSynchro PRIMO法は利用しています。

 もう一つの特徴は、 A = A0 * exp( -t/T ) のような形でも、上記のLs3/Ls1 は定数になってしまいます。

「リサージュ外積」なるもの、振幅の自乗と  sin^2(Ω) の積になっています。 元々のリサージュ図での挙動にもどって考えてみると、 リサージュ図形の隣接点Pn-1, Pnと原点がつくる微小な三角形の面積でもあるのです。

 悩みの種は「非線形演算」であること。2成分(Xn = U[n] + W[n] )でやろうろすると、「クロス項」 :

C = U[n]W[n] - U[n-1]*W[n+1] などが登場してしまい、ここが難しく足止めをくらっているところです。。。

 


電力系統周波数の謎

2016-10-06 10:40:53 | 信号処理

 当社で開発したSynchro PRIMO法によるPMUで、電力系統の「観測」を行っています。計測ではなく、電力系統を対象とする「観測・観察」だと思っています。

 一応、給電指令所のほうから「受給バランス」の調整をし周波数を維持していますが、実際の挙動は、電力系統工学の教科書に書かれているようなものではありません。ひたすら測定値を「深夜のライブ」で観察して分かってきたことの一つとして、夜間は周波数が不安定なこと(とは行っても数mHz)程度の一種の周波数変調がおこっており、変調周波数は、1〜4Hz程度でした。こんなに短周期に同期発電機の回転数が変動するわけはなく、理由をずっと考えていましたが・・・答えは簡単でした。

他には、電圧の変動があるのですが、コンセントからモニターしているため、近隣の負荷が原因かもしれません。が、毎日観察していると、一定の時間帯に、同じような電圧形状の変化が発生しているのです・・・

 1秒間に100回測定した値をWAVファイルに記録し、それを音として聞いてみると・・・サンプリング周波数の関係で、実際の変動周期の441倍の高さで聞こえます。・・・ここにも奇妙な音が突然聞こえたり、スイープ状の音もあります。

 瞬時周波数・瞬時電圧・位相角をリアルタイムに、1秒間に100回の測定をしていますが、「奇妙な現象」が他にもいくつも発見されました。

 系統内のどこかの変電所/開閉所で、切り替えをしている様子が観測されました。深夜から早朝にかけて「この現象」がみられます。

 以上のような観察結果をまとめ、電気学会 計測研究会で発表します。 

https://workshop.iee.or.jp/sbtk/cgi-bin/sbtk-showprogram.cgi?workshopid=SBW00004527

 2015年9月に国際会議のため、欧州に行った際、宿泊した各地で「電源の波形」を収集してきました。ニースで収集した波形からは、非常に「奇異な電圧・周波数変動」が測定され、波形を直接観察しても「ゆらゆら」揺れる波形でした。・・・この原因についての考察も発表します。

 精密測定してみると、予想しなかった知見が得られるものです。


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 をご覧ください。