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

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

KT 2023-9-21

返信が遅れて大変申し訳ありません。

トレンド除去の妥当性や面積補正の方法などを考えておりました。

非常に有用なコード作成までご教示ありがとうございました。

请先登录，再进行评论。

Answer 2

Hiro Yoshino 2023-9-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2017621-#answer_1303876

在 MATLAB Online 中打开

回転角を考えて、全体のデータを回転させると良いかと思います。

元のデータ

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 个评论
显示 5更早的评论隐藏 5更早的评论

Hiro Yoshino 2023-9-8

在 MATLAB Online 中打开

今のベースで少し修正するとこんな感じになります。ご評価あればいただけると、ここからスタートして変更することもできるかと。

% データ作成

x = (0:36)';

y = 2*randn(size(x)) + 10*sin(x*2*pi/20) - 1.97*x + 116.08;

plot(x,y,'o');

[lfit, gof] = fit(x,y,'poly1') % 一次関数をデータにフィット

lfit =

Linear model Poly1: lfit(x) = p1*x + p2 Coefficients (with 95% confidence bounds): p1 = -2.142 (-2.358, -1.926) p2 = 119.8 (115.2, 124.3)

gof = struct with fields:

sse: 1.6722e+03 rsquare: 0.9205 dfe: 35 adjrsquare: 0.9182 rmse: 6.9122

plot(lfit,x,y);

coeffs = coeffvalues(lfit) % モデルの係数を取得

coeffs = 1×2

-2.1421 119.7550

%theta = 2 * (pi/2 - atan(coeffs(1))) % 回転角を計算

theta = -atan(coeffs(1)) % 回転角を計算

theta = 1.1340

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,'bo-',x0y0(1,:),x0y0(2,:),'ro-'); % プロット

Hiro Yoshino 2023-9-8

ベースライン補正で何となくどんなことをされたいのかは分かりました。

(校正のことなのかなと理解しました)

その上で、ベースライン補正で利用されるアルゴリズムなどあれば、ご紹介いただけると取り掛かり易いです。

KT 2023-9-8

在 MATLAB Online 中打开

数学的に議論したいのですが、知識不足で他文献もないため、分かりにくく大変申し訳ありません。

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値が補正できればと考えています。

请先登录，再进行评论。