clc %emol
clearvars %emol
close all %emol

load beam_experiment % signals recorded for 20s (emol)

fs=256; % sampling rate (emol)
T=4; % length of time segments (original value: 4) (emol)
TT=20; % length of total time (emol)

% spectral density functions:            (emol)
% Hann window, 50% overlap (9 averages), (emol)
% no zero padding, physical frequency.   (emol)
[Gxx,f]=cpsd(x,x,hanning(T*fs),T*fs/2,T*fs,fs);  
[Gyy,f]=cpsd(y,y,hanning(T*fs),T*fs/2,T*fs,fs);
[Gxy,f]=cpsd(y,x,hanning(T*fs),T*fs/2,T*fs,fs);
[Gyx,f]=cpsd(x,y,hanning(T*fs),T*fs/2,T*fs,fs);

% FRF estimators %emol
H1=Gxy./Gxx; 
H2=Gyy./Gyx;
HT=(Gyy-Gxx + sqrt((Gxx-Gyy).^2 + 4*abs(Gxy).^2))./(2*Gyx);

% 1st stage: time signals (emol)
figure(1) %emol
t=0:1/fs:TT-(1/fs); %emol
plot(t,x,'r-',t,y,'b-') %emol
xlabel('Time (s)'); %emol
ylabel('x(t), y(t)') %emol
legend('x(t)','y(t)') %emol
grid on %emol

% 2nd stage: magnitude spectra of estimators (emol)
figure (2)
plot(f,20*log10(abs(H1)),'r-')
xlabel('Frequency (Hz)'); 
ylabel('|\itH\rm_1(\itf\rm)| (dB)')
axis([5 90 -45 25])
grid on %emol

figure (3)
plot(f,20*log10(abs(H2)),'g-')
xlabel('Frequency (Hz)'); 
ylabel('|\itH\rm_2(\itf\rm)| (dB)')
axis([5 90 -45 25])
grid on %emol

figure (4)
plot(f,20*log10(abs(HT)),'b-')
xlabel('Frequency (Hz)'); 
ylabel('|\itH_T(\itf\rm)| (dB)')
axis([5 90 -45 25])
grid on %emol

% 3rd stage: magnitude spectra zoomed at resonances (emol)
figure(5) %emol
subplot(1,2,1) %emol
plot(f(37:54),20*log10(abs(H1(37:54))),'r-',...
     f(37:54),20*log10(abs(H2(37:54))),'g-',...
     f(37:54),20*log10(abs(HT(37:54))),'b-') %emol
grid on %emol
xlabel('Frequency (Hz)'); ylabel('Magnitude spectrum (dB)') %emol
title('First resonance') %emol
legend('H1','H2','HT') %emol
subplot(1,2,2) %emol
plot(f(224:277),20*log10(abs(H1(224:277))),'r-',...
     f(224:277),20*log10(abs(H2(224:277))),'g-',...
     f(224:277),20*log10(abs(HT(224:277))),'b-') %emol
grid on %emol
xlabel('Frequency (Hz)'); ylabel('Magnitude spectrum (dB)') %emol
title('Second resonance') %emol
legend('H1','H2','HT') %emol

% 4th stage: phase spectra of estimators (emol)
figure(6) %emol
plot(f(21:361),unwrap(angle(H1(21:361))),'r-', ...
     f(21:361),unwrap(angle(H2(21:361))),'g-', ...
     f(21:361),unwrap(angle(HT(21:361))),'b-');
xlabel('Frequency (Hz)'); ylabel('Phase spectrum (rad)')
axis([5 90 -7 0.5])
grid on %emol
legend('H1','H2','HT') %emol

% 5th stage: coherence function (emol)
% (drops at resonances could be reduced by improving frequency resolution;
%  for instance, change T=4 by T=10 in line 8, to see just that)
figure(7) %emol
Gamma=abs(Gxy).^2./(Gxx.*Gyy);
plot(f,Gamma,'r-')
axis([5 90 0 1.1])
xlabel('Frequency (Hz)'); ylabel('Coherence function')
grid on %emol
