実車とシミュレーションの比較にはFFT(高速フーリエ変換)を使おう
はじめに
実車データとシミュレーション(SIM)データを比較する際、「ぱっと見で合っていそう」という判断でも良いですが、数値で定量的に表すときはFFT(高速フーリエ変換)を使用すると良いでしょう。MATLABではfft関数がこれに相当します。
FFTが有効なデータ
ある走行シーンで実車データとSIMデータがあるとき、以下の条件を満たすデータであればFFTによる比較が有効です。
- 実車データでなくても、すでにデータ取得を完了しているもの
- 0を中心に変化しているもの
- 特定の値(〇〇)を中心に変化しているもの
これらの条件を満たせば、どんなデータでも比較する際に使用できます。
FFT活用の3つの重要ポイント
ポイント1: 最初の点と最後の点を一致させる
FFTを掛ける範囲を設定する際、最初の点と最後の点がほぼ一致していなくてはなりません。
実装のコツ
0中心のデータであれば、0を通過するタイミングの2点でFFTをかけるのが簡単です。
% 0クロス点を見つけてFFT範囲を設定
idx_start = find(data >= 0, 1, 'first');
idx_end = find(data(idx_start+1:end) >= 0, 1, 'first') + idx_start;
fft_data = data(idx_start:idx_end);
なぜ重要なのか
開始点が0で終了点が100のとき、FFTは周期的なデータを仮定するため、100から0まで急激な変化が発生してしまうからです。
必ず最初と最後の点を繋げるようにしてください。
ポイント2: ナイキスト周波数を理解する
データのサンプリング周波数の半分までの周波数しか見ることはできません。
具体例
- サンプリング周波数が100Hzならば、50Hzまでの算出結果しか信用できない
理由
これは直感的に理解できると思います。100Hzのデータで100Hzの波の山と谷を捉えることができないからです。
たとえ谷や山の1点が取れたとしても、100Hzではない別の周波数だと判定されてしまうでしょう。
実務上の注意点
100Hzのデータがあるとき、50Hz以下の成分しかないのであれば正確に算出できるでしょう。
しかし、現実のデータの場合、50Hz以上のノイズも含まれているため、単純に「100Hzだから50Hz以下は絶対に正しい」というわけではありません。
エイリアシング対策
サンプリングした後で50Hz LPF(ローパスフィルタ)などをして高周波を削除しても、すでに遅いのです。
% 悪い例: サンプリング後にフィルタリング
sampled_data = sample(raw_data, 100); % 100Hzでサンプリング
filtered_data = lowpass(sampled_data, 50); % 後からLPF → 手遅れ
% 良い例: サンプリング前にアンチエイリアシングフィルタ
filtered_data = lowpass(raw_data, 50); % 先にLPF
sampled_data = sample(filtered_data, 100); % その後サンプリング
サンプリングする際に50Hz以下しか取らないようにする必要があります。(アンチエイリアシングフィルタの適用)
ポイント3: 窓サイズは2の累乗にする
窓サイズは2の累乗を強く推奨します。
理由
圧倒的に計算速度が変化するため、2の累乗の窓サイズにするようにしてください。
% 推奨される窓サイズ
window_sizes = [256, 512, 1024, 2048, 4096, 8192];
% 例: 2048サンプルでFFT
N = 2048;
fft_result = fft(data(1:N));
窓サイズと結果の関係
常に一定の周波数のデータであれば、どんな窓サイズにしても同じ結果になります。
もし窓サイズを変えて結果が違うということは、1つ目のポイントである「最初と最後の点の数値が違う」ことが予想されます。
まとめ
実車データとSIMデータをFFTで比較する際の重要ポイント:
- 最初の点と最後の点を一致させる → 0クロス点などを利用
- ナイキスト周波数を理解する → サンプリング周波数の半分まで信頼可能
- 窓サイズは2の累乗にする → 計算速度の最適化
これらのポイントを押さえることで、定量的で信頼性の高いデータ比較が可能になります。
参考: MATLAB実装例
% サンプリング周波数
fs = 100; % Hz
% 0クロス点を見つける
zero_crossings = find(diff(sign(data)));
start_idx = zero_crossings(1);
end_idx = zero_crossings(2);
% 2の累乗に調整
N = 2^nextpow2(end_idx - start_idx);
fft_data = data(start_idx:start_idx+N-1);
% FFT実行
Y = fft(fft_data);
P2 = abs(Y/N);
P1 = P2(1:N/2+1);
P1(2:end-1) = 2*P1(2:end-1);
% 周波数軸
f = fs*(0:(N/2))/N;
% プロット
plot(f, P1);
xlabel('周波数 (Hz)');
ylabel('振幅');
title('FFTスペクトル');