マイコンと加速度センサでタコメータを作る - 14. メータモード実装

回転数計算処理

回転数測定のスペックをちゃんと考えてみると、思ったより刻みが粗かった(187.5[rpm]で64ステップ程度)ので、ちょっと検討してみます。

実データ確認

FFTのグラフで、ピーク前後のスペクトルを使ったりでなんとか周波数をより細かく出せないか…

Octaveでちょっと加速度の時間波形を出してみます。

>> plot(gData(800*512+1:801*512,2))

f:id:nokixa:20201210002853p:plain:w600

この波形を見ると、21個の山があるように見えます。
正弦波を重ね書きしてみます。

>> hold on
>> plot(sin(2*pi*21*[1:512]/512))

f:id:nokixa:20201210004613p:plain:w600

hold onで前のグラフを消さずに上書きできます。
この正弦波とだいたい形は合っていますが、ちょっとずつずれていっています。

FFT結果を再確認します。ピーク周波数周辺を見てみます。

>> figure(2)
>> plot(abs(gFFT(1:100,800)), "*-")
>> xticks([1:1:800])
>> axis([11 30 100 300])

f:id:nokixa:20201210005850p:plain:w400

このグラフでは、横軸の値1がDC成分を表すので、横軸の値22が先ほどの正弦波周波数にあたります。確かにここがピークになっています。
1つ前後の周波数もそれなりの値を持っているので、使えそうな気も。

周波数計算検討

FFT対象の波形が本当に純粋な単一周波数、ただしサンプリング周波数の整数倍でないとします。

$f(t) = A\exp(j(\omega_n + \omega_f)t)$

$\omega_n = n\Delta\omega$ , $\Delta\omega=\cfrac{2\pi}{N\Delta t}$ , $\Delta t$ : サンプリング周期, $n$ : 整数, $N$ : サンプリング点数

$| \omega_f | \lt \Delta\omega / 2$

ここで、角周波数 $\omega_m = m\Delta\omega$ の成分を取り出すべく、 $\exp(-j\omega_m t)$ を掛けてサンプリング時間積分、サンプリング時間で割ると、

$\cfrac{1}{N\Delta t} \int_0^{N\Delta t}\exp(j(\omega_n+\omega_f -\omega_m)t)dt$

$=\cfrac{1}{N\Delta t (\omega_n+\omega_f -\omega_m)}[\exp(j(\omega_n+\omega_f -\omega_m)t) ] _0 ^{N\Delta t}$

$=\cfrac{1}{N\Delta t (\omega_n -\omega_m +\omega_f)}(\exp(jN\omega_f \Delta t)-1)$

この絶対値 $r$ は、

$r = \cfrac{1}{N\Delta t|\omega_n -\omega_m +\omega_f|} \sqrt{(\cos(N\omega_f \Delta t)-1) ^2 + \sin ^2(N\omega_f \Delta t) }$

$=\cfrac{1}{N\Delta t|\omega_n -\omega_m +\omega_f|} \sqrt{2-2\cos(N\omega_f \Delta t)}$

$1-\cos\theta = 2\sin ^2\cfrac{\theta}{2}$ より、

$r = \cfrac{2|\sin(N\omega_f \Delta t)|}{N\Delta t|\omega_n -\omega_m +\omega_f|}$

$n$ , $m$ の組み合わせについて、以下の3パターンを考えます。

$m = n$ ・・・ $m$ は、 $\omega_m$ が実際の角周波数に最も近くなる整数値
$m = n-1$ ・・・ $m$ は、a. の1つ前
$m = n+1$ ・・・ $m$ は、a. の1つ後

それぞれについて、FFTスペクトル絶対値 $r$ は、

$r_a = \cfrac{2|\sin(N\omega_f \Delta t)|}{|N\omega_f\Delta t|}$
$r_b = \cfrac{2|\sin(N\omega_f \Delta t)|}{|N\omega_f\Delta t + 2\pi|}$
$r_c = \cfrac{2|\sin(N\omega_f \Delta t)|}{|N\omega_f\Delta t - 2\pi|}$

これについて考えてみます。

まず、a. が最も大きい値となります。
b. とc. については、 $\omega_f$ の符号により大小関係が変わります。

$\omega_f \lt 0$ の場合・・・b. のほうが大きい
$\omega_f \gt 0$ の場合・・・c. のほうが大きい

これらをグラフにしてみると、こうなります。横軸に $N\omega_f \Delta t$ を、縦軸に上の3つのスペクトルそれぞれの $r$ をプロットしています。

>> omega_f = [-pi:pi/1000:pi]
>> r_a = 2*sin(omega_f)./omega_f
>> r_b = 2*sin(omega_f)./(omega_f + 2*pi)
>> r_c = 2*sin(omega_f)./(omega_f - 2*pi)
>> plot(omega_f, r_a, omega_f, r_b, omega_f, r_c)

f:id:nokixa:20201211005733p:plain:w400

まず、ピーク周波数前後のスペクトルの大小関係により、ピーク周波数よりどちらに本来の周波数があるかわかります。
また、ピーク周波数前後のスペクトルの大きさの比から、 $\omega_f$ の値が推定できるのでは、と考えられます。この比をプロットすると、

>> plot(omega_f, r_c ./r_b)

f:id:nokixa:20201212004244p:plain:w400

