Synchro PRIMO 法は、簡単な方程式を用いて、５点のサンプル値から瞬時周波数・瞬時振幅を計算するアルゴリズムであり数式です。数学的な意味は、三角関数 Y[n] = A*sin (Omega * n ) のYn を５つ用意すれば、周波数であるΩ,振幅A を求めることができるというものです
さて、多くの方から質問・コメントをいただいているのですが、精度向上の方法の考え方を説明しようと思います。
Pythonのサンプルcode 付きの記事を投稿していますので、ご参照ください。
（１） 瞬時周波数
入力信号を x[n] としたとき、Ls1,Ls3 は下記のように「サンプル間隔=k」を「広げる」ことで、より安定な計算ができます。
- Ls1[n] = x[n]*x[n] - x[n-k]*x[n+k]
- Ls3[n] = x[n-k]*x[n+k] - x[n-2k]*x[n+2k]
この場合、最終的な周波数計算値は
freq = fs / (2*PI* k ) *arcsin( sqrt( 3 - Ls3 / Ls1 )) となります。逆SINのあと k で補正します。
（２） 平均化
理想的なサイン波であれば、Ls1, Ls3は定数になります。そのため、 1/2 sqrt ( 3 - Ls3/Ls1 ) の計算は非常に安定となります。
ノイズを含む波形の場合は、Ls1,Ls3を、XY平面にプロットすると、ノイズである程度拡散されたエリアにプロットが現れます。計算式中にあらわれる Ls3/Ls1 は原点からプロットまでの傾きとなるので、プロットを平均化により一点に集約することで解も一点に収束します。
（３）ノイズ耐性
計算式を見てわかるかと思いますが、 Ls3 =0 の近傍が、最終的な計算結果に及ぼす影響が最小になるのは明らかです。 Ls3=0 のとき、 F= 0.16666 (=1/6) となりますが、逆にこの性質を応用し、前述の「サンプル間隔」を人為的にチューニングし、F=0.16666 あたりを利用するすることも可能です。
（４）Ls1, Ls3 がノイズにより受ける影響
適当な正弦波波形にノイズを加えた信号を用意し、Ls1,Ls3を計算してみます。結果を、XY平面にプロットすると「ある一定の範囲」に収まることが確認できます。
このとき、原点からP=（Ls1,Ls3）を結ぶ直線の傾きRが瞬時周波数値を決定してしまいます。つまり、Rの値をいかに安定に求めるかが実装的なテクニックになるわけです。
計算式からわかる通り、周波数を実数の範囲で計算するためには、 -1 ≦ R ≦ 3 が条件となります。この範囲に、前述のプロットがおさまらないと、計算エラーを発生します。しかし、Ls1, Ls3 を平均化することで、 R のとる範囲はぐっと小さくなります。興味のある方はシミュレーションしてみてください。
（５）外れ値の解釈
Ls1, Ls3 の計算値によっては、 R < -1 となったり、３＜R となったりします。この時の数学的意味については複素関数の知識が必要です。
（まとめ）
・ サンプリングインターバルを適切に設定することにより、R値が変化することで計算は安定する。
・ Ls3 =0 近傍になるように条件設定すると、ノイズに強い
・ Ls1, Ls3 を平均化してから、 sqrt ( 3 - Ls3 /Ls1 ) を計算するとばらつきが減少する