How to find and fit the objective function for a damped oscillation?
19 ビュー (過去 30 日間)
古いコメントを表示
Piyumal Samarathunga
2022 年 11 月 28 日
コメント済み: Star Strider
2023 年 2 月 16 日
The excel file contains force data of a damped oscillation system. In order to find the natural frequency of the system I am trying to fit the graph with an objective function as follows.
I have defined the function of the osillation.
function F = forcesim(t,A,omega,phi,C,K)
F=A*cos(omega*t+phi)*exp(-t*C) + K;
end
I have tried with a code to find the estimated frequency as follows.
dataset = readtable ("F300.xlsx");
Force = dataset.Force2;
Time = dataset.Time2;
plot(Time,Force)
xlabel("Time (s)")
ylabel("Force (N)")
grid on
ind = Time>0.3 & Time<1.3;
[pks,locs] = findpeaks(Force(ind),Time(ind),'MinPeakProminence',1);
peaks = [pks locs]
% visualizse the selected peaks and label them.
hold on
plot(locs,pks,'v')
text(locs+.02,pks,num2str((1:numel(pks))'))
hold off
% calculate the approximate frequency
meanCycle = length(pks)/(locs(end)-locs(1))
Hence, that frequency is an estimated one, and therefore I am planning to find the frequency of the graph by matching the objective function mentioned above.
It would be really appreciated if you could tell me how to do that and find the frequency ω.
0 件のコメント
採用された回答
Star Strider
2022 年 11 月 29 日
T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1211598/F300.xlsx')
x = T1.Time2;
y = T1.Force2;
start = find(y > 1, 1, 'first') % Find Offset
y = y(start:end);
x = x(start:end);
y = detrend(y,1); % Remove Linear Trend
yu = max(y);
yl = min(y);
yr = (yu-yl); % Range of ‘y’
yz = y-yu+(yr/2);
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0); % Returns Approximate Zero-Crossing Indices Of Argument Vector
zt = x(zci(y));
per = 2*mean(diff(zt)); % Estimate period
ym = mean(y); % Estimate offset
fit = @(b,x) b(1) .* exp(b(2).*x) .* (sin(2*pi*x./b(3) + 2*pi/b(4))) + b(5); % Objective Function to fit
fcn = @(b) norm(fit(b,x) - y); % Least-Squares cost function
[s,nmrs] = fminsearch(fcn, [yr; -10; per; -1; ym]) % Minimise Least-Squares
xp = linspace(min(x),max(x), numel(x)*10);
figure(3)
plot(x,y,'b', 'LineWidth',1.5)
hold on
plot(xp,fit(s,xp), '-r')
hold off
grid
xlabel('Time')
ylabel('Amplitude')
legend('Original Data', 'Fitted Curve')
xlim([min(x) 1])
text(0.3*max(xlim),0.7*min(ylim), sprintf('$y = %.3f\\cdot e^{%.1f\\cdot x}\\cdot sin(2\\pi\\cdot x\\cdot %.0f%.3f)%+.1f$', [s(1:2); 1./s(3:4); s(5)]), 'Interpreter','latex')
format short eng
grid
mdl = fitnlm(x(:),y(:),fit,s) % Get Statistics
.
8 件のコメント
Star Strider
2023 年 2 月 16 日
As always, my pleasure!
The Global Optimization Toolbox has functions such as ga that can make the parameter space search easier and more efficient.
その他の回答 (1 件)
Matt J
2022 年 11 月 28 日
編集済み: Matt J
2022 年 11 月 28 日
Using fminspleas from the File Exchange,
[t,F]= readvars ("https://www.mathworks.com/matlabcentral/answers/uploaded_files/1211598/F300.xlsx");
region=t>=0.2 &t<=0.8;
F=F(region);
t=t(region);
flist={@(p,x)cos(p(1)*t+p(2)).*exp(-t*p(3)),1 };
p0=[0.15,0,1];
[p,AK]=fminspleas(flist,p0,t,F);
f=@(t) AK(1)*cos(p(1)*t+p(2)).*exp(-t*p(3)) +AK(2);
plot(t,F,'o',t,f(t)); legend('Data','Fit')
参考
カテゴリ
Help Center および File Exchange で Surrogate Optimization についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!