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

　時系列波形 x[n], y[n]から位相差θを求める
簡易的な「位相差計算方法」を説明します

　サンプリング周波数fs, 周波数ωと振幅Ax,Ay がわかっている
ものとします。

（方法）

1) 周波数を正規化します。 Ω = ω/fs とします. (ω=2πf)

2) 以下のLsを計算します。

　　Ls = x[n] * y[n-1] - x[n-1] * y[n]

3) 位相差 sinθは

sin θ = Ls / (Ax * Ay *sin Ω)

ともとめられます。
　
　上記 Ls は 「波形のどの部分で計算しても」同じ値と
なります。必要あれば、広い区間で平均化して使うといいです。

振幅、周波数が既知である２つの波形の位相差は、上の計算式を
用いて位相差が ±90度の範囲で求めることができます。

（周期の大きな波形の場合）

１周期のサンプル数が大きい場合、Lsの計算が安定しません。
この場合、2)の処理で、「隣のサンプル」を広くとることができます。

　　Ls = x[n] * y[n-K] - x[n-K] * y[n] のようにLsを計算するなら
　　sin θ = Ls / (Ax * Ay *sin KΩ)

で計算できます。

（証明）

　x[n] = Ax * sin (Ωn)
　y[n] = Ay * sin(Ωn -θ）　

とおいて三角関数の各種公式で展開すると、計算式を導出できます

（問題）

　しかし、ちょっと不満がありますね。
・　振幅と周波数が「既知」でないといけない
・　位相差が ±π/2　の範囲でしか扱えない。

簡易的には利用価値が高いと思いますので、興味あるかたはどうぞ。

もし、周波数は分かっていて、振幅は未知の場合、

波形データから求める簡易的な位相差計算方法 （２）　をお試しください。

遅延器を使うIQ信号作成方法

周波数既知の信号が時系列波形 x[n] で与えられているとします。
この信号をIQ化したい。　どうするか？

周波数がわかっているなら、１周期あたりのサンプル数を計算し、
その1/4 の段数の遅延を加えると 90度の移相は完成します。

しかし、１周期のサンプル数が４の倍数にならなければどうするか。
または、整数にならない場合はどうするか？

（２つの遅延器）

　ユニークな方法を紹介します。

1) 遅延器　D1,D2　を直列につなぎます。
　それぞれK段の遅延。計2K段。
2)　D1への入力を x0, 　　D1の出力を x1,
　　D2の出力を x2 とします。

3) I 信号は D1出力をそのまま取り出します。
4) Q信号はちょっとトリッキーな方法で合成します。

　入力された信号の周波数をΩで表現すると、

Q= (x2- x0 ) * A
ただし A = 0.5 * 1/cosθ ; θ =KΩ-π/2

　となります。 周波数を「相対角周波数」で表現するのがコツです。

（９０度の移相器）

　この式を見ると、x2-x0 は90度の移相ができていることを
示しています。係数Aは振幅をあわせるものです。

こうやって、「周波数がわかってさえいれば」　ベクトル化(IQ化)　
できます。　そうすると・・・・

前回説明した、ベクトル演算（内積、外積）で位相差は簡単に
求められます。　今回説明した方法では周波数を知っておく必要が
ありますが、これも、すでに説明した方法で周波数Ωをもとめる
ことができます。

（位相差測定の手順）

まとめると、

１）信号１、信号２の周波数Ωをはかっておき、
２）今回説明した方法で IQ信号に変換、
３）ベクトル演算で、位相角がcosφ, sinφ　で得られます。

信号２に「基準信号」を設定すると、その基準からの位相角が
得られます。　同一の基準信号から測定した位相角の差分を
とれば位相差になります。

離散系における周波数の表現方法

離散系における周波数の表現方法には便利なやりかたがあります。

信号 y(t) が、 y(t)=A * sin (ωt + φ )　などと定義されていても、
サンプリング周波数を意識してその都度 t を計算しながら波形を
求めるのは面倒です。

実は、サンプリング周波数fs を基準とする「相対周波数」や
「相対角周波数」を使うと計算が非常に楽になります。


（相対周波数）

大文字 F をつかって, F = f / fs とします。サンプリングの限界である
ナイキスト周波数は F= 0.5 と書くことができます。

角周波数は ω =2πf と同様にω ->　Ωと大文字にして、
Ω = 2πF　とします。

（離散系表記）

　相対周波数をつかうと、各サンプル n における値は

y[n] = A * sin( Ωn - φ ) と書くことができます。　この式なら
Excelなんかでも簡単に波形が計算できそうですね。

（例）

fs = 44100 Hz 、 f=441 Hz なら、F= 0.01 となります。　
実は、Fの逆数は、「１周期あたりのサンプル数」と等しくなります。

ーーー

