センシングした波形の立ち上がりからピーク（最大値）に到達するまでの時間を知りたい

Question

実現したいこと

ESP32を使ってセンシングしたのちに、波形データの立ち上がりからピーク（最大値）に達するまでの時間を算出したいと考えています（ESPのプログラム内で）

立ち上がりとピーク位置については、下記グラフの青の波形に示したピンク矢印の位置をイメージしています（低周波波形はピークが頭打ちのような形になるので定義が難しいですが、傾きがほぼゼロとなり、最大値をとる最初の点になるかと思います）

ベースが毎回一定でない（波形ごとに多少変動する可能性がある）、波形の周期が一定でないこともあり、どのようにアルゴリズム（プログラム）を組めばいいか苦慮しております。
アドバイスいただけましたらありがたいです。どうぞよろしくお願いいたします。

※生データはばらつきが大きいため、メディアン処理、移動平均処理を施しています。
今回時間データは保存していないため、便宜的にデータ取得順に１，２，３、、、とすることとしています。

処理後データ抜粋（上段：低周波波形データ（１波形ぶん切り出し）、下段：高周波波形データ（２波形ぶん切り出し）

[64,64,64,64,64,64,64,64,64,64,64,64,64,64.1,64.1,64.1,64.1,64.1,64.25,64.35,64.4,64.5,64.6,64.75,64.75,64.85,64.95,65.05,65.25,65.45,65.7,65.8,66,66.1,66.1,66.3,66.5,66.65,66.7,66.85,66.95,67.05,67.15,67.2,67.45,67.55,67.65,67.8,67.9,67.95,67.9,68.05,68.05,68.15,68.4,68.5,68.7,68.75,68.95,69.05,69.2,69.35,69.55,69.75,69.9,70,70.05,70.1,70.2,70.35,70.6,70.75,70.95,71.05,71.2,71.35,71.35,71.55,71.55,71.6,71.7,71.75,71.75,71.8,71.8,72,72.4,72.55,72.65,72.7,72.65,72.75,72.8,72.9,73,73.05,73.15,73.2,73.3,73.4,73.5,73.6,73.6,73.75,73.95,74,73.9,74,74,74.15,74.3,74.5,74.65,74.75,74.85,74.95,75.05,75.2,75.3,75.3,75.4,75.5,75.7,75.9,75.9,76.05,76.1,76.35,76.55,76.55,76.7,76.8,76.8,76.8,76.8,76.8,77,77.2,77.35,77.45,77.45,77.45,77.65,77.6,77.8,77.9,78,78,78.1,78.2,78.3,78.3,78.4,78.5,78.7,78.9,78.8,78.8,78.85,79.05,79.25,79.45,79.45,79.4,79.4,79.4,79.5,79.5,79.5,79.6,79.6,79.6,79.7,79.8,79.8,79.6,79.7,79.7,79.7,79.7,79.7,79.7,79.7,79.8,79.8,79.8,79.8,79.8,79.8,79.8,79.8,79.8,79.8,79.8,79.8,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,80,79.85,79.45,78.85,78.1,77.3,76.5,75.7,74.9,74.1,73.3,72.5,71.7,70.9,70.1,69.3,68.5,67.7,66.9,66.1,65.3,64.65,64.25,64.05,64,64,64,64]

[64,64,64,64,64,64,64,64,64.1,64.35,64.95,65.75,66.55,67.35,68.15,68.95,69.8,70.75,71.7,72.7,73.7,74.65,75.45,76.25,76.85,77.15,77.25,77.3,77.2,76.95,76.35,75.55,74.75,73.95,73.15,72.35,71.5,70.55,69.6,68.6,67.6,66.65,65.85,65.05,64.45,64.15,64.05,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64,64.2,64.7,65.5,66.3,67.1,67.95,68.9,69.9,70.9,71.8,72.8,73.8,74.8,75.6,76.4,76.8,77,77.1,77.1,77.1,76.9,76.4,75.6,74.8,74,73.15,72.2,71.2,70.2,69.3,68.3,67.3,66.3,65.5,64.7,64.3,64.1,64,64,64,64,64,64,64,64,64,64,64]

試したことと問題点

