本日の給油(CRF1100Lアフリカツイン/2BL-SD10)

給油日 オドメーター (km) 給油量 (L) 単価 (円/L) 燃費 (km/L) 距離単価 (円/km)
2022-08-27 3401 16.63 157 16.42 9.56
2022-11-03 3727 17.20 158 18.95 8.34
2023-02-04 3987 16.75 157 15.52 10.12
2023-04-20 4165 10.37 157 17.16 9.15
2023-04-30 4563 18.22 160 21.84 7.32
2023-06-17 4982 16.41 159 25.53 6.23
2023-09-02 5296 19.46 178 16.14 11.03
2023-10-01 5700 18.89 180 21.39 8.42
2024-01-28 6003 18.62 165 16.27 10.14
2024-03-22 6292 14.10 165 20.50 8.05

千葉でのキャンプに向けて給油。中距離くらいのお出かけに使っていたので燃費はやや回復したかな。

週末キャンプでScaniverseのGaussian Splattingを試してみました。右側のスクリーンショットは写真じゃなくて3Dモデルなのよ。とてもリアルでびっくり。

ランク付けにおける重複なし確認のやりかた

複数個の対象についてランク付けをしてもらいたいことがあります。例えば、6種類のソーダ水をおいしさ順に並べてもらい、ソーダ水どうしのおいしさの差を見るなどが考えられます。今回は、そのデータをとるときに「順位の重複を許さないように判定する方法」を考えたいと思います。6種類を「5位・2位・2位・3位・6位・1位」などと評価してしまうと、2位が重複するとともに4位がなくなってしまいます。順序はさておき1位から6位までが一回ずつ含まれていると確認するにはどうすればよいでしょうか。

いちばん簡単な方法は集合の関数を使うことです。setdiffは二つの集合の差をとってくれるので、numbersに入っていてhogeに入っていない数が返ってきます。

numbers <- c(1, 2, 3, 4, 5, 6)
hoge <- c(5, 2, 2, 3, 6, 1)
setdiff(numbers, hoge)   #=> 4

集合が使いにくいプログラミング言語で同様のことをしたいときに、「総和と総積が一致すればよいのでは」と考えてみました(この後、この方法ではうまくいかないことが分かります)。前述の例だと、総和は1+2+3+4+5+6=21、総積は1×2×3×4×5×6=720です。一方でsum(hoge)は19、prod(hoge)は360と、さっきの値と一致しないので評価値の重複や不足が推測できます。

この方法を実用するためには、本当にすべての場合でそうなるかを確かめなければなりません。すべての場合を調べるのは大変なので、以下のように「総和と総積の一致」が差集合と同じようにチェック法として使えるのかをモンテカルロ的に調べてみました。Nは対象の個数とします。総和と総積が想定したものと一致しているのに差集合が空集合でないときには、その数列を表示するというプログラムです。

N <- 8
numbers <- seq(from=1, to=N)
p <- prod(numbers)
s <- sum(numbers)

for (k in 1:1000000) {
  hoge <- sample(numbers, replace=TRUE)
  if (prod(hoge)==p && sum(hoge)==s) {
    if (length(setdiff(numbers, hoge)) != 0) {
      print(hoge)
    }
  }
}

数列をランダムに生成して100万回の確認をしますが、何度実行してもNが8以下の場合は何も表示されないので、おそらく「総和と総積が一致すれば、1~Nは一度ずつ登場している」といえそうです。

ところが、Nが9以上になると話が変わってきます。たとえばN=10のときの「2 7 10 5 4 9 4 4 1 9」という数列は

> sum(c(2, 7, 10, 5, 4, 9, 4, 4, 1, 9))
[1] 55
> prod(c(2, 7, 10, 5, 4, 9, 4, 4, 1, 9))
[1] 3628800

のように1~10が一回ずつ出てくる時と同じ値になります。N=20までは確認し、いずれも何かしらの数列が表示されましたので、少なくともNが9~20のときには「総和と総積の一致」は重複・不足の確認には使用できないことがわかりました。21以上のNについても同様だと予想したのですが、これ、どこかで証明されてるんですかね。

本日の給油(スーパーカブ110プロ/JA07)

