「５点の直から周波数を求める計算式」を以前から紹介しております。
この方法は Synchro PRIMO法と、筆者は呼んでおります。
今回、python によるデモを公開します。波形を作成し、Synchro PRIMOで周波数を計算、その結果をグラフ表示するだけのものです。計算プロセスをソースから読み取っていただければ幸いです。
計算結果がグラフに表示されますが、非常に小さな「演算のゆらぎ」が見えますが、これは浮動小数点計算で誤差が蓄積した結果です。サンプルでは、周波数＝ 50.123456789 Hz と設定していますが、計算結果をみると小数点以下9桁までぎりぎり正しいようです。
６７行目の interval を小さくすると、計算するデータ幅は小さくなりますが、精度は落ちます。
ソース中の、 freq = 50.123456789 # Hz を修正し、他の周波数も試してください。周波数が低い時のみ、interval 値は高めにセットし、結果がどうなるか試してみください。
実行すると、波形が表示されます。このウィンドウを閉じると計算を再開し、次に計算結果を表示します。
#========================
import numpy as np
import random
from matplotlib import pyplot as plt
PI2 = np.pi * 2.0
#
# Generate Wave
#
fs = 44100 # sampling freq.
freq = 50.123456789 # Hz
init_phase = 10 # deg.
NumSample = 1000
N = range(0,NumSample)
FREQ = freq / fs
Omega = PI2 * FREQ
Phi = init_phase/ 180.0 * np.pi
# amplitude
Amp_sin = 1.0
# sinusoid
sineWave = np.array( [Amp_sin * np.sin(Omega * n - Phi) for n in N ])
# show wave
waveForm = sineWave
plt.plot(waveForm)
plt.show()
#
# Synchro PRIMO
#
def synchro_primo(x, pos, k=1):
# samples
x0 = x[pos - 2*k]
x1 = x[pos - 1*k]
x2 = x[pos]
x3 = x[pos + 1*k]
x4 = x[pos + 2*k]
# Lissajous Products
Ls1 = x2 * x2 - x1 * x3
Ls3 = x1 * x3 - x0 * x4
R = Ls3 / Ls1
sinOmega = 0.5 * np.sqrt( 3 - R )
Omega = np.arcsin( sinOmega )
FREQ = Omega / PI2 / k
freq = FREQ * fs
return freq
#
# calculate Instantaneous frequency
#
Num = 500
Interval = 21 # sampling interval
# output array
calc_freq = []
calc_amp = []
for pos in range(2 * Interval, Num - 4 * Interval):
freq = synchro_primo(waveForm, pos, Interval)
calc_freq.append(freq)
# show Average
print( np.average(calc_freq))
# graph (freq)
plt.plot(calc_freq)
plt.show()