と、 $\omega_f$ の値にしたがって単調増加する形になります。
なので、2つの比で $\omega_f$ が一意に決まります。

※ここで、./は配列の要素ごとの除算になります。

$r_a$ と $r_b$ または $r_c$ の比でもいいかもしれません。
下のようになります。

>> plot(omega_f, r_b ./ r_a)
>> hold on
>> plot(omega_f, r_c ./ r_a)
>>

f:id:nokixa:20201211235802p:plain:w400

こっちだと、単純な単調変化とならないので、最初に $r_b$ , $r_c$ の大きさ比較を行う必要が出てきます。
前者のもののほうがシンプルに思うので、こちらでいきたいと思います。
ただし、 $r_b$ , $r_c$ の値があまり小さかったら中心の周波数をそのまま採用することとします。大きさの基準としては、スペクトル絶対値の平均の2倍、というくらいかなと。

実データで検証

元のデータに戻って、前後のスペクトルの大きさの比を調べてみます。

>> r_21 = abs(gFFT(21,800))
r_21 =  141.03
>> r_22 = abs(gFFT(22,800))
r_22 =  278.40
>> r_23 = abs(gFFT(23,800))
r_23 =  126.31

1つ前のスペクトルのほうが大きいので、ピーク周波数より少し小さい周波数が本来の周波数のようです。
前後の比およびそこから $\omega_f$ を逆算すると、

>> r_23 / r_21
ans =  0.89568
>> [w, iw] = min(abs(r_c ./ r_b - r_23/r_21))
w =  0.000056458
iw =  891
>> omega_f(iw)
ans = -0.34558
>> omega_f(iw)/pi
ans = -0.11000

となります。だいたい周波数単位の0.11倍小さいところが本来の周波数らしい。
加速度の時間波形にこれを重ねてみます。

>> plot(gData(800*512+1:801*512,2))
>> hold on
>> plot(sin(2*pi*(21-0.11)*[1:512]/512))

f:id:nokixa:20201211014848p:plain:w600

結構いい線行ってる！周期がばっちり。
これでやってみたいと思います。

時系列FFTグラフから回転数列を出してみたいと思います。
以下のOctave スクリプトを作りました。

function rpm = calc_rpm(fft_abs, df, amp_ratio)

  nt = size(fft_abs)(2);
  nf = size(fft_abs)(1);
  fft_mean = mean(fft_abs(2:end));
  
  for i = 1:nt
    [v, iv] = max(fft_abs(2:end, i));
    v_m1 = 0; v_p1 = 0;
    if(iv >= 2)
      v_m1 = fft_abs(iv-1, i);
    endif
    if(iv != nf)
      v_p1 = fft_abs(iv+1, i);
    endif
    
    if( (v_m1 > fft_mean) && (v_p1 > fft_mean) )
      [r, ir] = min(abs(amp_ratio - v_p1/v_m1));
      rpm(i) = iv + (ir / size(amp_ratio)(2) - 0.5);
    else
      rpm(i) = iv;
    endif
    
  endfor
  
endfunction

これを動かしてみます。

まず、今回の処理を追加する前の時系列FFTグラフ。上記のスクリプトで、 rpm(i) = iv + (ir / size(amp_ratio)(2) - 0.5); の部分をrpm(i) = iv; で置き換えて表示しました。

>> rpm = calc_rpm(gFFT_abs(1:64,:), 1, r_c./r_b);
>> plot(rpm)

f:id:nokixa:20201212022046p:plain:w600

処理を追加すると、

f:id:nokixa:20201212021725p:plain:w600

うーん…それほどギザギザ感がなくなったかというと…
多少はよくなったぐらいか。
とりあえずこれでやってみようと思います。

マイコン用テーブル作成

マイコンでの実装としては、前後周波数成分比のテーブルを用意して、実際の測定値での前後周波数成分比が入る範囲を探す、ということをしようと思います。
テーブルの要素数分だけ分解能が上がることになります。
今回の設定では187.5[rpm]の分解能なので、20要素用意すれば10[rpm]くらいの分解能を持つことができます。

先ほどの計算から、テーブルの値として使うものを取り出しておきます。

>> r_c(1:100:end)./r_b(1:100:end)
ans =

 Columns 1 through 17:

   0.33333   0.37931   0.42857   0.48148   0.53846   0.60000   0.66667   0.73913   0.81818   0.90476       NaN   1.10526   1.22222   1.35294   1.50000   1.66667   1.85714

 Columns 18 through 21:

   2.07692   2.33333   2.63636   3.00000

>>

$\omega_f = 0$ のところでは $r_b$ , $r_c$ いずれも0になるので、NaNが出てしまっています。
ここは1.0で埋めておけばいいかと思います。

ここまで

一旦ここまでにします。

今回分かったことはどういうことかというと、

FFTを使って周波数解析を行う場合、周波数分解能は $\cfrac{1}{N\Delta t}$ ( $\Delta t$ : サンプリング間隔、 $N$ : サンプリング数)となる
周波数分解能を上げたい場合、 $N\Delta t$ を大きくする必要があるが、そうするとサンプリング時間が長くなり、結果処理時の更新レートが下がってしまう
ピーク周波数の前後の周波数成分を使うことで、の非整数倍の周波数成分検出ができ、短いサンプリング時間のままで周波数分解能を上げることができる
- ただし対象の波形には単一周波数しか含まれていない前提