瞬時周波数を求める方法（その２）を紹介します。
５サンプルの離散波形による瞬時周波数計算 の変形です。
”２倍角法”と命名しました。
１） 周波数を計算したい波形から、4サンプルをとりだします。
１周期がちょうど１００サンプルとなる波形を用意しておいて、
任意の位置から4サンプルとりだしてみましょう。
もし、サンプリング周波数が 44100 Hzとするとこの波形の実周波数は441Hzになります。
時系列波形は；
x[n] = sin ( ２* PI * 441 * (n/44100) ) と計算できます。
位相差項をいれても構いません。
２） とりだしたサンプルルを、順番に x1, x2, x3, x4, x5 とします。
３） 以下のような Ls1 と Ls2を計算します。
Ls1 = x3*x3 - x2*x4
Ls2 = x2*x3 - x1*x4
単純な一種の「たすき掛け演算」です。
４） つぎに、Ls1, Ls2 から、 r = Ls2 / Ls1 のように
比をとります。「リサージュ外積比」と呼んでいます。
５）以上のように求めた r から次の計算で周波数（瞬時周波数）が
計算できます。
f = (fs / 2 PI) * arccos( 0.5 * r)
fs はサンプリング周波数です。 ”２倍角法”を命名しました。
x5 までの５点をつかう方法も可能です。 ”３倍角法”と命名しました。
リサージュ外積：Ls1, Ls3を以下のように計算し、
Ls1 = x3*x3 - x2*x4
Ls2 = x2*x4 - x1*x5
f = (fs / 2 PI) * arcsin( 0.5 * sqrt(3-Ls3/Ls1))
となります。
リサージュ外積:Ls1,Ls2,Ls3は限られたサンプルから定義にしたがった計算で機械的に求めることができます。あとは２つのリサージュ外積の比をとって、上記の逆三角関数に代入するだけです。
リサージュ外積 Ls_n を任意によった時どうなるかですが、 リサージュ外積 Ls_n, Ls_m (n>m. m≧0)を与えたとき、
Ls_n / Ls_m = sin((2n-1)Ω) / sin((2m-1)Ω)
となることが分かっています。この方程式を解いて目的の周波数Ωを計算するわけです。倍角の公式を使って変形することで、説明した公式が得られます。２倍角法の例では n=2 , m=1 とする例です。
さらに解法の一般化を試みると、
リサージュ外積比は x=sinΩ または x=cosΩとしたときの、チェビシェフ級数 Tn(x), 第二種チェビシェフ級数 Un(x) で表現できそうです。
現在のところ、「２倍角法」「３倍角法」のみわかっています。もっと高次なものでのやりかたはよくわかっていません。