% Amplitude spectrum of a sdof system response %emol

clc %emol
clearvars %emol
close all %emol

fs=100; t=[0:1/fs:5-1/fs];
A=200; zeta=0.01; wn=2*pi*20; wd=sqrt(1-zeta^2)*wn;
x=(A/wd)*exp(-zeta*wn*t).*sin(wd*t);

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

H=A./(wn^2-(2*pi*f).^2+i*2*zeta*wn*(2*pi*f)); % true FT %emol

figure(1) %emol
plot(t,x,'b-',t,xh,'r-') %emol
xlabel('t(s)') %emol
ylabel('x(t)') %emol
grid on %emol
legend('rect','hann') %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(sqrt(8/3)*abs(Xh(1:N/2+1)/fs)),'r-')
plot(f(1:N/2+1), 20*log10(abs(H(1:N/2+1))),'g--') %emol
axis([0 50 -150 0])
xlabel('Frequency (Hz)'); ylabel('Modulus (dB)')
legend('rect. window','hann window','true spectrum') %emol
hold off
grid on %emol

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