5

私は、減衰された高調波振動を含むデータセットに指数曲線を当てようとしています。データは下図のように正弦波振動が多くの周波数が含まれているという意味で少し複雑です:MATLABの減衰された高調波振動データに指数曲線をフィットさせる方法は?

enter image description here

私はデータで減衰率を見つける必要があります。私が使っている方法はhereです。それがどのように機能するかは、定常状態の値より上のy値のログをとり、次に使用します:

lsqlin(A,y1(:),-A,-y1(:),[],[],[],[],[],optimset('algorithm','active-set','display','off')) 

これを適合させます。

しかし、これは次のようなデータに収まる結果: enter image description here

私はそれが平均を取ったので、明らかに動作しませんでした線形回帰適合を使用してみました。私はまた、RANSACが山の近くにもっと多くのデータがあると考えました。線形回帰より少し良く機能しましたが、間違った領域に多くの点が存在するため、この方法には欠陥があります。

このデータのピークにちょうど合う良い方法を知っている人はいますか?

現在、私は500のデータ点を10の異なる領域に分け、各領域で最大の値を見つけることを考えています。最後に、上記の指数関数フィッティング法を使ってフィットできる50ポイントを得てください。この方法についてどう思いますか?

答えて

1

私はすべての人に、可能性のある潜在的な解決策のアップデートを提供したいと思います。先に述べたように、データは正弦波の周波数が変化すると複雑になります。そのため特定の方法が機能しない可能性があります。下記の方法は、データと関連する周波数に応じて良好にすることができます。それは言われていると

y = 290 + b*e^-(c*x) 

を、のに飛び込むてみましょう:私の場合は

y = average + b*e^-(c*x) 

、平均は私たちが持っている290です:

まず第一に、私は、データの形式を持っていることを前提としてい私が試し異なる方法:

findpeaks()メソッド

これはAlexanderBüseが提案した方法です。これはほとんどのデータにとってはかなり良い方法ですが、私のデータでは、複数の正弦波の周波数があるので、間違ったピークが得られます。赤いxはピークを示します。あなたがピークで、あなたのデータのほとんどを持っている場合

% Find Peaks Method 
[max_num,max_ind] = findpeaks(y(ind)); 
plot(max_ind,max_num,'x','Color','r'); hold on; 
x1 = max_ind; 
y1 = log(max_num-290); 
coeffs = polyfit(x1,y1,1) 
b = exp(coeffs(2)); 
c = coeffs(1); 

enter image description here

RANSAC

RANSACが良いです。私の場合、複数の周波数のために、頂上近くにピークが存在することがわかります。しかし、私のデータの問題は、すべてのデータセットがこのようなものではないということです。したがって、時にはうまくいった。

