はじめに

以前書いたnote記事「動かして学ぶバイオメカニクス#5　〜フィルタリングと加速度計算〜」の平滑化スプライン（smoothing spline）のコードを見直しました．

平滑化スプラインはオフラインで計算を行いますので，最適化を用いて計算するという意味では，リアリタイム処理には使用できません．しかし筆者は，バイオメカニクスなどの力学計算では，最も優れたフィルタで，力やトルクなどの計算における二階微分による力学計算では，計算困難な理由がない限りできるだけ平滑スプラインを使用すべきと考えています．

平滑化スプラインは，観測データをスプライン関数（spline function）で近似します．ただし，これは「オリジナルのデータを必ず通る」補間（interpolation）ではありません．スプライン関数による関数近似は行いますが，曲線を二階微分しても滑らかになるような４次の関数を選択し，かつオリジナルのデータの「近傍」を通る点に最適化することで，平滑化と微分を同時に実現することができます．微分は差分ではなく，スプライン関数という数学関数で記述していますので，厳密に微分値を得ることができます．つまり，もとのデータとその微分データは，完全に数学的な微分関係が成立しています．一方，差分計算は近似ですから誤差を含みます．

さて，ここではバイオメカニクスでも一部の人が使用しており，比較的性能がよいと言われている，Singular Spectrum Analysis (SSA)という特異値分解（SVD）を用いたフィルタリング（特異スペクトル解析）との比較も行いました．この記事ではSSAについては簡単な説明しか行っていませんが，前回の記事「Numerical Recipes in Biomechanics #1 ~回転行列のキャリブレーション（平均化）~」でも取り上げた特異値分解（SVD）が理解できれば，フィルタリングの意味がわかるでしょう．

結論を述べると，SSAはバターワースなどのフィルタよりは性能は良さそうですが，多くの部分で平滑化スプラインのほうが優れていると思います．ただし，平滑化スプラインも，使いづらいところもありますので，使い分けるとよいでしょう．

平滑化スプライン

以前の記事で紹介したコードと機能的にはほぼ違いはないですが，コードをAIを頼りにかなり書き直しています．若干コンパクトなコードになりました．

加えた機能としてはデータが１次元の配列でも，ベクトル（２次元配列）にも対応しています．データがベクトルの場合に，重みを配列で個別に与えられます．ノイズの性質がデータのチャンネルによって，異なることもあるからです．

なお，もともとの関数でも欠損値（NAN）に対応していましたが，引き継いでいます．

あまり詳しく解説しませんが，使用例などを参考にしてください．

def smooth_spline(data, weights=1.0, sf=1000.0, order=0, degree=4):
    """
    UnivariateSplineを用いて時系列データを平滑化（および微分）する。

    Parameters
    ----------
    data : array-like
        平滑化したい時系列データ。1次元または2次元配列
    weights : float or array-like 
        平滑化の重み。スカラーまたは配列．配列の場合データの列数とし，各列ごとに個別に重みを与えることができる
    sf : float
        サンプリング周波数（Hz）
    order : int, optional
        微分の次数（0なら平滑化のみ）。デフォルトは0
    degree : int, optional
        スプラインの次数。デフォルトは4．加速度まで計算する場合はデフォルトの4とする．

    Returns
    -------
    smoothed : ndarray
        平滑化（または微分）されたデータ
    """

    data = np.asarray(data)
    if data.ndim == 1:
        data = data[:, None]
    
    # weightがスカラーまたはデータの列数と一致する配列であることを確認
    if not np.isscalar(weights):
        assert len(weights) == data.shape[1], f"weights must be scalar or array of length {data.shape[1]}, but got length {len(weights)}"
    
    # data = np.asarray(data)
    if data.ndim == 1:
        data = data[:, None]

    n_samples, n_series = data.shape
    t = np.arange(n_samples) / sf

    # 出力配列
    result = np.full_like(data, np.nan, dtype=np.float64)

    for i in range(n_series):
        y = data[:, i]
        valid = ~np.isnan(y)
        if np.sum(valid) < degree + 1:
            continue  # 有効な点が少なすぎる場合はスキップ

        x_valid = t[valid]
        y_valid = y[valid]

        # 重み処理
        w_valid = np.full_like(y_valid, 1, dtype=np.float64)
        

        # スプラインフィッティング
        if np.isscalar(weights):
            spline = UnivariateSpline(x_valid, y_valid, w=w_valid, k=degree, s=weights)
        else:
            spline = UnivariateSpline(x_valid, y_valid, w=w_valid, k=degree, s=weights[i])
            
        result[:, i] = spline(t) if order == 0 else spline.derivative(order)(t)

    return result.squeeze()

'''計算例'''
import numpy as np
import matplotlib.pyplot as plt

