凛世の自転周期を計測する

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">■緊急告知■<br>杜野凛世フィギュア 100時間連続配信をYoutubeで実施します！！<br><br>GW中、回転台の上に乗せた凛世のフィギュアを100時間放送し続けるだけの配信です。<br><br>開始時間は明日2日(土)の11時！<a href="https://twitter.com/hashtag/%E3%82%B7%E3%83%A3%E3%83%8B%E3%83%9E%E3%82%B9?src=hash&amp;ref_src=twsrc%5Etfw">#シャニマス</a> <a href="https://twitter.com/hashtag/%E5%87%9B%E4%B8%96100%E6%99%82%E9%96%93%E9%85%8D%E4%BF%A1?src=hash&amp;ref_src=twsrc%5Etfw">#凛世100時間配信</a><br><br>⬇️待機所⬇️<a href="https://www.youtube.com/watch?v=kZZt3OCrbDU&feature=youtu.be">https://www.youtube.com/watch?v=kZZt3OCrbDU&feature=youtu.be</a> <a href="https://t.co/khhkrddmTc">pic.twitter.com/khhkrddmTc</a></p>&mdash; amiami_product (@amiami_products) <a href="https://twitter.com/amiami_products/status/1256146340423389185?ref_src=twsrc%5Etfw">May 1, 2020</a></blockquote>