ベースの値（今回は64くらい）が１波形ごとに多少変動することを考慮し、生データを微分してゼロになる箇所から立ち上がりおよびピーク位置を割り出すことを考えました。高周波の波形ではこの方法でうまく解析できそうですが、低周波波形では、波形の立ち上がりにかけてデータが多少乱れることと、なだらかに上昇していくために、微分による立ち上がりとピークの判別が難しいということがネックになっています。

※グラフ横軸は、データ取得順に１，２，３、、、としており、これを便宜的に時刻としています。

Comments

[y_waiwai]

で、現状の処理ではどういう問題があるんでしょうか。
そしてそれをどうしたいという話なんでしょう

[Uka]

y_waiwaiさん
コメントありがとうございます。要点がわかりづらく申し訳ございません。見出しを微修正させていただきました。現状の処理と問題については見出しの「試したことと問題点」に示したとおりです。自分が考えた微分の方法ではとくにピークの判別が難しいと思ったため、微分の方法以外にもほかにもっといいアルゴリズムがあればご教示いただければありがたいです。

[y_waiwai]

波形データの立ち上がりとピークの時間、というやつが意味不明です。
あなたにとって、これはどういう意味なんでしょうか。
問題が見えません。

繰り返しますが、現状の問題は何でしょうか

[ozwk]

(「時間」は時刻の差ですからある一点を指したいときは「時刻」です)

ピーク及び立ち上がりそれぞれの時刻とは何でしょうか？

回答者が適当に決めて回答しても「それは私が思っていたものと違う」と言われかねないので回答しにくいです

立ち上がり時刻とはピークのx%を通過した時刻とする
とかある程度明確な説明をください
(それが決められないからこの質問なのかもしれませんが…)

[SaitoAtsushi]

質問者は「立ち上がり」「ピーク」が何であるか定義してください。
と言いたいところですが、おそらくそれをどう定義してよいかわからないということなんでしょう。
まずはきちんとした定義でなくてもよいですから、左の波形のどこが立ち上がりでどこがピークだと考えているのかを図示したほうがよろしいですね。

[thkana]

プログラミングは全然関係ないですよね。まぁいいですけど。
で、「低周波」の方はどこを「ピーク」とみなすのが目的に対して妥当なのですか?
(その波形で「ピーク」を論じること自体に疑問を感じてしまうのです)

[Uka]

>ozwkさん、SaitoAtsushiさん、thkanaさん
コメントいただきありがとうございます。とくに立ち上がりとピークについて説明足らずで大変失礼しました。説明文およびグラフを修正させていただきました。
まず知りたいことは、「立ち上がりからピークに到達するまでの時間」です。
立ち上がりとピークの定義については、SaitoAtsushiさんのおっしゃるように一義的に定義できないところもあるため、イメージとしてグラフに矢印で図示させていただきました。
立ち上がりについては、たとえばベースから５％上昇したときの時刻、
ピークについては、最大値となる時刻（低周波波形については最大値をとる最初の時刻）
などを考えています。

[jbpb0]

> 生データはばらつきが大きいため、メディアン処理、移動平均処理を施しています。

移動平均だと波形が歪んでピークがなまりやすい(高周波成分が変わりやすい)ので、Savitzky-Golay法の方がいいと思います
平滑化と微分を同時に行うこともできます
https://izadori.net/math-savitzky-golay/
https://datachemeng.com/wp-content/uploads/preprocessspectratimeseriesdata.pdf

[Uka]

>jbpb0さん、
アドバイスいただきありがとうございます。データ前処理の方法について大変参考になります。
こちらの方法で試してみたいと思います。

[thkana]

「低い周波数」なのだから、「高い周波数」との比に応じてサンプリング周波数も同じだけ下げてやれば同じような処理で同じような結果が得られそうな気がしますが、そういうわけにはいかないのでしょうか。

[jbpb0]

> 立ち上がりについては、たとえばベースから５％上昇したときの時刻

変曲点を「立ち上がり」とすることがあります
https://www.mlit.go.jp/river/shinngikai_blog/shaseishin/kasenbunkakai/shouiinkai/kihonhoushin/060719/pdf/add1.pdf
http://www.mekatoro.net/digianaecatalog/valco-sougou/book/valco-sougou-P0241.pdf

微分値のノイズを減らす必要があります
https://hanafsky.com/tips/savgol/

質問者さんのケースが該当するか分かりませんが、ご参考までに

[Uka]

