clear all
clc 
close all

% Initialization: let's start by chosing the series of measurements 
% (1.Passive or 2. Active) and save the 8 time signals into a matrix "sig"
[filename,pathname]=uigetfile('Pick the reference signal (velocity) of interest','*.mat');  % Chose the first time record (velocity on channel 1), it will define the 7 other microphone signal
ind1=find(filename=='e');
ind2=find(filename=='.');
Id=eval(filename(ind1+1:ind2-1)); % identify the first index of the time measurement of room transfer function H1
measname=filename(1:ind1)

for ii=Id:Id+7
    eval(['load(''' measname num2str(ii) '.mat'')'])
    eval(['sig(ii-' num2str(Id-1) ',:)=S' measname '_' num2str(ii) ';']);
end

eval(['Fs=S' measname '_' num2str(Id) '_SamplingFreq;' ])
t = [0:size(sig,2)-1]/Fs;
duration=t(end);

% Now we just show a "brute force" estimation of H1: let's apply
% "tfestimate" to the whole time signal (with a frequency resolution of
% 51200 Hz) aiming at a resolution of about 0.1 Hz==> we chose N=2^19 so
% that Fs/N~0.1

N=2^19; % number of points
for ii=2:8
    [H(:,ii-1),f] = tfestimate(sig(1,:),sig(ii,:),kaiser(N,15),[],N,Fs,'Estimator','H1');
    leg(ii-1,:)=['H_1(mic_' num2str(ii-1) '/vel) (Pa.s/m)']     % legend  
end

figure
semilogx(f,20*log10(abs(H)));set(gca,'xlim',[10 150]);
grid
legend(leg)