% RANSAC Method 
ind = (y > avg); 
x1 = x(ind); 
y1 = log(y(ind) - avg); 
iterNum = 300; 
thDist = 0.5; 
thInlrRatio = .1; 
[t,r] = ransac([x1;y1'],iterNum,thDist,thInlrRatio); 
k1 = -tan(t); 
b1 = r/cos(t); 
% plot(x1,k1*x1+b1,'r'); hold on; 
b = exp(b1); 
c = k1; 

enter image description here

Lsqlin方法

この方法はhereを用いるものです。 Lsqlinを使用してシステムを制約します。しかし、それは真ん中のデータを無視しているようです。あなたのデータセットに応じて、元の投稿の人物と同じようにうまくいくかもしれません。

% Lsqlin Method 
avg = 290; 
ind = (y > avg); 
x1 = x(ind); 
y1 = log(y(ind) - avg); 
A = [ones(numel(x1),1),x1(:)]*1.00; 
coeffs = lsqlin(A,y1(:),-A,-y1(:),[],[],[],[],[],optimset('algorithm','active-set','display','off')); 
b = exp(coeffs(2)); 
c = coeffs(1); 

enter image description here

が期間

これはでピークを探すには、私は、各地域でピークを取得し、私の記事で述べた方法です。この方法はかなりうまくいくので、私のデータが実際に完全な指数関数的フィットにならないかもしれないことに気付きました。最初は大きなピークにフィットすることができないことがわかります。最初の150のデータポイントを使用し、定常状態のデータポイントを無視するだけで、これを少し改善することができました。ここで私は25のデータポイントごとにピークを見つけました。最初の150個のデータポイントを使用して Using all 500 data points

:すべての500個のデータポイントを使用して

% Incremental Method 2 Unknowns 
x1 = []; 
y1 = []; 
max_num=[]; 
max_ind=[]; 
incr = 25; 
for i=1:floor(size(y,1)/incr) 
    [max_num(end+1),max_ind(end+1)] = max(y(1+incr*(i-1):incr*i)); 
    max_ind(end) = max_ind(end) + incr*(i-1); 
    if max_num(end) > avg 
     x1(end+1) = max_ind(end); 
     y1(end+1) = log(max_num(end)-290); 
    end 
end 
plot(max_ind,max_num,'x','Color','r'); hold on; 
coeffs = polyfit(x1,y1,1) 
b = exp(coeffs(2)); 
c = coeffs(1); 

私はそれをしたいので enter image description here

は、Bは期間中のピークを探すには、

を拘束しました最初のピークから始めて、私はb値を制限しました。私はシステムがy=290+b*e^-c*xであることを知っていて、私はそのようなことを拘束します。b=y(1)-290。そうすることで、私はちょうどc=(log(y-290)-logb)/xのCのために解決する必要があります。次に、cの平均値または中央値を取ることができます。このメソッドはかなり良いですが、それは最後の値にも合わないが、変更が最小であるため、それほど大きな取引ではない。

% Incremental Method 1 Unknown (b is constrained y(1)-290 = b) 
b = y(1) - 290; 
c = []; 
max_num=[]; 
max_ind=[]; 
incr = 25; 
for i=1:floor(size(y,1)/incr) 
    [max_num(end+1),max_ind(end+1)] = max(y(1+incr*(i-1):incr*i)); 
    max_ind(end) = max_ind(end) + incr*(i-1); 
    if max_num(end) > avg 
     c(end+1) = (log(max_num(end)-290)-log(b))/max_ind(end); 
    end 
end 
c = mean(c); % Or median(c) works just as good 

ここで私はすべての25個のデータポイントのピークを取り、その後、C enter image description hereここ

の平均を取る、私はすべての25個のデータポイントのピークを取り、その後、C enter image description here

の中央値を取ります

は、ここで私はすべての10個のデータポイントのピークを取り、その後、C enter image description here

0

主要な目的がフィットからダンピングパラメータを抽出することである場合は、ダンピングされたサインカーブをデータに直接フィッティングすることを検討してください。 (カーブフィッティングツールで作成)このような何か:

[xData, yData] = prepareCurveData(x, y); 
ft = fittype('a + sin(b*x - c).*exp(d*x)', 'independent', 'x', 'dependent', 'y'); 
opts = fitoptions('Method', 'NonlinearLeastSquares'); 
opts.Display = 'Off'; 
opts.StartPoint = [1 0.285116122712545 0.805911873245316 0.63235924622541]; 
[fitresult, gof] = fit(xData, yData, ft, opts); 
plot(fitresult, xData, yData); 

は特に、あなたの例のデータの一部以来、本当に(ノイズ以上)の興味深い地域の多くのデータポイントを持っていません。

しかし、実際に実験データの最大値に直接フィットする必要がある場合は、findpeaks関数を使用して最大値のみを選択してからフィットさせることができます。 MinPeakProminenceパラメータを少し再生して、必要に応じて調整することができます。

+0

おかげでアレクサンダーの平均を取ります。したがって、このデータのトリッキーなことは、それが単なる正弦波の周波数ではないことです。したがって、findpeaksを使用すると、実際には領域の最大値ではない波のピークが得られます。オリジナルのポストを実際の信号で更新し、ポイントを接続しました。私はまだカーブフィッティングツールボックス(コンピュータ上にはない)にフィッティングを試みてきましたが、私が学校にいるときに試してみましょう。 –