clc %emol
clearvars %emol
close all %emol

fs=100; % sampling frequency (emol)
t=0:1/fs:2.5-1/fs; % t truncated at 2.5s (emol)

% impulse response of a sdof system (emol)
A=100; zeta=0.03; fn=10; wn=2*pi*fn; wd=sqrt(1-zeta^2)*wn;
h=(A/wd)*exp(-zeta*wn*t).*sin(wd*t); 

% white noise signal as input (emol)
T=100; % 100 and 2000
rng('default') %emol
  % generation of white noise input, band-limited up to fs/2,
  % with variance equal to 4 (emol)
x=2*randn(1,T*fs); 

% output by convolution (emol)
y=conv(h,x); % sequence y not divided by fs (emol)
y=y(1:end-length(h)+1);
% y=filter(h,1,x);

% correlation functions (emol)
maxlag=2.5*fs; % maximum lag at 2.5s (emol)
[Rxx, tau]=xcorr(x,x,maxlag,'unbiased');
[Ryy, tau]=xcorr(y,y,maxlag,'unbiased');
[Rxy, tau]=xcorr(y,x,maxlag,'unbiased');
tau=tau/fs;

% spectral density functions (emol)
% (different scaling factors for relating to the continuous functions, in
%  relative scale, due to the convolution operation in line 20) (emol)
N=2*maxlag;
f=fs*(0:N-1)/N;
Sxx=fft(Rxx(1:N)); 
Syy=fft(Ryy(1:N))/(fs^2);
Sxy=fft(Rxy(1:N))/fs;

% frequency response functions (emol)
H1=Sxy./Sxx; % estimation by spectral densities (emol)
H=fft(h,N)/fs; % Fourier transform of impulse response

figure(1)
plot(tau,Rxx)
xlabel('Lag (\it\tau)'); ylabel('\itR_x_x\rm(\it\tau\rm)')
grid on %emol

figure(2) 
plot(tau,Ryy)
xlabel('Lag (\it\tau)'); ylabel('\itR_y_y\rm(\it\tau\rm)')
grid on %emol

figure(3) 
plot(tau,Rxy)
xlabel('Lag (\it\tau)'); ylabel('\itR_x_y\rm(\it\tau\rm)')
grid on %emol

figure(4)
plot(f(1:N/2+1), 20*log10(abs(H1(1:N/2+1)))); hold on
xlabel('Frequency (Hz)'); ylabel('|\itH\rm_1(\itf\rm)| (dB)')
plot(f(1:N/2+1), 20*log10(abs(H(1:N/2+1))),'r-'); hold off %emol
grid on %emol
legend('|H_1(f)|','|H(f)|'); %emol

figure(5)
plot(f(1:N/2+1), unwrap(angle(H1(1:N/2+1)))); hold on
xlabel('Frequency (Hz)'); ylabel('arg\itH\rm_1(\itf\rm) (rad)')
plot(f(1:N/2+1), unwrap(angle(H(1:N/2+1))),'r-'); hold off %emol
grid on %emol
legend('|H_1(f)|','|H(f)|'); %emol

% Alternative record time: use T = 2000 in line 13 (emol)

% % Comments 1 (still with T = 2000) (minor interest) (emol)
% 
% % The cross-correlation Rxy has a shape similar to the impulse response h
% % (in fact, Rxy=k.h, where k is the variance of the input white noise)
% figure(6) %emol
% plot(t,4*h,'r-'); hold on
% plot(tau(maxlag+1:end), Rxy(maxlag+1:end), 'b-'); hold off %emol
% xlabel('Time (s) and lag (\it\tau\rm)'); ylabel('Amplitude')
% grid on %emol
% legend('k.h(t)','R_x_y(\tau)'); %emol
% 
% % The autocorrelation Ryy is the scaled version of the autocorrelation Rhh
% % (that is, Ryy = k.Rhh)
% figure(7) %emol
% Rhh=xcorr(h,h,maxlag);
% plot(tau,4*Rhh,'r-'); hold on
% plot(tau, Ryy, 'b-'); hold off %emol
% xlabel('Lag (\it\tau\rm)'); ylabel('Amplitude')
% grid on %emol
% legend('k.R_h_h(\tau)','R_y_y(\tau)'); %emol

% % Comments 2 (with T = 100 in line 13) (major interest) (emol)
% 
% % Above, spectral density functions are estimated by the Fourier transform
% % of correlation functions.
% % Below, they are estimated by the segment averaging (Welch's) method.
% Sxx_w=cpsd(x,x, hanning(N),N/2, N, fs, 'twosided');
% Sxy_w=cpsd(y/fs,x, hanning(N),N/2, N, fs, 'twosided');
% H1_w=Sxy_w./Sxx_w;
% 
% figure(6) %emol
% plot(f(1:N/2+1), 20*log10(abs(H1(1:N/2+1))), 'b-'); hold on %emol
% plot(f(1:N/2+1), 20*log10(abs(H1_w(1:N/2+1))),'r-'); hold off %emol
% xlabel('Frequency (Hz)'); ylabel('|\itH\rm_1(\itf\rm)| (dB)')
% grid on %emol
% legend('by correlation functions','by method of Welch') %emol
% 
% figure(7) %emol
% plot(f(1:N/2+1), unwrap(angle(H1(1:N/2+1))), 'b-'); hold on %emol
% plot(f(1:N/2+1), unwrap(angle(H1_w(1:N/2+1))),'r-'); hold off
% xlabel('Frequency (Hz)'); ylabel('arg\itH\rm_1(\itf\rm) (rad)')
% grid on %emol
% legend('by correlation functions','by method of Welch') %emol