clear all;close all;
syms A dr w1 r  t;
A=0.02;w1=3;r=0;dr=0.06;
x=A.*exp(-dr*w1.*t).*sin(sqrt(1-dr^2)*w1.*t+r);
t=linspace(1,20/w1,5120);
subx=subs(x);
figure('Name','theoretical free decline');
subplot(211);
plot(t,subs(x),'LineWidth',2);
for i=2:(length(t)-1)
    if (((subx(i)-subx(i-1))/(t(i)-t(i-1)))*((subx(i+1)-subx(i))/(t(i+1)-t(i))))<0
        text(t(i),subx(i),[num2str(t(i)),'s,',num2str(subx(i))], 'HorizontalAlignment','center','Edgecolor','blue', 'FontSize',8, 'LineStyle',':');
    else
    end
end
hold on;
grid on;
for n=1:length(t)
x2(n)=sum(subx(1:n).*(t(2)-t(1)));
end
x22=x2-mean(x2);
plot(t,x2,'c');plot(t,x22,'m','LineWidth',1);xlabel('t/s');ylabel('magnitude/m');title(['theoretical free decline OmegaD:',num2str(sqrt(1-dr^2)*w1),' OmegaN:',num2str(w1),' epsilon:',num2str(dr)]);
axis([1 20/w1 -0.02 0.03])
for i=2:(length(t)-1)
    if (((x22(i)-x22(i-1))/(t(i)-t(i-1)))*((x22(i+1)-x22(i))/(t(i+1)-t(i))))<0
        text(t(i),x22(i),[num2str(t(i)),'s,',num2str(x22(i))], 'HorizontalAlignment','center','Edgecolor','m', 'FontSize',6, 'LineStyle',':');
    else
    end
end
legend('acceleration','integral velocity','zeroed velocity','location','best','Orientation','horizontal');
Y1=fft(subx);
Ts1=t(2)-t(1);
Ws1=2*pi/Ts1;
Wn1=Ws1/2;
w1=linspace(0,Wn1,length(t)/2);
Ya1=abs(Y1(1:length(t)/2));
i1=find(Ya1>=mean(Ya1));
ii1=find(Ya1==max(Ya1));
subplot(212);
plot(w1(1:i1(length(i1))),Ya1(1:i1(length(i1))),'LineWidth',2);xlabel('frequncy/Hz');ylabel('manitude');title('focused frequncy spectrum >mean value');
text(w1(ii1),Ya1(ii1),[num2str(w1(ii1)),'Hz,',num2str(Ya1(ii1))], 'HorizontalAlignment','center','Edgecolor','blue', 'FontSize',8, 'LineStyle',':');
text(w1(length(i1)),mean(Ya1)*3,[num2str(mean(Ya1))], 'HorizontalAlignment','center','Edgecolor','blue', 'FontSize',8, 'LineStyle',':');
grid on;hold on;
Y2=fft(x22);
Ts2=t(2)-t(1);
Ws2=2*pi/Ts2;
Wn2=Ws2/2;
w2=linspace(0,Wn2,length(t)/2);
Ya2=abs(Y2(1:length(t)/2));
i2=find(Ya2>=mean(Ya2));
ii2=find(Ya2==max(Ya2));
subplot(212);
plot(w2(1:i2(length(i2))),Ya2(1:i2(length(i2))),'m');axis([0 50 0 30])
text(w2(ii2),Ya2(ii2),[num2str(w2(ii2)),'Hz,',num2str(Ya2(ii2))], 'HorizontalAlignment','center', 'FontSize',8,'Edgecolor','m', 'LineStyle',':');
text(w2(length(i2)),mean(Ya2)*3,[num2str(mean(Ya2))], 'HorizontalAlignment','center','Edgecolor','m', 'FontSize',8, 'LineStyle',':');