給油日 オドメーター (km) 給油量 (L) 単価 (円/L) 燃費 (km/L) 距離単価 (円/km)
2022-08-31 13546.7 2.92 156.16 48.53 3.22
2022-11-05 13699.4 3.21 157.01 47.57 3.30
2023-02-06 13848.1 3.26 157.06 45.61 3.44
2023-03-02 13973.9 2.60 156.92 48.38 3.24
2023-04-01 14095.9 2.99 156.86 40.80 3.84
2023-07-13 14218.0 2.73 163.00 44.73 3.64
2023-09-25 14380.5 3.60 176.94 45.14 3.92
2023-11-25 14529.5 2.93 162.12 50.85 3.19
2024-02-07 14644.2 2.80 165.00 40.96 4.03
2024-03-04 14775.6 2.93 164.85 44.85 3.68

二拠点間をよく乗っていたので、1ヶ月で給油となりました。まだ燃費は悪いままですが、今月は3速を多様したライディングが多かったからかもしれません。

信号の包絡線を計算する (Common Lisp)

信号の包絡線 (signal envelope) は「信号の頂点にシーツをかぶせたときのような概形」を持つ波形です。

MatlabRで書いた信号包絡線の計算プログラムCommon Lispに移植しました。信号包絡線の計算方法は振幅の絶対値に対して、(1)低域フィルタをかける、(1')Savitzky-Golayフィルターがよい、(2)ピーク検出をしてスプラインでつなぐ、(2')いやいや秋間補間がよい、(3)ヒルベルト変換を使う、(3')結局それに低域フィルタかけるじゃろ、(4)スペクトル包絡には線形予測符号も使われてるよね、などなどいくつも提案できるのですが、ここではヒルベルト変換を用いて解析信号を計算する方法でいきます。解析信号の絶対値を計算すると信号包絡線になります。

まずは解析信号の実装。フーリエ変換をして、負の周波数成分をゼロにして、逆フーリエ変換します。

(defun hilbert (x)
  "Compute analytic signal of real signal vector X using Hilbert transform."
  (let* ((xx (napa-fft:fft (zeropad-pow2 x)))
         (yy (make-array (length xx) :initial-element #C(0 0)))
         (zz nil)
         (out (make-array (length x))))
    (setf (aref yy 0) (aref xx 0))
    (loop for n from 1 below (/ (length xx) 2)
          do (setf (aref yy n) (* 2 (aref xx n))))
    (setf (aref yy (/ (length xx) 2)) (aref xx (/ (length xx) 2)))
    (loop for n from (1+ (/ (length xx) 2)) below (length xx)
          do (setf (aref yy n) 0))
    (setq zz (napa-fft:ifft yy))
    (loop for n from 0 below (length out)
          do (setf (aref out n) (aref zz n)))
    out))

Matlabhilbert()と比較してみましょう。

>> hilbert([1 2 -3 4])
ans =
   1.0000 + 1.0000i   2.0000 + 2.0000i  -3.0000 - 1.0000i   4.0000 - 2.0000i
CL-USER> (hilbert #(1 2 -3 4))
#(#C(1.0d0 1.0d0) #C(2.0d0 2.0d0) #C(-3.0d0 -1.0d0) #C(4.0d0 -2.0d0))

同じ結果になっていますね。ヨシ。

次に信号包絡線を計算します。

(defun envelop (x)
  "Calculate amplitude envelop of a signal vector X."
  (let* ((xh (hilbert x))
         (xe (make-array (length xh))))
    (dotimes (n (length xh))
      (setf (aref xe n) (abs (aref xh n))))
    xe))

こちらもMatlabと比較してみましょう。

>> abs(hilbert([1 2 -3 4]))
ans =
    1.4142    2.8284    3.1623    4.4721
CL-USER> (envelope #(1 2 -3 4))
#(1.4142135623730951d0 2.8284271247461903d0 3.1622776601683795d0
  4.47213595499958d0)

本物の音響信号でやってみましょう。スネアドラムを一発叩いた音の録音があるので、信号包絡線を計算してみます。wavreadは以前作ったプログラムです。(→Common Lisp用のwavread - 丸井綜研

(ql:quickload :clgplot)
(let* ((snd (wavread "snare04-short.wav"))
       (x (first snd))
       (fs (second snd))
       (xe (envelop x))
       (xel (filt (biquad-lpf 100 (sqrt 1/2)) xe))
       (tt (list->vector (mapcar #'(lambda (v) (/ v fs)) (range (length xe))))))
  (clgp:plot x :x-seq tt)
  (clgp:plot xel :x-seq tt))

結果は以下のようになりました。左側が元の信号で右側が包絡線です。横軸は時間(秒)で縦軸は振幅(プラスマイナス1の範囲に正規化)。少しスムーズにするために包絡線にカットオフ周波数100 Hzの低域フィルタをかけてあげています。

Aseprite買った

以前から欲しいと思っていたAsepriteというドット絵作成ツールが半額になっていたので、思わず購入してしまった。半額セールはたぶん2月13日まで。

www.aseprite.org

で、初めて描いてみたのが以下の絵。16×16ピクセルBrüel & Kjær Type 4128-Cさんを想像で描きました。

本物の写真を見てみたら、チョーカーみたいなのつけてるし、目はなかった。ドット絵のほうは緑色のゾンビか、サボテンみたい。

昔ね、素焼きの陶器でできた人形の頭かなんかにカイワレの種が埋まってて水をかけるとニョキニョキ育つというおもちゃがあったのよ。カイワレの芽が髪の毛みたいに見えて、ヘアカットするとカイワレが収穫できるというもの。それを思い出してしまいましたよ。(後日、調べてみたらChia PetとかChia Headというものだったらしい)

www.youtube.com

本日の給油(スーパーカブ110プロ/JA07)

給油日 オドメーター (km) 給油量 (L) 単価 (円/L) 燃費 (km/L) 距離単価 (円/km)
2022-07-27 13405.0 2.21 157.92 57.29 2.76
2022-08-31 13546.7 2.92 156.16 48.53 3.22
2022-11-05 13699.4 3.21 157.01 47.57 3.30
2023-02-06 13848.1 3.26 157.06 45.61 3.44
2023-03-02 13973.9 2.60 156.92 48.38 3.24
2023-04-01 14095.9 2.99 156.86 40.80 3.84
2023-07-13 14218.0 2.73 163.00 44.73 3.64
2023-09-25 14380.5 3.60 176.94 45.14 3.92
2023-11-25 14529.5 2.93 162.12 50.85 3.19
2024-02-07 14644.2 2.80 165.00 40.96 4.03

ストップ・アンド・ゴーばかりの短距離運用だったので、燃費は過去最悪に近いです。寒い時期は燃費が悪くなりがちなのはしょうがないですね。

本日の給油(CRF1100Lアフリカツイン/2BL-SD10)

給油日 オドメーター (km) 給油量 (L) 単価 (円/L) 燃費 (km/L) 距離単価 (円/km)
2022-08-27 3401 16.63 157 16.42 9.56
2022-11-03 3727 17.20 158 18.95 8.34
2023-02-04 3987 16.75 157 15.52 10.12
2023-04-20 4165 10.37 157 17.16 9.15
2023-04-30 4563 18.22 160 21.84 7.32
2023-06-17 4982 16.41 159 25.53 6.23
2023-09-02 5296 19.46 178 16.14 11.03
2023-10-01 5700 18.89 180 21.39 8.42
2024-01-28 6003 18.62 165 16.27 10.14

前回アフリカツインに火を入れたのが11月23日のデイキャンプだったので、じつに2か月以上も放置してしまいました。前に乗っていた某大型バイクだと1か月エンジンをかけないとバッテリーが上がるという謎があったのでヒヤヒヤしましたが、さすがは安定のホンダ車。一発始動です。高価なリチウムイオン電池を積んでいると聞いたのと、交換が難しい箇所にあるとのことなので、工賃もかかるバッテリー交換はしたくない……(パーツカタログで調べたらエリーパワーのHY110https://amzn.to/49r3e6H)だそうです。5万円近くしますね……)

前回の単価は180円でしたが、会員価格で今回は165円。

今日は東京都~埼玉県を走り回り、ようやく6000 km超えました。