交互に正負反転して現​れるパターンの一方の​みを残したい(削除し​たい)

[Matlab初心者です]
「順回転→静止→逆回転」を繰り返す振り子の順回転時の加速度データのみ得ようと試みています.
データは図(上)のように,逆方向へのスパイクに続いて回転方向への加速度が現れるパターンとなっています(縦軸:加速度,横軸:時間).
スパイクを含めて,順(正)方向への回転時の加速度のみのデータ図(下)にする良い方法があればご教示願いたいです.
直近で,自身で下述の方法
⓪加速度±5を0に置き換え,静止時を0にする
①findpeaksで順回転時のピークサーチ,peak Widthを取得
②ピークの平均的な間隔(avet (s))の取得
③ピークから avet / 2 前タイミング(delRef)を取得
④ delRefから± peak Width (s)を0に置き換える
を試したのですが,ディテクトされるピークの位置が山ごとに違うので,どうしても削除しきれなかったり削り過ぎてしまったりします.
浅慮かもしれませんが,今後似たようなデータを多量に取得する予定なので,なにか良い方法をご存じな方がいましたらご教示ください.
回転の加速度

2 comentarios

Hernia Baby
Hernia Baby el 17 de Jun. de 2022
Editada: Hernia Baby el 17 de Jun. de 2022
すみません。数点教えてください。
①findpeaksで順回転時のピークサーチ,peak Widthを取得
 順回転は負の方向ですか?
 peaks Widthとは何でしょうか?
②ピークの平均的な間隔(avet (s))の取得
 これはピーク間の横軸距離の平均という意味ですね
③ピークから avet / 2 前タイミング(delRef)を取得
④ delRefから± peak Width (s)を0に置き換える
  xを横軸とすると各ピークに対して±Δ(x/2)だけ0に変更するであってますか?
kozo
kozo el 17 de Jun. de 2022
ご丁寧にご質問ありがとうございます.
返答させて頂きます.
①順回転は,正の方向です.
 peak Widthは,ピークの幅(プロミネンスの半分を基準としたピーク幅)のことを示したつもりでした.
②おっしゃるとおりです.ピークとピークの横軸距離(時間)です.
④いえ,おおよそ逆回転の各ピークのタイミング±Δ(ピークの幅/2)をゼロにしました.
 書き忘れておりましたが,ピーク幅を取得した際に平均ピーク幅を算出しました.
なので修正しますと,
①findpeaksで順回転時のピークサーチ
②ピークの幅(プロミネンスの半分を基準としたピーク幅)を取得し,平均のピーク幅を算出.
 ピークのタイミングを取得し,ピーク間の横軸距離(avet (s))の平均を算出.
③逆回転のピークは,順回転のピークとその次のピークの半分あたりに現れるかなと考え,順回転の各ピークのタイミングからavetの1/2を差し引いたもの(delRef)を作成.
④delRef±Δ(ピークの幅/2)をゼロに置き換える.
を行いました.
返信させて頂く過程で...
①負のピークをサーチ
②負の各ピークのタイミング±Δ(ピークの幅/2)をゼロに置き換える
を行えば,簡単にほぼ同じ事ができる気づいてしまいました(汗).
しかし,実際に試してみたのですが,やはりピークの位置が山ごとに違うので,どうしても削除しきれなかったり削り過ぎてしまったりするのが現状です...

Iniciar sesión para comentar.

Respuestas (2)

Hernia Baby
Hernia Baby el 17 de Jun. de 2022

2 votos