[【シャニマス】杜野凛世フィギュア 100時間耐久配信！！ 〈amiami〉](https://www.youtube.com/watch?v=kZZt3OCrbDU)

体験版: [【シャニマス】杜野凛世フィギュア 30分間ミニ耐久配信！！ 〈amiami〉](https://www.youtube.com/watch?v=Lpb5xxKDALM)

---

> いえ、プロデューサーさまに、お喜びいただく……
> それだけが、凛世の幸せなのです

いや、なんていうか……
ちょっと聞いてみたいんだけど

> はい。どのような、ことでも

なんで 100 時間回転してるんだ?

# 『なんのために……』

> まわる……まわる……

ふと思いたって差分をとってみたら時間帯によって自転周期が変化していることに気付いてしまったのだ……。

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr"><a href="https://twitter.com/hashtag/%E5%87%9B%E4%B8%96100%E6%99%82%E9%96%93%E9%85%8D%E4%BF%A1?src=hash&amp;ref_src=twsrc%5Etfw">#凛世100時間配信</a><br>だんだん速くなってた <a href="https://t.co/rD9Subk0DD">pic.twitter.com/rD9Subk0DD</a></p>&mdash; イショティハドゥス (@Ishotihadus) <a href="https://twitter.com/Ishotihadus/status/1256889957462929408?ref_src=twsrc%5Etfw">May 3, 2020</a></blockquote>

いやまさかループじゃなくて本当にライブ放送だったとは……。

しかもリアルタイムである。

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">地震で揺れる凛世 <a href="https://twitter.com/hashtag/%E5%87%9B%E4%B8%96100%E6%99%82%E9%96%93%E9%85%8D%E4%BF%A1?src=hash&amp;ref_src=twsrc%5Etfw">#凛世100時間配信</a> <a href="https://t.co/38xkl8NG54">pic.twitter.com/38xkl8NG54</a></p>&mdash; イショティハドゥス (@Ishotihadus) <a href="https://twitter.com/Ishotihadus/status/1257717314683760645?ref_src=twsrc%5Etfw">May 5, 2020</a></blockquote>

[地震情報](https://www.jma.go.jp/jp/quake/20200505170115393-06015718.html)

> どうぞ……お付き合いくださいませ……
> プロデューサーさま……

# streamlink / youtube-dl

> ……
> 『恋（ライブ hls）は野の鳥』……

> 『来たと思えば（expire して）飛んでいく』……

## [youtube-dl](https://github.com/ytdl-org/youtube-dl)

何も考えなくてもリンク貼るだけで動いて便利。

ライブ HLS はどうしても ffmpeg に通されるらしい。

ライブ時は `--hls-use-mpegts` というおまじないがある。異常終了しても再生可能なデータが残る。

## [streamlink](https://github.com/streamlink/streamlink)

出力ファイル名はもちろん画質すら与えないといけない。何か考えないと動かない。

ts そのままダウンロードできる。便利。ts にしか乗らない情報がほしい場合にはこっちを使うとよい。

## youtube-dl & streamlink V.S. 6 時間の壁

> 3 つめ（のログ）は……
> 何色だったでしょう……

> 紅だったかも……
> しれません……

YouTube のダウンローダは、[get_video_info](https://kitoko552.hatenablog.com/entry/20180311/1520743836) にアクセスして、player_response にある json の streamingData/hlsManifestUrl に書いてある m3u8 のマニフェストをダウンロードしてきて、さらにその中にある m3u8 プレイリストを読みまくることで実現されている。

```sh
curl 'https://www.youtube.com/get_video_info?video_id=kZZt3OCrbDU' | 
ruby -r json -r uri -e "puts JSON.parse(URI.decode_www_form(gets).to_h['player_response'])['streamingData']['hlsManifestUrl']"
# => https://manifest.googlevideo.com/api/manifest/hls_variant/expire/ ... /file/index.m3u8
```

このマニフェストやプレイリストはちょうど 6 時間で失効してしまうので、同じマニフェストは 6 時間しか連続して使えないことになる。
もちろん再発行すれば済む話だが、どうも youtube-dl も streamlink もこれに対応しておらず、残念ながらこれらを使うかぎり連続 6 時間までしか保存できない。
streamlink は一時期対応していたらしいが、どうやらこの挙動は最近の挙動らしく、コケる。すかぽんたん。

どちらも負けでございます。

あたらしいツールを作ったり streamlink をいじったりする気力はさすがにないので、（連続している必要がないなら）シェルで while 文を回したほうがはやい。今回はめんどうくさいのでそうした。

ゆうて ts つなげるだけだしたぶん 3 時間くらいで書けそう。誰かやって。

まあバッファの長さが 30 秒なので、短時間で切り替えられるなら後処理でも実はなんとかなる。

ちなみに Web 版の YouTube 本体はどうしているのかというと、どうやら 6 時間ごとにリロードがはいる。ちらっとリロード表示が出る。凛世ずっと回してたからわかる……。

# 『どうやって……どこから……』

## 手法

凛世はずっとずっと同じ場所で同じように回り続けるので、あるフレームの凛世を適当にキャプチャして（リファレンス凛世）、毎フレームの凛世とリファレンス凛世の差分をとって、凛世の差異を計算すればよい。

![凛世はいつも回っている。地球が北極星の方角から見ると反時計周りに回っているように、凛世も上から見ると反時計周りに回っている。](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78668/d261f3dc-00ee-bce0-6c3f-f7716f2974dd.png)

各フレームの PSNR を計算しようと思い MATLAB を取りだし Computer Vision Toolbox をせこせことインストールしたが、MATLAB にはそもそも VideoReader があるし、VideoReader はクソ遅いし、ありがとう MATLAB、いいやつだったよ。手法の挙動確認には便利なんだ。この恩は忘れない。

> 『身の回りの MATLAB Toolbox が、少し増えました』
> 『大変愉快なものばかりです』

[ffmpeg でできる](https://ffmpeg.org/ffmpeg-filters.html#psnr)。後処理がめんどうくさいので sed をかました。凛世は Darwin の上で回転してしまったので、BSD sed が使われている。ところでこの「かます」は俗語なので中高型ではなく平板型である。

```sh
ffmpeg -i りんぜ.ts -i reference.png -filter_complex psnr=- -an -f null /dev/null | sed -Ee 's/n:([0-9]+) .* psnr_y:([0-9.]+).*/\1 \2/' > りんぜ.txt
```

## セグメントロスの罠

YouTube は 5 秒のセグメントごとにわけて送ってくるが、このセグメントが落ちることがあるらしく（受信側の問題なのか配信者の問題なのか YouTube 側の問題なのかよくわからない）、5n 秒凛世が瞬間移動してしまう場合がある。これは外れ値検出を上手にやれば回避できそうだ。

> セグメント……
> セグメント…………

> 凛世を……
> 慰めて……くれるのですか……

どうでもいいが、YouTube はどうも完全に連続したライブ配信をするのは苦手らしく、何十時間と連続して見ていると挙動が怪しくなったりするときも訪れるらしい。YouTube の動画配信サーバが応答を返してくれなくなったりする。ほとんどの視聴者が追いだされたこともあった。これはもうその都度 Ctrl-C するしかない。

## 実装

ファイルの作成日時を stat コマンドで確認しながら ffmpeg をちまちまと叩いて psnr を計算しまくる（これは途中の画像）。

![回りながら分析される凛世の図](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78668/6361d550-3f74-097b-7d0f-7efb886d63e1.png)

次にこれを MATLAB スクリプトでどん。

<details><summary>誰も使わないであろうスクリプト</summary><div>

```matlab
function [x, y, v] = script
% x: ピークの時刻のリスト
% y: そのときの凛世の自転周期のリスト
% v: そのときの PSNR
files = dir('rinze_*.txt');
paths = arrayfun(@(e) fullfile(e.folder, e.name), files, 'un', 0);
names = arrayfun(@(e) e.name, files, 'un', 0);
tokens = regexp(names, '^rinze_(\d{4})(\d{2})(\d{2})_(\d{2})(\d{2})(\d{2})\.txt$', 'tokens');
ctimes = cellfun(@(t) datetime(cellfun(@str2double, t{1})), tokens, 'un', 0);
[x, y, v] = cellfun(@process_file, paths, ctimes, 'un', 0);
x = cell2matdatetime(x);
y = cell2mat(y);
v = cell2mat(v);


function [x, y, v] = process_file(path, ctime)
FPS = 30 * 1000 / 1001;
LEN = FPS * 30;

file = fopen(path);
list = fscanf(file, '%d %f', [2 Inf]);
fclose(file);

% 初期位置は 30 秒すすめたところ
s = LEN;

peaks = [];
while true
    e = s + LEN * 2;
    if e > size(list, 2)
        break
    end
    % 基準点から 60 秒間で最も PSNR が高い時刻をとってくる
    [~, idx] = max(list(2, s+1:e));
    idx = s + idx - 1;
    peaks = [peaks; idx]; %#ok<AGROW> 
    s = idx + LEN;
end

% リスト生成
x = ctime + seconds(peaks / FPS);
y = diff(peaks / FPS);
v = list(2, peaks + 1)';
sz = numel(y);

% 外れ値の除去
y2 = y;
for i = 1:sz
    indices = i + [-3:-1 1:3];
    indices = indices(0 < indices & indices <= sz);
    score = mean(y2(indices)) - y(i);
    if score > 2.5
        warning('Outlier - date:%s, period:%f, score:%f', x(i), y(i), score);
        y(i) = NaN;
    end
end

% y は diff をとっている都合上 1 つ増やす必要がある
y = [y; NaN];


function output = cell2matdatetime(input)
% datetime に cell2mat が使えなかったので苦肉の策
output = NaT(sum(cellfun(@numel, input)), 1);
s = 0;
for i = 1:numel(input)
    c = input{i};
    sz = numel(c);
    output(s + 1:s + sz) = c;
    s = s + sz;
end
```

cellfun とか arrayfun とかいっぱいでてきてアレだが MATLAB っていうのはこういう言語である。

> ……
> 『オマエさ、そんなコード書くの……反則』

MATLAB はプロプライエタリで有料のプログラミング言語だけど、その分ドキュメントとか挙動とかが非常に洗練されいてるので、ああもうこういうとこがプロプライエタリ最高なんだよな。GNU Octave とは計算の洗練度が違う。探索とか苦手な処理もあるけど。

![好きな数値解析ソフトウェアについて聞く凛世](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78668/3ef3aa31-9ccb-0ec5-a48f-8b13505fb7bb.png)

</div></details>

# 計算結果と考察

まず最初に書いておくと、残念ながら凛世の自転周期を調べだしたのが放送開始からかなり時間が経ってからだったので、100 時間分データがあるわけではない。

> プロデューサーさま……
> こんな凛世を、お許しくださいますか?

## 凛世は少しずつ変化している

![回転したことによる凛世の変化](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78668/233688ee-c59e-24c9-cb23-004df8268d73.png)

※ 実際の開始時刻は 11:37 だったので、開始からの時間はもう少し短い。

実装が面倒だったのである時刻（5/4 1時ごろ）の凛世を常にリファレンスとして用いていたが、PSNR を見ると凛世が少しずつ変化していることが明らかになった。

![少しずつ大人になっていく凛世](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78668/853234cb-3504-f2b2-56a8-d6d76ce82e4a.png)

<font color="#89c3eb">実際のところ全体的な色味が少し変わってるっぽい。</font>

## 本当に自転周期はわかるのか

ある時刻での PSNR を見ると、ちゃんと心電図のように規則正しくピークが立っている。いや心電図とは全然違うけど。

![周期性が観測される凛世](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78668/98d9afc2-c382-e019-808e-23a548b1e8f4.png)

赤丸は検出されたピーク。およそ 1 分ごとにピークがあるので、目で見た感じと矛盾していなさそうだ。

## で、自転周期は変化しているんですか?

**1 日の中での変動も、長期的な変動もある。凛世は常に変わり続けている。**

### 1 日の中での変動

5 月 5 日（開始から 61〜85 時間）のグラフを示す。オレンジ色の線はガウス加重移動平均。

![こどもの日の凛世。凛世はこいのぼり好きそう。知らんけど。](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78668/7a37eca0-17a8-8673-783a-f1c918776c38.png)

かなりトゲトゲして見えるが、実際には 1 分に 1 回転しているので、結構ゆるやかに変化している。たぶん見た目だけで「今はやくなってるなー」と感じることはほとんどできないと思う。

といってもたった 6 周で自転周期が 3 秒以上速くなったり、1 周で 2 秒以上遅くなったりするなど、常にゆるやかに変化しているわけではない。

![凛世の変化の中では大きい方](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78668/8632293f-5b23-4b44-99e6-3786485a0d26.png)

![唐突にやる気を出した凛世](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78668/92b2f32c-a263-f9c7-4199-4cd3bad10653.png)

なお、どうやら「朝は元気で夜は元気がない」みたいな傾向はないらしい。

繊細な凛世のことだから 1 日の間で変動する重力の影響を受けたりするんじゃないかと思ったが、9.8 m/s<sup>2</sup> の重力加速度に対して 2 μm/s<sup>2</sup> くらいしか変動しないらしい [^gal]。<font color="#89c3eb">周りに人が通ったほうが影響受けそうだな。</font>
[^gal]: [地球の重力 - Wikipedia](https://ja.wikipedia.org/wiki/%E5%9C%B0%E7%90%83%E3%81%AE%E9%87%8D%E5%8A%9B#%E3%81%9D%E3%81%AE%E4%BB%96%E3%81%AE%E8%A6%81%E5%9B%A0)

### 長期的な変動

記録していた全期間のグラフを示す。オレンジ色の線はガウス加重移動平均。

![凛世が少しずつ大人になっていった記録](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78668/717ad1f2-54e0-b00f-b704-451d180bd856.png)

どうやら、だんだん遅くなっている期間と、だんだん速くなっている期間があるように見える。

はじめは「少しずつ速くなってる」「少しずつ遅くなってる」のどちらかだと思ったが、どうも速くなったり遅くなったりするらしい。きっと凛世の心もこの 100 時間でめまぐるしく変化しているのだろう。

### 最速凛世と最遅凛世

観測していた凛世の中で、最も元気だった凛世と最も落ち込んでいた凛世を比較してみた。

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">観測史上最速凛世と最遅凛世 <a href="https://t.co/1Im6fvzfVd">pic.twitter.com/1Im6fvzfVd</a></p>&mdash; イショティハドゥス (@Ishotihadus) <a href="https://twitter.com/Ishotihadus/status/1257947006128939008?ref_src=twsrc%5Etfw">May 6, 2020</a></blockquote><script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

## ヒストグラムをとる

ヒストグラムを見れば、全体的な自転周期の分布がある程度把握できる。

### 自転周期のヒストグラム

凛世が自転したときの周期の長さの分布を調べた。

![凛世の回転のヒストグラム。概念が意味不明すぎる](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78668/95309ded-8e3b-fdd2-9319-5e0c262793e3.png)

パッと見短い側に裾が伸びているが、およそ正規分布に従うらしい。標準偏差は 0.93 秒だった。

そこそこ綺麗な分布を描いているので、凛世を乱数生成器にした擬似乱数が作れるかもしれない。

### 自転周期の変化のヒストグラム

今の自転と次の自転の、周期の差（周期の回転数微分）をヒストグラムにした。x 軸の単位は「s/回」が正しい（なおすの面倒くさかった）。

![凛世の回転周期の速度のヒストグラム。およそ角加速度にあたる。さらに意味不明になってきたな。](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78668/4f4fe842-e71e-05d5-5493-d656c3c344f2.png)

こちらは、平均は 0.0019 秒とほぼ 0 秒だが、中央値は 0.033 秒と、少し偏った分布になっている。

**遅くなるときはゆっくり遅くなり、速くなるときは一気に速くなる**らしい。

ただ、平均がたった 0.0019 秒でも、（この平均がもし正確だったとして）100 時間積み重なれば 10 秒以上の差になるので、凛世は少しずつ遅くなっている、と結論してもいいのかもしれない。

# 今後の展望

## 角周波数の計測

凛世の角周波数がわかれば、離散的にしか計算できなかった凛世の回転が、より連続的に評価できそう。

どうやって実現するんだ。

## 凛世の「揺れ」

長時間の回転により、凛世は回転角度によって揺れてしまう現象が起きてしまっている。

これを検出できたらおもしろそう。凛世の自転周期とも関連がありそうだ。

# 結論

![こひめ](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78668/da65a882-830a-bc1a-427c-cacbffa0aded.png)

ありがとう凛世、あなたが 100 時間回転していたことを僕はきっと忘れない。

いつかこの「100 時間回った凛世」の実物が見られるといいですね。

<script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>