clear
close all

SOFAstart;

% Subject index of the file to convert
if ~exist('subject','var'), subject=1; end;
% File name of the RIEC file
RIECfn='./data/RIEC_hrir_subject_';

%SourcePosition
azim=45; %azimuth
elev=0;  %elevation

%% Path definitions
SOFAfile=[RIECfn sprintf('%03d',subject) '.sofa'];

%% Loading the full object
disp(['Loading full object: ' SOFAfile]);
Obj=SOFAload(SOFAfile);

%% Get index of measurements with the same directions
idx=find(Obj.SourcePosition(:,1)==azim & Obj.SourcePosition(:,2)==elev);

idx2=find(Obj.SourcePosition(:,1)==270 & Obj.SourcePosition(:,2)==elev);

%% Extract and plot the fully loaded data
RIEC_hrirL=squeeze(Obj.Data.IR(idx,1,:)); %left ear
RIEC_hrirR=squeeze(Obj.Data.IR(idx,2,:)); %right ear

RIEC_hrtfL=abs(fft(RIEC_hrirL));
RIEC_hrtfR=abs(fft(RIEC_hrirR));

fs=Obj.Data.SamplingRate; %Sampling Frequency
N=length(RIEC_hrirL);
f=fs*(0:N-1)/N;
t=(0:N-1)/fs*1e3;

ss=get(0, 'ScreenSize');
figure('Position', [ss(3)*(1/4) ss(4)*(1/4) ss(3)*0.5 ss(4)*0.5 ])
set(gcf,'defaultAxesFontSize',12);

subplot(2,2,1);
plot(t,RIEC_hrirL);
axis([0 t(end) -1.2 1.2]);
xlabel('Time [ms]');
ylabel('Amplitude');

subplot(2,2,2);
plot(t,RIEC_hrirR);
axis([0 t(end) -1.2 1.2]);
xlabel('Time [ms]');
ylabel('Amplitude');

subplot(2,2,3);
semilogx(f(1:end/2),20*log10(RIEC_hrtfL(1:end/2,1)));
axis([90 fs/2 -60 20]);
xlabel('Frequency [Hz]');
ylabel('Relative level [dB]');

subplot(2,2,4);
semilogx(f(1:end/2),20*log10(RIEC_hrtfR(1:end/2,1)));
axis([90 fs/2 -60 20]);
xlabel('Frequency [Hz]');
ylabel('Relative level [dB]');
