２か月間　BNO085のyaw角誤差をRTK　GPSの高精度パラメータを使って補正する作業をしてきましたが、
最終的に、RTKで大きなオフセット補正を合わせてから、波形内の誤差の最大の外れ値の周波数が4-5Hzである点に注目して、１Hzのローパスフィルタで、元波形をつぶすことで、誤差のレンジを半減近くにすることで終わりとしました。
あとは、２０２５シーズンで、いろいろな滑走データと比較検証して、妥当な処理か検証します。

●誤差外れ値の波形

代表的なターンを切り取って、フィルター有無で比較してみると　黒が元波形、赤がフィルタ後波形
上のグラフがHeading角とyaw角のプロットで、２００度から８０度まで回転してます。CCWなので左ターンです。
下のグラフがheading角とyaw角の差分＝誤差をプロットしてあります。
横軸は、時間で、１目盛10msecです。下グラフの黒の波形をみると25-40目盛周期で大きな誤差が発生してます。
これは、Heading角が揺れているのですが、IMUのyaw角はこの揺れを感知してませんので、誤差となってでてきてます。どちらが正しいのかというとHeading角が正しいはずですが、RTKの周期が120msecと遅いため、周期変動が
２－３サンプリングで７度も揺れていることを示してますが、スキー板で方向がそんなに速い周期で変わるのも不自然なので、この大きな外れ値は、基準となるRTK GPS側のノイズだろうと判断して基準波形にローパスフィルタを１Hzカットオフでいれました。そうすると黒線が赤線までなだらかになって、誤差が減少するという結果になります。

●全体の統計値
抜けが多いのですが、抜け以外のデータをみると５０－８０％レンジまで縮小してます。
角度範囲でいうと、だいたい±３．５度以内におさまっているものが多いということで、±５－８度だった元波形を
フィルターでつぶして、誤差を圧縮したという結果です。

　●ローパスフィルタのやり方はscipyライブラリを使えばバターワースフィルタが簡単にできます。
コードはgistにあります。https://gist.github.com/dj1711572002/f781fc4eb451b36e230616341e6ffc39

シンプルな説明が判り易かったので、こちらのサンプルプログラムをいただきました。作者様に感謝します。

【Python】scipy.signal.butterでローパス（ハイパス）フィルタを作成する

コードは。def LFP(df,coln)にdfとローパスフィルタをかけるカラム名を引数として渡します。
戻り値は、dfに新たに、coln+”_f”とサフィックスがついたフィルター済みカラムが追加されて戻りdfとして受け取れます。フィルタの詳細はわからないので、いろいろな参考サイトがあります。
https://watlab-blog.com/2019/04/30/scipy-lowpass/

from scipy import signal

def LPF(df,coln):

print(“df_colnにLowPassFilter処理して、df_coln_fというコラムをつけて戻す”)

#————butter worse————————–

# パラメータ設定

sample_freq = 100

cutoff_freq = 1#1hz固定

filter_order = 4

sec = len(df)*0.01#dfの時間

# —-colnをフィルタ

t = np.linspace(0, sec , int(sample_freq * sec))

y=df.loc[:,coln].to_numpy()

sos = signal.butter(filter_order, cutoff_freq, ‘lowpass’, output=’sos’, fs=sample_freq)

yf = signal.sosfiltfilt(sos, y)

print(yf)

new_coln=coln+”_f”

df[new_coln]=pd.DataFrame(yf)

return df

これを使ったlpf処理の機能
ターン毎に分割されたcsvデータファイル群を一個ずつ読んでフィルター処理して、プロット新たにデータとグラフと統計値保存する機能です。

elif v==’lpf’:#既存のファイル群をまとめてフィルター処理を追加

print(“low pass filter “)

print(“>>>> SELECT sITP_ _ FILE “)

#headMOtのピークを真の値としてHeading360とyawのピーク高さを合わせることで補正とする

root = tkinter.Tk()

root.withdraw()

iDir = r’C:/RTK_Log/hokan/’

fTyp = [(“BNOデータファイル”, “*.csv;*.xlsx;*.xls”), (“すべてのファイル”, “*.*”)]

filename = FileDialog.askopenfilename(parent=root, initialdir=iDir, filetypes=fTyp)

basenameITP = os.path.basename(filename)

#シリーズファイル群まとめて処理