# サンプルデータ生成：2つのsin波にノイズを加える（2チャンネルの時系列）
np.random.seed(0)
t = np.linspace(0, 10, 100)
sf = 10  # サンプリング周波数（Hz）

# 2系列の信号（列ごとに異なるsin波）
signal1 = np.sin(t) + 0.2 * np.random.randn(len(t))
signal2 = np.cos(t) + 0.2 * np.random.randn(len(t))
data = np.column_stack((signal1, signal2))  # shape: (100, 2)

# 一部NaNを挿入（欠損のある実データを模倣）
data[30:35, 0] = np.nan
data[50:55, 1] = np.nan

w = [2.5, .01]

# 関数の使用（平滑化）
smoothed = smooth_spline(data, weights=w, sf=sf, order=0, degree=4)

# 関数の使用（1階微分）
diff1 = smooth_spline(data, weights=w, sf=sf, order=1, degree=4)

# 可視化
fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True)

axes[0].plot(t, data[:, 0], 'k.', label='original ch0')
axes[0].plot(t, smoothed[:, 0], 'b-', label='smoothed ch0')
axes[0].plot(t, diff1[:, 0], 'r--', label='1st derivative ch0')
axes[0].legend()

axes[1].plot(t, data[:, 1], 'k.', label='original ch1')
axes[1].plot(t, smoothed[:, 1], 'b-', label='smoothed ch1')
axes[1].plot(t, diff1[:, 1], 'r--', label='1st derivative ch1')
axes[1].legend()

axes[1].set_xlabel("Time [s]")
axes[0].set_ylabel("Signal")
axes[1].set_ylabel("Signal")
plt.tight_layout()
plt.show()

図１：２種類のデータの比較（上図：重み2.5．下図：重み.01）

Singular Spectrum Analysis

バイオメカニクスでは，一部の人達が，Singular Spectrum Analysis (SSA)という特異値分解を用いたフィルタリング（特異スペクトル解析）を行っていますが，ここで平滑化スプラインはと比較します．

SSAの基本的なステップを簡単に説明します．前章「Numerical Recipes in Biomechanics #1 ~回転行列のキャリブレーション（平均化）~」の特異値分解を復習してください．

埋め込み（Embedding）
・長さ$N$の時系列データ$\boldsymbol{X}=[x_1,x_2,,…,x_N]$を使って、ラグ（window length）$L$ による遅延座標行列$\boldsymbol{X}$を作成．
　窓の幅$K$を定め，窓をずらしながら時系列データを切り出し，軌道行列 (trajectory matrix)と呼ばれるの行列を作ります．
・これは、Hankel（ハンケル）行列（対角が同じ値を持つ行列）といいます．
・$K=N-L+1$とすると
$$
\boldsymbol{X}=\begin{bmatrix}
x_1 & x_2 & \cdots & x_K\\
x_2 & x_3 & \cdots & x_{K+1}\\
\vdots & \vdots & \ddots & \vdots\\
x_L & x_{L+1} & \cdots & x_N
\end{bmatrix}
$$
特異値分解
・$\boldsymbol{X} = \boldsymbol{U}\boldsymbol{\Sigma} \boldsymbol{V}^T$と分解
・ここで，$\boldsymbol{\Sigma}$は特異値の行列，$\boldsymbol{U}, \boldsymbol{V}$は直交行列です．
・特異値の大きい方から並べると，トレンド・周期・ノイズを意味します．
再構成（Reconstruction）
・大きな特異値から選択した成分で行列を再構成し，対角平均（diagonal averaging）によって時系列に戻します．
・ウィンドウ長$L$は，周期より大きく，時系列長の1/2程度に設定するのが一般的のようです．

周波数ベースのフィルタリングと似ていますが，ノイズが周波数で定まるのではなく，特異値の小さい成分で構成されている点が異なり，位相遅れも発生しません．また，筋シナジー解析のように非負（>0）の拘束がない限り，これはシナジー解析です．そういう意味では，シナジー解析はローパスフィルタやスペクトル解析に過ぎません．

ここでは示しませんが，バターワースフィルタなどと比較し，比較的性能はよいという印象です．

from numpy.linalg import svd

def embed_series(series, L):
    """時系列をHankel行列に埋め込む"""
    N = len(series)
    K = N - L + 1
    X = np.column_stack([series[i:i+L] for i in range(K)])
    return X

def diagonal_averaging(X_recon):
    """再構成された行列から時系列への対角平均"""
    L, K = X_recon.shape
    N = L + K - 1
    ts_recon = np.zeros(N)
    counts = np.zeros(N)
    for i in range(L):
        for j in range(K):
            ts_recon[i + j] += X_recon[i, j]
            counts[i + j] += 1
    return ts_recon / counts

def ssa_decompose(series, L):
    """SSA: 埋め込み + SVD"""
    X = embed_series(series, L)
    U, s, Vt = svd(X, full_matrices=False)
    return U, s, Vt, X