■データの準備
勝手ながら画像からデータを作成しました。
なのでこの工程は無視してください。
clear,clc;
load data.mat;
x = Data(:,1)';
y = Data(:,2)';
figure
hold on
plot(x,y,'o-','Color', [.4,.4,.4]);
x1 = linspace(0,x(end),1e3)';
y1 = interp1(x,y,x1);
plot(x1,y1,'r:')
hold off
■加速度±5を0に置き換え,静止時を0にする
加速度が大きい気もしますが、5を閾値とします
y1(abs(y1)<=5) = 0;
figure
hold on
plot(x,y)
plot(x1,y1)
hold off
■findpeaksで順回転時のピークサーチ,peak Widthを取得
値が20以上のものからピークサーチします
[pks,locs,w,p] = findpeaks(y,x,'MinPeakHeight',20);
% 確認用
figure
findpeaks(y1,x1,'MinPeakHeight',20,'Annotate','extents')
■ピークの平均的な間隔(avet (s))の取得
w_mean = mean(w)
w_mean = 0.3817
■ピークから avet / 2 前タイミング(delRef)を取得
論理値をとります
for ii = 1: length(locs)
idx(:,ii) = double((x1(:) <= locs(ii) + w_mean/2) & (x1(:) >= locs(ii) - w_mean/2));
end
idx_any = any(idx,2);
■ delRefから± peak Width (s)を0に置き換える
y2 = y1;
y2(~idx_any) = 0;
■確認
figure
plot(x,y,'-','Color',[.6 .6 .6]);
hold on
plot(x1,y2,'r')
hold off

3 comentarios

Atsushi Ueno
Atsushi Ueno el 17 de Jun. de 2022
画像からデータを起こす謎の技術すごい!
Hernia Baby
Hernia Baby el 17 de Jun. de 2022
ふぁいるえくすちぇんじ〜
kozo
kozo el 18 de Jun. de 2022
ご丁寧に記して頂き感謝です.
勉強になります.
実際試してみて,④delRefを基準に値を0に置き換える際,± peak Width (s)をX0.6~0.8程度に調整すると綺麗に逆回転時の値を削除することができ現状の課題は解決することができました.
ご教示ありがとうございます.

Iniciar sesión para comentar.

Atsushi Ueno
Atsushi Ueno el 17 de Jun. de 2022

1 voto

「パルス幅(pulsewidth関数)でスパイクを捉える作戦」を実行してみました。
fs = 100; f = 1; w = 2 * pi * f; t = 0:1/fs:4; % サンプルデータの時刻t[sec]
y = (sin(w*t) > 0.7) * 50 - (sin(w*t) < -0.7) * 50 + (rand(size(t)) - 0.5) * 20;
y(12:fs:end) = -30 - rand(size(y(12:fs:end))) * 30; % サンプルデータのy座標
plot(t,y); % サンプルデータ。ちょっとキレイ過ぎるかな
y2 = y;
y2(abs(y-mean(y2)) < 0.5*std(y2)) = 0; % 変化の小さい数値を0にする
[posw,posinicrs,posfincrs] = pulsewidth(max(0,y2),fs); % 順方向(0以上の正極性)パルス幅
[negw,neginicrs,negfincrs] = pulsewidth(y2,fs,'Polarity','negative'); % 逆方向とスパイク(負極性)パルス幅
posfincrs = floor(posfincrs(posw > 0.1) * fs); % 順方向パルスの最終遷移時刻を得る
neginicrs = floor(neginicrs(negw < 0.05) * fs); % スパイクの初回遷移時刻を得る
mask = zeros(size(y2)); % スパイクの初回遷移時刻⇒順方向パルスの最終遷移時刻のマスクを作成
for i = 1:numel(posfincrs)
mask(1,neginicrs(i):posfincrs(i)) = 1;
end
figure
plot(t,y,'-','Color',[.6 .6 .6]);
hold on
plot(t,y.*mask,'r','LineWidth',2)
hold off

1 comentario

kozo
kozo el 18 de Jun. de 2022
Editada: kozo el 18 de Jun. de 2022
ご回答頂きありがとうございます.
密かに「パルスに形が似てるから近似できないか?,正のパルスとその周辺を残したら...?」などと考えつつ,着手までの前提知識が皆無なので喉元に留めていたところ,このようなご回答頂き大変感動しています.
「パルス幅(pulsewidth関数)でスパイクを捉える作戦」大成功しました.
凄くシンプルに解決に向かえました.
柔軟な発想,とても勉強になります.
ありがとうございます.

Iniciar sesión para comentar.

Preguntada:

el 17 de Jun. de 2022

Editada:

el 18 de Jun. de 2022

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!