Effacer les filtres
Effacer les filtres

一次式によるベースライン補正について

68 vues (au cours des 30 derniers jours)
KT
KT le 6 Sep 2023
Commenté : KT le 21 Sep 2023
グラフで得られているx軸0-36までの平均値は変えずに、
最小二乗近似法で得られたオレンジ点線の直線(y = -1.97x+116.08)と「同じ傾きで右肩上がりの直線」を求め、その直線に従いy軸を変化させるにはどのようにしたらいいでしょうか?
もしくは、
x軸が大きくなるにつれて、y軸が高いところを低く、低いところを高く、平均値は変えず全体的に右肩上がりに補正できる方法はないでしょうか?
スペクトルではないので、ベースライン補正なのか分かりませんが、簡易的な一次式で補正できれば考えています。

Réponse acceptée

Hiro Yoshino
Hiro Yoshino le 9 Sep 2023
面積の補正が必用ということで、いただいたコードを少し変えてみました。
利便性を考慮して最初のデータでも動作するようにしました。y が before, y2 が after です。
% Original data
x = (1:36)';
y = [0;134.7407;205.9038;210.5938;148.8053;101.2817;94.7168;89.2617;92.2141;100.3457;103.1461;98.5113;72.87;72.7713;71.103;60.4682;53.5974;58.5912;51.7347;51.1305;50.0097;44.1407;49.0088;45.1054;52.8028;49.5283;68.3816;58.6178;75.5166;67.1995;82.3074;84.7901;103.3022;165.2237;0;0];
plot(x,y);
title("元のデータ");
トレンド除去
idx = find(y>0); % y>0 のインデックス
idx0 = find(y<=0); % y<=0 のインデックス
d = detrend(y(idx),1); % y>0 だけに限定してトレンド除去後を計算
% 全てのデータポイントでトレンド除去を考慮
y2 = zeros(size(y));
y2(idx) = d;
y2(idx0) = y(idx0); % y<=0 の部分には元のデータを保持
トレンド除去後の面積補正
% 面積の補正量
sComp = mean(y(idx)-y2(idx))
sComp = 86.9007
y2(idx) = y2(idx) + sComp; % 面積補正後
% 比較
figure;
tiledlayout(2,2);
nexttile
scatter(x,y,'b');
lsline;
title("Before と LS線")
nexttile
scatter(x,y2,'m');
lsline
title("After と LS線, '0'のためLS線が下がる");
nexttile
plot(x,y,x,y2);
title("2つの比較");
nexttile
plot(x,y-y2);
title(sprintf("面積が等しい証明: mean(y-y2)=%1.3f",mean(y-y2)));
legend("y-y2");
  1 commentaire
KT
KT le 21 Sep 2023
返信が遅れて大変申し訳ありません。
トレンド除去の妥当性や面積補正の方法などを考えておりました。
非常に有用なコード作成までご教示ありがとうございました。

Connectez-vous pour commenter.

Plus de réponses (1)

Hiro Yoshino
Hiro Yoshino le 7 Sep 2023
回転角を考えて、全体のデータを回転させると良いかと思います。
元のデータ
x = 0:0.1:36;
y = -1.97*x + 116.08;
plot(x,y,'b:');
回転させます
theta = 2 * (pi/2 - atan(-1.97)); % 回転角を計算
R = [cos(theta) -sin(theta); sin(theta) cos(theta)]; % 回転行列
center = repmat([mean(x); mean(y)], 1, length(x)); % 中心座標
s = [x;y] - center; % 原点に移動
s0 = R*s; % 回転
x0y0 = s0 + center; % 元に戻す
plot(x,y,'b:',x0y0(1,:),x0y0(2,:),'r:');
  7 commentaires
Hiro Yoshino
Hiro Yoshino le 8 Sep 2023
ベースライン補正 で何となくどんなことをされたいのかは分かりました。
(校正のことなのかなと理解しました)
その上で、ベースライン補正で利用されるアルゴリズムなどあれば、ご紹介いただけると取り掛かり易いです。
KT
KT le 8 Sep 2023
数学的に議論したいのですが、知識不足で他文献もないため、分かりにくく大変申し訳ありません。
yデータに対し一次式でのベースライン補正を考えたときに、detrend(x,1)が使えるのではないかと考えました。
% x, yデータ
x = (1:36)';
y = [0;134.7407;205.9038;210.5938;148.8053;101.2817;94.7168;89.2617;92.2141;100.3457;103.1461;98.5113;72.87;72.7713;71.103;60.4682;53.5974;58.5912;51.7347;51.1305;50.0097;44.1407;49.0088;45.1054;52.8028;49.5283;68.3816;58.6178;75.5166;67.1995;82.3074;84.7901;103.3022;165.2237;0;0]
y = 36×1
0 134.7407 205.9038 210.5938 148.8053 101.2817 94.7168 89.2617 92.2141 100.3457
% 線形近似 データの傾向として右肩下がりなのか確認のため
X = [ones(length(x),1) x ];
b = X\y
b = 2×1
116.0829 -1.9689
yCalc = X*b;
plot(x,y,'-')
hold on
plot(x,yCalc,'--')
% 線形トレンド
D = detrend(y,1);
plot(x,D)
plot(x,y-D,":k")
legend("Input Data","Detrended Data","Trend")
y値は正のため、0を含めてしまうとトレンド除去の際にy値0が負に補正されてしまうため、以下では0は除去することを考えました。
x = (1:33)';
y = [134.7407;205.9038;210.5938;148.8053;101.2817;94.7168;89.2617;92.2141;100.3457;103.1461;98.5113;72.87;72.7713;71.103;60.4682;53.5974;58.5912;51.7347;51.1305;50.0097;44.1407;49.0088;45.1054;52.8028;49.5283;68.3816;58.6178;75.5166;67.1995;82.3074;84.7901;103.3022;165.2237]
y = 33×1
134.7407 205.9038 210.5938 148.8053 101.2817 94.7168 89.2617 92.2141 100.3457 103.1461
% 線形近似 データの傾向として右肩下がりなのか確認のため
X = [ones(length(x),1) x ];
b = X\y
b = 2×1
122.2141 -2.0773
yCalc = X*b;
plot(x,y,'-')
hold on
plot(x,yCalc,'--')
% 線形トレンド
D = detrend(y,1);
plot(x,D)
plot(x,y-D,":k")
legend("Input Data","Detrended Data","Trend")
ただし、これでは積分値が変化することになります。これを変化させないようにトレンド除去する方法はありますでしょうか?
>「右肩上がりの直線に変換し」 - 変換する方法 (回転) を示しました。この部分はこの方法で大丈夫ですか?
よく分かりました。回転を利用することで、近似曲線を右肩上がりに変換できるので納得しました。ただ、変換した近似曲線にしたがってy値が補正できればと考えています。

Connectez-vous pour commenter.

Produits


Version

R2021a

Community Treasure Hunt

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

Start Hunting!