def select_rank(s, energy_threshold=0.9):
    """自動次数選択（累積寄与率が閾値以上の成分数）"""
    cumulative_energy = np.cumsum(s ** 2) / np.sum(s ** 2)
    r = np.searchsorted(cumulative_energy, energy_threshold) + 1
    return r

def ssa_reconstruct(U, s, Vt, r):
    """指定された次数で再構成"""
    return sum(s[i] * np.outer(U[:, i], Vt[i, :]) for i in range(r))

def ssa_interpolate_s(series, L=40, energy_threshold=0.9):
    """主処理関数：SSA補間＋フィルタリング"""
    series = pd.Series(series)
    # is_nan = series.isna()
    ts = series.fillna(0).to_numpy()

    # SSA分解
    U, s, Vt, X = ssa_decompose(ts, L)

    # 自動ランク選択
    r = select_rank(s, energy_threshold)

    # 再構成
    X_recon = ssa_reconstruct(U, s, Vt, r)
    ts_recon = diagonal_averaging(X_recon)
    return ts_recon

def ssa_interpolate(series, L=40, energy_threshold=0.9):
    """主処理関数：SSA補間＋フィルタリング"""
    if series.ndim == 1:
        return ssa_interpolate_2(series, L, energy_threshold)
    else:
        # 2次元配列の場合は各次元に対してSSA補間を適用
        result = np.zeros_like(series)
        for i in range(series.shape[1]):
            result[:,i] = ssa_interpolate_s(series[:,i], L, energy_threshold)
        return result

以下は平滑化スプラインとの比較例です．

smoothed2 = ssa_interpolate(data[:,0], L=5, energy_threshold=0.9329)
fig, axes = plt.subplots(1, 1, figsize=(10, 6), sharex=True)
plt.plot(t, data[:, 0], 'k.', label='original')
plt.plot(t, smoothed2, 'g--', label='smoothed by SSA')
plt.plot(t, smoothed[:, 0], 'b-', label='smoothed by smoothing spline')
# plt.plot(t, diff1[:, 0], 'r--', label='1st derivative ch0')
plt.legend()

図２：図１の上の例との比較（青：平滑化スプライン．緑破線：SSA．上と下で平滑化スプライン
　　の重み係数が異なる）

この比較はノイズの多い例なので，良し悪しを判断するには微妙です．

しかし，平滑化スプラインは，データを関数で記述するので，微分も可能で，かつ関数の微分を計算するので，もとのデータに対しても整合性のある微分値を返してくれます．

力学は二階微分まで行うので，少しでも性能の良い平滑化スプラインのほうが優れています．

ただし，速度（一階微分）しか計算しない場合はSSAも悪くはないでしょう．

おわりに

今回の記事の目的は，平滑化スプラインのコードを見直しです．前のコードが汚かったので訂正をしたかっただけですが，ついでにSSAとの比較も行いましたが，改めて平滑化スプラインの性能の良さを認識しました．

平滑化スプラインには，苦手な計算もあります．たとえば，データ長が大きいと，計算時間が非線形的に長くなります．それらの問題を克服するためには，その背後の数理を理解することで，多くの問題は解決できます．ただし，生データの精度が悪すぎる場合は，どのようなフィルタでも対処しようがありません．たとえばDLTで計測したデータを用いて力学計算を行うことは，ほとんどの場合，平滑化スプラインを使用しても無謀な試みとなるでしょう．

また，ノイズ分布が一様ではないとき，そのノイズに合わせる必要があります．また重み係数の設定が，試行錯誤的になります．平滑化スプラインは万能ではありません．

基本はできるだけ精度の高い計測を行うことです．

なお，空間と時間の分解能を上げ，データの精度を向上することで，平滑化スプラインを使用しなくてもある程度計算できますが，二階微分の演算を伴う力学計算では可能な範囲で平滑化スプラインを使用すべきでしょう．微分演算を含むわけですから．

なお，MatLabの平滑化スプラインは，３次のスプライン関数までしか対応していないので，加速度計算はできません．４次関数は二階微分によって２次関数になりますが，３次関数を二階微分すると１次関数になりますので，加速度が滑らかな波形になりませんので．ご注意ください．

繰り返しになりますが，微分演算を行うなら平滑化スプラインを使うことを強くおすすめします．もともとのデータの精度がよほど良くない限り，二階微分演算によって算出した力の信頼性が問われます．

二階微分演算を通したトルクなどを計算しても，トルクの波形を示さず，統計データだけを示し，それがどの程度の信頼性があるかよくわからない研究も多いのではないでしょうか．微分演算はフィルタリングの性能に著しく影響を受けます．バターワースフィルタで計算された力やトルクは，条件が整わない限り，どの程度の信頼性があるのかと疑いの目で見たほうがよいでしょう．