簡易式の位相差計算方法で、どこまで「小さな位相角」がはかれるか試してみました。
つい先日、発表されたブラックホールの視角は 1マイクロ秒角。 1度の 1/3600 さらに 100万分の1というオーダーです。
本Blogで紹介した方法で計算してみると、 「１マイクロ秒角」で有効数字３桁程度で計算できていました。
２つの信号の「位相差」を極小の状態でも計算できるので、そのうち・・・だれか効果的なアプリケｰションを開発してくれることを夢見ています。現実問題は、いかにPureな正弦波をとりだすか。または、含まれるノイズの影響をどこまで取り除けるかだろうと思います。
フーリエ変換使用すると、「成分」を分離することができます。 Synchro PRIMOでは、少ないサンプル数で「単一成分」を精密に計算できますが、「複数成分」の同時計測は可能なのでしょうか。
以前、２成分でシミュレーションを行い理論式も導出してみましたが、意味がさっぱりわかりません。しかし最近になってひとつ重要なことに気づきました。
本ブログで紹介しましたが、「２成分の混合波形」は、それぞれの周波数、振幅によって定義される「瞬時周波数」の振る舞いを呈します。この計算式は L.O.コーエンの文献にも解説されていることを紹介しました。 逆に言うと、２成分ということがわかっているなら、瞬時周波数の振る舞いから、それぞれの成分の周波数・振幅を計算できることになります。それぞれの成分の振幅比が大きい場合（＝どちらかの成分が優位）瞬時周波数は脈動し、脈動する周期の逆数が「周波数差」として、中心周波数からの変動量で、おそらく・・・振幅が計算できるはずです。
では３成分以上になるとどうか・・・数学的には、未知変数の数＋１のデータが要求されます。成分数が無限になると・・・必要なデータ量も増えていきます。実現のためには「時間窓」を大きくし、瞬時性を犠牲にするか、さらに高精度のサンプリングを行うしか方法はありません。
ここで現実的な実装を考えてみます。サンプリング周波数を高くしてみましょう。隣接サンプルとの値の差がどんどんなくなり、振幅方向のさらに高い分解能が必要になります。振幅方向の分解能が高いと、今度はノイズが問題になってきます。
以上をまとめると、成分数（＝帯域？）、瞬時性、ノイズの３つの制約が存在し、「超えられない壁」が存在するような予感がします。
実は「時間・周波数・ノイズ」の制約は、なにかシャノンの定理に絡んでいるな・・・と博士課程の段階で気づいていました。シャノンの「通信路容量」は、C = B* log(1+SNR) となります（＝シャノン・ハートレーの定理） 。
Synchro PRIMO でもシャノンの定理が見え隠れしているのです。（余談：この関連について博士論文の審査の際、別の先生からコメントがあったらしいです）
周波数は普通は０以上の実数です。しかし周波数変換・変復調をおこなうと計算上「負の周波数」が登場します。
では、複素数になったらどうなるか？ 一応「複素周波数」というものはあり z = ρ + j*ωt のように表現できます。実際の波形には１）振動項、２）減衰項が存在します。
波形を、f(t) = 1/2 { exp(z) - exp(z) } と考えると振動項・減衰項もあらわれ、なんとなく見えてきます。ただし、この式のzをどう呼べばいいのでしょうか。とりあえず「複素位相角」としておきます。
Synchro PRIMO では、指数的に減衰しながら振動する波形の周波数は正しく求めることができます。上記の複素周波数がうまく処理できているからです。しかし、一般の振幅変調波形ですと、一定の値が計算できません。さあここからどうするか。
振幅変調された信号の基本周波数は一定なのでしょうか？ ここは直感的な理解と複素解析から得られる結果が一致しないことも覚悟して検討しなくてはなりません。 上記の「複素位相位相角 z」 の微分が「複素周波数」であるとすればよさそうです。そう仮定すると、σ＝一定の減衰をともなう振動では、瞬時周波数の計算には影響しないことがなんとなく分かります。
振幅変調波形の瞬時周波数を複素解析的にもとめる前に複素平面で位相角の動きをイメージしてみます。 複素位相 z は普通は虚数軸方向に「速度ωで走って」いますが、振幅変調がはいると、実数軸方向に「ぶれる」ことになり、瞬時的な速度＝瞬時周波数は「一定でない」と考えられます。ここが盲点で、我々が無線工学でイメージしている周波数と変調された信号の瞬時周波数は厳密には異なるということが示唆されているのです。
Synchro PRIMOの計測（計算）では、指数減衰する三角関数の周波数は厳密に求めることができます。（指数形式で入力信号を与えると、計算式のLsを展開することで容易に証明できます）
では正弦波で振幅変調された波形の瞬時周波数はどうなるか。log (1+m*sin(ωt)) (mは変調指数）に比例した変動ではないかと考えています。
２つの、異なる周波数成分を足し合わせても、「平均周波数」は振幅の大きい側に支配されることは以前述べました。今度は、成分が無限に含まれているとされる「雑音」について、周波数はどう表現できるかを考えてみたいと思います。
私がいま考えている方法は、現波形と１サンプルずらした波形をリサージュ図に描画したとし、その「回転数」をみる方法です。厳密には確率的信号処理の考えで検討しなければならないのですが、ノイズのスペクトルの包絡線の違いでさまざまな計算値がでると思います。
複数の周波数成分が「入り乱れている」信号の周波数領域の解析において伝統的にF0が重要視されますが、調波の構造が比較的単純でも「ミッシング・ファンフダメンタル」のような現象がおこると、高調波をすべて計測し、その構造からF0を推定するという厄介な作業となります。
今回提唱している「平均周波数」（名前を工夫して、代表周波数みたいなものも可能）の定義と数理が明確になれば、いちいちフーリエ変換しなくても、かなりのことがわかるのではないかと思います。
フーリエ変換は便利な手法であるため、広い領域で使われています。時間の区間が±無限であれば）任意の波形は成分に分けられる（＝順変換）ということはたぶん成立するのでしょうが、逆に、任意の成分を用意すれば任意の波形は合成できる（逆変換）は 微分不能な波形を生成することがあります（リーマン関数など）。連続波形なおのに導関数が定まらないというケースです。 どう解釈してよいのか今の私には分かりません。
私見ではありますが、フーリエ変換の原点にかえって「正しい使い方」をもっと知ってもらいたいなと思っています。ＦＦＴが考案されて以降、濫用されているような気がしてなりません。
「連続スペクトル」「スペクトル包絡」などという言葉がありますが、さらに意味を考えてみます。
連続系の信号で、明らかに連続スペクトルになるものには、一例として、「インパルス」「ガウス関数（ガウシアン）」があります。インパルスの時間幅はゼロ・振幅は無限大・積分すると１という定義はともあれ、現実的には、一定の時間幅をもった矩形波で考えてもいいでしょう。
では、インパルスは平坦な連続スペクトルなので、周波数がすこしずつ異なるサイン波を無限に重ねればインパルスを合成できそうです。。次の思考実験ですが、パルスをつかって三角関数は合成できるでしょうか。離散系では、sin関数であっても、δ関数の和として表現できます。ここで注意ですが、δ関数のスペクトルは平坦。そんなδ関数をつかって合成すると、１本の線スペクトルになります。この理由は、線スペクトルといっても、各周波数で位相角はきまっていて、位相まで含めて合計すると１本の線になってしまうということです。δ関数の連続スペクトルを合成するとさまざななスペクトルが表現可能、こういう見方も可能です。
今度はガウス関数。フーリエ変換した結果もガウス関数という、面白い性質をもっています。 ガウス関数をつかって任意の波形が合成できるなら、フーリエ変換も同じような重ね合わせとなってしまうでしょう。この方法は・・・考え中です。
現時点で悟っていることは、「厳密な連続スペクトルはパルスでしか得られない」、「線スペクトルはいくつ足し合わせても線スペクトルの塊」・・・有限和と無限和の境目がむずかしいです。
ちょっとした思考実験です。 フーリエ級数（変換）のように、任意の波形が「基底波形」の無限和で表現できるのはご存じのとおりですが、無限積ではきないか。
まずは有限和で一例を。 f(t) = sin(ω1 t), g(t) = sin(ω2 t ) として y = f + g を考えてみます。三角関数の和積変換を用いると、 y = 2 * sin ( (ω1+ω2)/2) * cos ((ω1-ω2)/2) と表現できます。 ２成分で振幅が等しい場合は、和を積の形に変形できます。
では、振幅が異なる場合；f(t) = sin(ω1 t), g(t) = 2*sin(ω2 t ) ならどうするか？
y = (f + (1/2 g)) + (1/2)g とすると、 y = 2 * sin ( (ω1+ω2)/2) * cos ((ω1-ω2)/2) +sin(ω2 t ) となります。振幅変調のかかった成分と、ω２成分の残骸。 さらにこの２つを統合できないものでしょうか。
次のケース。振幅はすべて等しいとして、３成分の合成ではどうなるか。 y = sin A + sin B + sin C のようなケース。高校でならった和積の３成分への拡張です。実は数学オリンピッククラスでは当たり前の公式だそうです。 https://mathtrain.jp/waseki ただし A+B+C = π の条件があります。この制約を取っ払うのも一苦労。３成分ができれば、Ｎ成分への拡張も・・・根性いれればできるでしょう。
あとは、振幅が違い周波数が同じ場合ですが、合成公式：https://mathtrain.jp/asinbcos を使えば１つに集約可能。振幅が違っても、うまく合成公式を活用して、すべての成分の振幅をそろえて、Ｎ成分の「和積」をつかえば目的のものは完成です。 その中ででてくる各項の周波数は何を意味しているのか。 我々が知ってる級数の和と同じことを表しているわけですから、その解釈に、いままで築かなかった「意味」があるのかもしれません。
今時点ではっきりしているのは、振幅の等しい２つの成分の和は、べつの２成分の積で表すことができること。ここをとっかかりに、なにかをブレークスルーしたいと思っています。
本ブログで提案している「瞬時周波数」の計算公式ですが、入力に正弦波を加えると波形の任意の場所で周波数が計算可能です。
この公式にパルス波形を入力するとどうなるか？ 例として 「ガウス関数」を計算してみましょう。あれれ・・・ 計算できてしまいます。しかも一定の値を示します。もちろん、時間＝ ±∞ではほとんど計算不能ですから、パルスが立ち上がっている箇所がどうなるかが重要です。
ガウス関数 は g(t) = exp ( - A・t^2) という形で、とりあえず定義できます。もともとは正規分布で活躍する関数なので統計量を意識した定義方法がありますが、ここでは、パルスの急峻さ=A のみにで考えることにします。
https://ja.wikipedia.org/wiki/ガウス関数
もし、瞬時周波数の計算結果が正しいならば、ガウス関数によるパルス波形の周波数が仮想的に定義できることになります。ちょっと計算してみれば分かるのですが、波形がなだらかなな場合は F=0に近づき、急峻な場合は、どうやら Ω=0.16666 付近に収束しているようにも思えます。今のところ意味はよくわかりません。
ちなみに、ガウス関数を ±∞の区間で積分すると（＝ガウス積分）・・・・ガウス関数の不定積分は存在しないのに、無限大区間での積分値は計算できてしまいます。
https://ja.wikipedia.org/wiki/ガウス積分
もう少し妄想をするなら、「任意のパルス波形」はガウス関数を用いて合成できるのか？ ということになりますが、いくつか制約を加えれば意味のあるものが可能なのかもしれません。また周波数が計算可能なら位相角も定義可能であり、同時に、時間軸との対応もピタリとあうはずです。パルスの瞬時値・瞬時周波数に・・・なにが隠されているか・・・・
Synchro PRIMOのアイディア「リサージュ外積」の計算に現れる数学を紹介します。
別に公式というほどでもないのですが、解法には大学受験で覚えた三角関数の公式が登場します。
（余弦自乗差）・・・と勝手に命名.
cos^2 nθ - cos^2 (n+1)θ = sinθ * sin (2n+1)θ
例として n=0 の場合、 1 - cos^2 θ = sinθ * sin θ となり正しいことが分かります。リサージュ外積 Ls = Xn-1 * Xn+1 - Xn-2 * Xn+2 をｚ領域にこじつけると、上の式が登場します。 （信号の積を使う非線形関数にｚ変換を使ってよいものか自信ないですが）
右辺の正弦積は、１）自乗差を "和差の積" に分解、２）それぞれを"和積"の公式で展開。３）項を入れ替え倍角の公式を再適用 という手順求めることができます。 Synchro PRIMOは、「リサージュ図の面積はどうなるか」からスタートしていますが、実は 200年ほどまえ Pronyという数学者が同じようなことをやっています。下記はProny法をのべている数少ない論文です。
https://www.jstage.jst.go.jp/article/jasj/42/11/42_KJ00001454017/_article/-char/ja/
Pronyの原著らしきものも入手したのですが、悲しいかなフランス語が読めないのと当時の数式の流儀が分からず手が出せていません。。。
Synchro PRIMO法では、複素数を用いた連立方程式を回避し、離散波形の隣接５点の値からProny法とそっくりな計算式を導出しているのが特徴となっています。その基本となっているのが、「余弦自乗差」なのです。数IIIの演習にも使えそうな式ですね。
１０ヶ月ほど前に下記：
How to calculate Frequency, Amplitude and Phase Difference from Discrete Data
という動画をを投稿しています。最近になってアクセスが微妙に増えてきました。なぜかインドと米国から集中。アクセスの分析はクリエータツールからできます。
動画の内容なのですが、２つの波形（正弦波）を定義したExcelシートからスタートして、１）それぞれの周波数の計算手順を示しながら瞬時周波数を算出、２）同様に瞬時振幅を計算、３）それぞれの信号を複素化して、４）位相差を計算 という手順のExcel操作をデモしています。計算式は国際会議で発表ずみのものです。
計算式の数学的意味はさておき、周波数・振幅・位相差がピタリとシートに出現するシーンは、今改めて見直してみると面白いです。（最近、細かいことを忘れかけているので、新鮮な目線で見ることができている？） 特にオリジナルのアイディアである「リサージュ外積」の計算はExcelのセル選択の様子で一目瞭然です。
画面を見ながら（全画面表示なら文字はくっきり）同じようにセルを選択・計算式を入力すると、全く同じことができますので興味あるかたはチャレンジしてみてください。
今まで、数機種の「オーディオインタフェース」（ＡＤＣ）を使用してきたのですが、概してサンプリング周波数が低めなのです。
以前からローランド製 UA-22 を携行用に使用。この機種はＵＳＢからの電源供給でコンパクトです。サンプリング周波数(基準 44.1kHz）をGPSの１秒パルスを使用して精密測定すると、 -38 ppm 。電源を入れて数時間はドリフトがありますが最後は安定します。
最近購入した UA-22の後継機種： Rubix-22 、これも同様に実験したところ、Fs = 44091.48 との結果がでました。つまり１秒間あたり、8サンプル分遅いのです。水晶発振器のクラスには SPXO, TCXO, OCXO ・・・さまざまあり、安価なものですと 100ppmくらいずれても普通。もしそれが問題なら、「外部同期」付きの機種を使用し、Word Clockで同期すればいいだけです。とはいえ、 50ppm程度におさまっているかと期待していたのが、今回の結果( -190 ppm )はちょっと納得できません。
とはいうもの、周波数安定度は短期・中期とも抜群で、さすが音楽用という印象です。
現在、電力用PMUを安価なＨＷを使用した構成を検討していますが、時刻の基準となるGPSからの 1pps信号とADCのサンプリング周波数の補正をどうするか悩んでいます。GPS 1pps信号受信機は2000円くらいで購入（秋月電子）できるので、これを利用しサンプリングクロック周波数を実時間で追跡できると使用するADCのコスト削減が可能かと思います。以前は 高級機種に「ワードクロック発振器」（精度は0.3ppm級）にくわえルビジウム発振器で外部同期していましたが・・・機材の値段、ゼロが一つ多いのです。
広域に接地する機器の「時刻同期」「周波数同期」は思いの外難しいです。
ちなみに、今回の精密測定では、１年前に考案した「パルスの精密時間差」アルゴリズムを使用しています。Synchro PRIMOで使用しているアイディアをPulseに拡張したものですが、数学が超難解で悶絶中です・・・１本の積分計算式の導出にレポート用紙４枚。結果はそのうちに。