　わたしが書いた論文では、多くの場合、この相対周波数表記を
使用しています。

当社の位相差測定方法

位相差測定・計算のための公式は知られておりません。

　例えば 時系列波形 x,y 　を与えて、φ =phase ( x ,y ) のような
位相差計算ができるとこれは画期的なことです。

しかし、位相差測定（計算）で、周波数が異なる信号の位相差は
どうなるか？リサージュ図を使って位相差を求める方法は教科書に
かかれていますが、周波数が異なり図形が安定しないと、リサージュを
使う方法は不可能です。

　DFT/FFTを使えば一発ではないか、という意見もあります。
しかし角度の引き算は以外に厄介で位相差が変化する場合（すなわち
周波数座が微妙にあある）は面倒です。

（ベクトル図による方法）

　ここでベクトル図の基本に立ち返ってみます。　２つのベクトル
p,q があったとすると、p, q のなす角度はどうやれば計算できるか？

　これは簡単で、p ,q の内積を求めて；
cos φ = (p・ q) / (|p| *|q| ) とやればOKです。
p =( x1, y1 ), q = (x2, y2) とすると、角度φは；

cos φ = (x1 * x2 + y1 * y2 ) / {sqrt( x1^2 + y1^2 ) sqrt(x2^2 + y2^2 )}
と計算できます。

外積を使うと、今度はsinφの形式で
sin φ = (p × q) / (|p| *|q| )
sin φ = (x1 * y2 - x2 * y1 ) / {sqrt( x1^2 + y1^2 ) sqrt(x2^2 + y2^2 )}
となります。

　もし、あつかう信号がベクトルであるならば、位相差は、
内積・外積・絶対値から求めることができます。90度付近はcos に
よるものが有利です。 0度付近では外積を使うものが計算しやすく
なります。

（課題）

　しかし、測定したい信号は通常実数です。IQ信号を扱えるところは
少ないです。　もし、実信号を「非常に簡単な方法」でベクトル化
（IQ信号）に変換できれば、以上説明した「ベクトル」による非常に
シンプルな位相差計算ができることになります。

（既知の方法）

　直交信号を作成するためには「ヒルベルト変換」を使うことができます。
時系列信号を X_re であたえると、X_im = Hilbert ( X_re )のように一発で
計算できます。(Matlab等では標準）
しかし、この中身はあるインパルス応答の畳み込み積分。計算量は少なく
ありません。

（当社の方法）

　当社で考案した方法では「遅延器」をつかって90度信号をシフトする
方法です。遅延時間と移相量は比例することを利用しますが、問題は
90度をどうやって正確にコントロールするか。

90度の移相の方法；
「遅延器を使うIQ信号作成方法」をご覧ください

５サンプルの離散波形による瞬時周波数計算

ここで当社で考案した方法を紹介します。


１）　周波数を計算したい波形から、５サンプルをとりだします。　
　１周期がちょうど１００サンプルとなる波形を用意しておいて、
　任意の位置から５サンプルとりだしてみましょう。
　　もし、サンプリング周波数が 44100 Hzとすると、この波形の
　実周波数は441Hzになります。

　時系列波形は；
　x[n] = sin ( ２* PI * 441 * (n/44100) ) と計算できますね。
　位相差項をいれても構いません。

２）　とりだしたサンプルルを、順番に x1, x2, x3, x4 ,x5
　　とします。

３）　以下のような　Ls1　と　Ls3を計算します。

Ls1 = x3*x3 - x2*x4
Ls3 = x2*x4 - x1*x5

単純な、一種のたすき掛け演算です。

４） つぎに、Ls1, Ls3 から、 r = Ls3 / Ls1 のように
　　比をとります。

５）以上のように求めた r から次の計算で周波数（瞬時周波数）が
　　計算できます。

　 f = (fs / 2 PI) * arcsin ( 0.5 sqrt (3- r) )
fs はサンプリング周波数です。

・　時系列波形が十分に有効数字のある浮動小数点でしたら
　　完全に一致します。
・　正弦波であれば、「どの部分を切り出しても」正しく
　　求められます。
・　５サンプルが十分に瞬間とみなせれば周波数の変化も
　　捉えることができます。

※　実際の応用では様々な工夫が必要です。
※　上の例は原理をご理解いただくためのものです。


興味あるかたはお試しください。

発表しました／発表します

　当研究所の最新成果についてです

■　電気学会　計測研究会 2015.6.25 金沢 で発表しました。

「遅延器による擬似直交変換とベクトル演算を用いる高時間分解能な
　同期位相・周波数・振幅計測方法」

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


■　EUSIPCO 2015 で発表します。
2015.8.31 - 9.4 Nice, France.
　"" Principle and Implementation of Vector-based PMU Algorithm using Delay Devices"
Suminori Nishie (PRIMO Signal Processing Laboratory, Japan)

http://www.eusipco2015.org/content/technical-programme

開業！ 当研究所について.

　当研究所は、洗練された独自開発のアルゴリズムによる
「精密周波数測定」「精密位相差測定」　を専門とする会社です。
すべて数式で記述することのできるアルゴリズムで、計算量が非常に
少なく、応用先の要求にあわせ多様な実装が可能なものです

（PRIMO とは）

　Primary Rotation and Instantaneous Movement Observation の略で、
　瞬時の「ある基本回転」と「その瞬時の観測」を基本とする手法で
　わたくしが命名したものです。

（位相差測定・計算について）

　オシロスコープの発明と「リサージュ」図形の組み合わせによって
電気工学、電子工学では重要である「位相差」の測定ができるように
なりました。

　しかし、現在のデジタル化時代となっても、当時の方法から進化して
いないとは感じないでしょうか。　その証拠の一つに、２つの離散波形
から位相差を計算するという、単純なことに未だに苦労しています。

　当社で提供する「位相差計算」は、ずばり、数個の数式で構成される
単純なものです。Excelでも実行可能な平易なもので、電気・電子工学、
その他の分野で、「位相差計算」でお悩みの方の問題を解決します。

　従来のように、オシロスコープの目盛りを読んで計算する必要は
なく、波形データから直接位相差を計算できます。
アルゴリズムの性質上、1)周波数、2)振幅、3)位相差をワンセットで
求めることができます。

