% Amplitude spectrum of a tdof system response %emol

clc %emol
clearvars %emol
close all %emol

fs=100; t=[0:1/fs:5-1/fs];
A=200; B=0.001*A; zeta1=0.01; zeta2=0.01;
wn1=2*pi*20; wd1=sqrt(1-zeta1^2)*wn1;
wn2=2*pi*30; wd2=sqrt(1-zeta2^2)*wn2;
x1=(A/wd1)*exp(-zeta1*wn1*t).*sin(wd1*t); %emol
x2=(B/wd2)*exp(-zeta2*wn2*t).*sin(wd2*t); %emol
x=x1+x2; %emol

N=length(x);
whan=hanning(N); xh=x.*whan';
X=fft(x); Xh=fft(xh);
f=fs*(0:N-1)/N;

H=A./(wn1^2-(2*pi*f).^2+i*2*zeta1*wn1*(2*pi*f))+...
  B./(wn2^2-(2*pi*f).^2+i*2*zeta2*wn2*(2*pi*f)); % true FT %emol

figure(1) %emol
subplot(3,1,1); plot(t,x,'b-'); xlabel('t(s)'); ylabel('x(t)'); %emol
axis([0 5 -2 2]); grid on %emol
subplot(3,1,2); plot(t,x1,'b-'); xlabel('t(s)'); ylabel('x_{1}(t)'); %emol
axis([0 5 -2 2]); grid on %emol
subplot(3,1,3); plot(t,x2,'b-'); xlabel('t(s)'); ylabel('x_{2}(t)'); %emol
axis([0 5 -2 2]); grid on %emol

figure(2) %emol
plot(f(1:N/2+1), 20*log10(abs(X(1:N/2+1)/fs)),'b-'); hold on
plot(f(1:N/2+1), 20*log10(abs(H(1:N/2+1))),'g--') %emol
axis([0 50 -60 0])
xlabel('Frequency (Hz)'); ylabel('Modulus (dB)')
legend('rect. window','true spectrum') %emol
hold off
grid on %emol

figure(3) %emol
plot(f(1:N/2+1), 20*log10(sqrt(8/3)*abs(Xh(1:N/2+1)/fs)),'b-')
hold on
plot(f(1:N/2+1), 20*log10(abs(H(1:N/2+1))),'g--')
axis([0 50 -160 0])
xlabel('Frequency (Hz)'); ylabel('Modulus (dB)')
legend('hann window','true spectrum') %emol
hold off
grid on %emol