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] などが登場してしまい、ここが難しく足止めをくらっているところです。。。