fsname=csv_Allname(filename)#ファイルシリーズのフル名のソート済みリスト

for i in range (0,len(fsname)):

print(i,”:”,fsname[i])

#==============LPG処理======================================

i=0

list_stat=[]

list_statAll=[]

while i<len(fsname):

print(“====>”,i,fsname[i])

fn=fsname[i]

basename = os.path.basename(fn)

df = pd.read_csv(fn, low_memory=True)

df=LPF(df,”Heading360″)#LowpassFilter関数

df=LPF(df,”yawHpg”)#LowpassFilter関数

print(df.columns)

print(“df[yasHg_f]=”,df[“yawHpg_f”])

df[“H-Ysapg_f”]=df[“Heading360_f”]-df[“yawHpg_f”]

df.to_csv(r”C:/RTK_Log/hokan/LPF/LPF_LOG_”+str(i)+”_”+basename+”.csv”, index=True) # CSV sav

#max min range

old_max=df[“H-Ysapg”].max()

old_min=df[“H-Ysapg”].min()

old_ave=df[“H-Ysapg”].mean()

old_std=df[“H-Ysapg”].std()

old_range=old_max-old_min

fltd_max=df[“H-Ysapg_f”].max()

fltd_min=df[“H-Ysapg_f”].min()

fltd_ave=df[“H-Ysapg_f”].mean()

fltd_std=df[“H-Ysapg_f”].std()

fltd_range=fltd_max-fltd_min

list_stat=[old_max,old_min,old_ave,old_std,old_range,fltd_max,fltd_min,fltd_ave,fltd_std,fltd_range]

list_statAll.append(list_stat)

#plot準備

fig = plt.figure()

x=np.arange(len(df))

yHD0=df[“Heading360”].to_numpy()

yHDf=df[“Heading360_f”].to_numpy()

yyw0=df[“yawHpg”].to_numpy()

yywf=df[“yawHpg_f”].to_numpy()

yHY0=df[“H-Ysapg”].to_numpy()

yHYf=df[“H-Ysapg_f”].to_numpy()

#plot

ax1 = fig.add_subplot(2,1, 1)

fig.suptitle(“LPF_1Hz_”+fsname[i])

ax1.set_title(“Heading,yaw plot”)

ax1.plot(x,yHD0,label=”Heading360″,color=”orange”)

ax1.plot(x,yHDf,label=”Filtered Heading360 “,color=”red”)

ax1.plot(x,yyw0,label=”yawHpg”,color=”blue”)

ax1.plot(x,yywf,label=”Filtered yawHpg”,color=”green”)

plt.legend()

ax2 = fig.add_subplot(2,1, 2)

ax2.set_title(“Error=Heading-yaw plot”)

txt1=”old:max=”+str(round(old_max,1))+”,min=”+str(round(old_min,1))

txt2=”Filtered:max=”+str(round(fltd_max,1))+”,min=”+str(round(fltd_min,1))

ax2.text(x[0],old_max, txt1, size=10,color=”black”)

ax2.text(x[0],old_max-1,txt2, size=10,color=”red”)

ax2.plot(x,yHY0,label=”source Error data”,color=”black”)

ax2.plot(x,yHYf,label=”Filtered Error data”,color=”red”)

#ax2.plot(t,yHY0,label=”HeadingFiltered data”,color=”green”)

plt.legend()

plt.show()

fig.savefig(“C:/RTK_Log/hokan/LPF/PLOT_LPF1Hz_”+basename+”.png”)

i+=1

df_stat=pd.DataFrame(list_statAll)

df_stat.columns=[“old_max”,”old_min”,”old_ave”,”old_std”,”old_range”,”fltd_max”,”fltd_min”,”fltd_ave”,”fltd_std”,”fltd_range”]

df_stat.to_csv(r”C:/RTK_Log/hokan/LPF/LPF_STAT_0-“+str(i)+”_”+basename+”.csv”, index=True) # CSV sav

print()

●以後
IMUのyaw角補正は、

０段目：ターンの開始から終わりまでの範囲で処理を行う。
1段目：RTK headMot角のピーク値とyaw角の差分のオフセット補正する
２段目：yaw角を１Hzのローパスフィルタを通す
以上で、誤差が±３度付近まで減少すると思われます。
以降、この補正方法を基準に、
いろいろな滑り方の実験をして、統計的な誤差がどの程度おさまるか検証していきます。