>thkanaさん、
コメントありがとうございます。できれば測定と解析（立ち上がりからピーク値までの時間）をリアルタイムで行いたいと考えていますが、高周波か低周波どちらの波形になるかは成り行きになってしまうので、測定中にサンプリング周波数を変えることは難しいと思っています。

[thkana]

リアルタイムという言葉は、少なくとも「リアルタイムシステム」という言葉は「入力からある時間内に処理を終える系」という意味で使います。システムの目的によっては、それがマイクロ秒なのか数時間なのか、定義次第(まぁ普通は秒単位まででしょうが)。逆に言えば、制限時間を明らかにしなければ、何も言っていないと同じことになってしまいます。
なので、何ミリ秒以内とか、あるいはフィルタ方法を考える観点からは何サンプル周期以内に検出できなければいけない、という形で定義するのもありでしょう。

> 測定中にサンプリング周波数を変えることは難しい
フィルタにとってサンプリング周波数が下がっていればいいわけで、つまりフィルタに与えるデータを間引きして処理してやればよいのではないかと思いますがいかが? で、フィルタは高周波低周波どちらも回しておいて、微分値かなにかでどちらの結果を採用するか決める、と。もちろん、低周波時の応答時間は遅くなるでしょうが、でも信号の変化が緩いんだから仕方ないです。

Answer 1

生データはばらつきが大きいため、メディアン処理、移動平均処理を施しています。

移動平均だと波形が歪んでピークがなまりやすい(高周波成分が変わりやすい)ので、Savitzky-Golay法の方がいいと思います
平滑化と微分を同時に行うこともできます

【C++】【Eigen】Savitzky-Golay法を実装する
の「スムージングの例」の「1次微分」のグラフを見ると、単純な1次微分よりもノイズが抑えられることが分かります
(平滑化を強くすれば、もっとノイズは減ります)

時系列及び波形データの平滑化３手法(smoothing)
の「2. Savitzky-Golay filter」の「2.2 １階微分」に、pythonのコード例があります

「scipy.signal.savgol_filter()」の二つ目の引数(上記Webページのコード例では「フレーム長」と書かれてる)の数値(奇数のみ)を大きくすると、平滑化が強くなりノイズがより減りますが、微分後データの形状(急激な立ち上がり・立ち下がり等)が徐々に崩れるので、そのバランスで数値を決めます

　
なお、データが等間隔ではない場合は、おそらく「scipy.signal.savgol_filter()」は使えないと思います
その場合は、
pythonで曲線の各点の傾きを表すグラフを作成したいです
の私の回答の【追記】のコードを使えば計算できます

Comments

[jbpb0]

質問に記載のデータでやってみました

高周波では

x = scipy.signal.savgol_filter(wave, 5, 2, deriv=1)

で計算した微分値が、
・0.1くらいを超える時刻を立ち上がり
・(ほぼ)0.0に戻る時刻をピーク
とすればいいようです

低周波ではノイズが多いので、平滑化をうんと強くして

x = scipy.signal.savgol_filter(wave, 51, 2, deriv=1)

で計算した微分値が、
・0.05くらいを超える時刻を立ち上がり
・0.01くらいより下がる時刻をピーク
とすれば何とかなりそうに思いますが、ピークの時刻は

> 傾きがほぼゼロとなり、最大値をとる最初の点

よりもちょっと前(グラフの左側)です
(かなり最大値には近いですが)

微分値の絶対値(ほぼ)0.0がある時間以上連続したら、的な条件も併用した方がいいかもしれません

他のデータに対しても同じ判定値でうまくいくのかは、データ形状がどれくらい似てるかによります
平滑化をうんと強くすると、Savitzky-Golay法でも微分後データの形状は崩れるので、平滑化を強くしてる低周波の方がデータ形状の変動の影響を受けやすいと思います

Answer 2

信号がLowからHighになるとき、理想的な状況では無限大の時間が掛かりますので、Highになった時刻というのは定義できません。
ですので一般的に使われている、信号がLowからHighになる速さを測るための概念、「立ち上がり時間」と呼ばれているものですが、これはふつう10%から90%の間の遷移にかかる時間で定義されています。
今回の場合、Low→Highの遷移は一定以上の変化があったときに認識するとして、あとはLowとHighの電圧を適宜最大値なり99%なりで定義してやれば10%から90%までの時間が求められるでしょう。