（周波数測定について）

　周波数カウンター、ソフトウェア的な方法（FFT,DFT）を使うことで
周波数測定は実用上問題ない技術水準に見えます。しかし、バイオリン
演奏のビブラートを解析できるだけの「精密周波数測定」は存在したで
しょうか。　ビブラートを測定するためには、一秒間に数十回の測定が
必要です。しかし、DFT/FFTによる方法では、時間分解能と周波数分解能
を同時に高めることはできません。

（提供する周波数測定技術の特徴）

　FFT/DFTは使用しません。したがって、フーリエ変換の「不確定性」
の影響をうけず、高い時間分解能と周波数分解能が可能です。

・　周波数の微小な変化の測定なら μHz、位相差での精密測定は μrad 級
　　のものを提供しております。

（使用機器）

　離散信号処理ではAD変換器を使用しますが、サンプリング周波数の誤差を
意識されたことがはあるでしょうか。サンプリングクロックは水晶発振器に
よって生成されているため、電子部品の等級にもよりますが、ppmオーダーの
偏差・変動は避けられません。

　当社で提供する技術はppm 級の不確かさを扱うものです。したがって、
AD変換器のサンプリングクロック周波数も問題となるケースが存在します、
1μHz級をめざす場合は、AD変換器のクロックはppb級に校正しておく必要が
あり、ルビジウム発振器で安定化されたWORDクロックを使用するAD変換期
を推奨しております。

・　通常の条件なら、使用機器は、普通のPCと、ちょっとだけ高級なオーディオ用
　　AD変換器をご用意くだされば大丈夫です。
・　余裕のある方は高精度のWORDクロック発振器で同期されたインタフェース
　　をご使用くだされば、安定性の高いデータが測定可能です。

（使用している機材）

• ADC             Steinberg社製 UR-824 (WORDクロック入力有り） 
• WORD Clock 城下工業社製　SWD-CL10 （TCXO級水晶を使用）
• Rb OSC        Symmetricom社　LPRO-101 （中古を2010年に入手）

（周波数測定の原理）

　基本演算子に「リサージュ外積」と命名した演算を用いています。
サイン波に対してリサージュ外積演算を行うと、波形の任意の場所で、この
リサージュ外積は一定値をとるという性質があります。

　また、条件の異なる、２種類のリサージュ外積を計算し、その「比」を
求めると、この「リサージュ外積比」は、"周波数のみの関数"となり、
振幅項は関与しません。この原理を利用して精密周波数測定を実現しています。

ーーー
　当研究所の最新成果、現場で使えるさまざまな「信号処理テクニック」を、
随時紹介していきたいと思います。

（精密周波数測定の例）

まずは、1/100 Hz 程度の「周波数偏差」を「実時間測定する技術」。

・　振り子の先に設置したスピーカからの音を床面のステレオマイクでキャッチします。
・　音源とマイクとの相対速度で、１Hz以下のドップラーシフが発生します。
・　検出するドップラーシフトによる周波数偏差は、1/100 ～1/10 Hz程度です
・　測定は 23回／秒と、非常に短時間で測定できます。

詳しくは,　この動画！
https://www.youtube.com/watch?v=59Fk8jCeU84

解説した論文は、ここ！。
http://www.eurasip.org/Proceedings/Eusipco/Eusipco2013/papers/1569746251.pdf