Shadron_Pena_2022/HRTF_Wagner_Lab/Scripts/HRIR_Berechnung_translated.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program evaluates the measured data resulted from the
% HRTF-Measurements. The HRIRs associated to all measured 
% positions are determined by correction with the 'Empty-
% Measurement'.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
clc 
close all

tic

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Setup owl and condition (see also annotations.xlsx)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Define Parameters for saving
owl = 'Merry' %#ok<NOPTS>
condition = 'ruff_intact' %#ok<NOPTS>
date_of_experiment = '2003-09-13' %#ok<NOPTS>

% Data directories (in doctorade_datadir) and filename
% Empty-Measurement:
datadir_empty = 'Merry\normal';
filepattern_empty = 'MerrynormalP0P0.dat';
% HRTF-Measurement:
datadir = 'Merry\normal';
filepattern = '*0.dat';

% Measured angles
AZIMUTHS = (-160:20:160); 
ELEVATIONS = (-60:20:80);
% AZIMUTHS = (-160:10:170);
% ELEVATIONS = (-70:10:80);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Load doctorade_config, e.g. doctorade_datadir
doctorade_datadir = 'D:\MATLAB\HRTF_Wagner_Lab\raw_data';
doctorade_resultdir = 'D:\MATLAB\HRTF_Wagner_Lab\reults2';

% Samplingrate
samplingrate = 100000; % in Hz

% Set parameters ("arbitrary" but fixed)
window_length = 3.6; % length of the part of the HRIRs that is non-zero [in ms]
hrir_length = 500; % corresponds to 5 ms
begin = 400; % lower cutoff frequency of the band-pass filter [in Hz]
ende = 9000;  % upper cutoff frequency of the band-pass filter [in Hz]
bandwidth = [num2str(begin) ' to ' num2str(ende) ' Hz'] %#ok<NOPTS>
filterordnung = 10; % half of the filter order


L_lang = 2 * 51000; % length of the FFT (twice as long as the measuremant)
frequencies = (0:L_lang-1)*samplingrate/L_lang;
    
% Matrices for saving the HRIRs 
hrirshifted_l = nan(numel(AZIMUTHS),numel(ELEVATIONS),L_lang);
hrirshifted_r = nan(numel(AZIMUTHS),numel(ELEVATIONS),L_lang);

    
% Highpass filter (top-hat filter with cutoff frequency 50Hz)
highpass = zeros(1,L_lang/2);
begin_high = 50;
highpass(1,ceil(begin_high/100000*L_lang)+1:L_lang/2) = 1;
highpass = [highpass 1 fliplr(highpass(1,2:length(highpass)))];

     
% Band-pass filter (Butterworth filter with cutoff frequencies 400 and 9000 Hz and order 20)
[z,p,k] = butter(filterordnung ,[begin/(samplingrate/2) ende/(samplingrate/2)], 'bandpass'); 
sos = zp2sos(z,p,k);
[h,f] = freqz(sos,L_lang/2,samplingrate);
filter = [h.' 0 conj(fliplr(h(2:length(h)).'))];

%% Empty-Measurement  
% The Empty-Measurement at (0°,0°) is used to determine the transfer functions of the loudspeakers, microphones, etc. These are needed later for correction.
namesempty = dir(fullfile(doctorade_datadir, datadir_empty, filepattern_empty)); % lists all Empty-Measurements at (0°,0°)
leftmeanallempty = zeros(length(namesempty),L_lang);
rightmeanallempty = zeros(length(namesempty),L_lang);
for i = 1:1:(length(namesempty)) % averages across five sweeps (every single measurement consits of five sweeps)
    filenameempty = namesempty(i).name;

    [rightempty , leftempty] = hrtf_dateilesen(fullfile(doctorade_datadir, datadir_empty, filenameempty)); 

    leftmeanempty = mean(leftempty(1:5,:));
    rightmeanempty = mean(rightempty(1:5,:));
    leftmeanempty = ifft(fft(leftmeanempty, L_lang) .* highpass);
    rightmeanempty = ifft(fft(rightmeanempty, L_lang) .* highpass);
    leftmeanallempty(i,:) = leftmeanempty;
    rightmeanallempty(i,:) = rightmeanempty;
end

leftmeanempty = mean(leftmeanallempty(1:i,:),1);% averages across five different measurements
rightmeanempty = mean(rightmeanallempty(1:i,:),1); 

time = (0:L_lang-1)*1/samplingrate;

fftleftempty = fft(leftmeanempty, L_lang);% is needed for corretion
fftrightempty = fft(rightmeanempty, L_lang); 

%% Owl-HRIRs gets determined and transformed into matrix-form
names = dir(fullfile(doctorade_datadir, datadir, filepattern)); % lists all owl-data
for i = 1:1:(length(names)) 
    filename = names(i).name;
    disp(filename); 

    [azi,ele] = filename2azel(filename); % determines azimuth and elevation of each measurement
    row = find(AZIMUTHS == azi);
    column = find(ELEVATIONS == ele);

    [right , left] = hrtf_dateilesen(fullfile(doctorade_datadir, datadir, filename));

    rightmean = mean(right(1:5,:)); % averages across five sweeps (every single measurement consits of five sweeps)
    leftmean = mean(left(1:5,:));
    rightmean = ifft(fft(rightmean, L_lang) .* highpass);
    leftmean = ifft(fft(leftmean, L_lang) .* highpass);

    fhr = fft(rightmean, L_lang);
    fhl = fft(leftmean, L_lang);

    fftIRright = (fhr./fftrightempty);% correction via division by Empty-Measurement in the frequency domain
    fftIRleft = (fhl./fftleftempty);

    fftIRleftband = fftIRleft .* filter; 
    fftIRrightband = fftIRright .* filter; 

    IRleftbandpass = ifft(fftIRleftband, 'symmetric');
    IRrightbandpass = ifft(fftIRrightband, 'symmetric');

    % Save HRIRs
    hrirshifted_l(row,column,:) = IRleftbandpass;
    hrirshifted_r(row,column,:) = IRrightbandpass;
end 

%% Further transformations of the HRIR-Matrices
% Correct the HRIR's offset (cyclic shift)
[hrir_l, hrir_r, az_num, el_num] = offset(hrirshifted_l, hrirshifted_r, window_length, 1); % third argument in ms; fourth argument switches smoothing on or off 

% Shorten the HRIRs to a length of 5 ms  (or p.r.n. a little bit more)
limit = 1; % to make sure that the 5 ms window does not cut away information
for i = 1:numel(AZIMUTHS)
    for j = 1:numel(ELEVATIONS)
        signal_l = squeeze(hrir_l(i,j,:)).';
        signal_r = squeeze(hrir_r(i,j,:)).';
        pos_l = find(signal_l ~= 0);
        pos_r = find(signal_r ~= 0);
        end_l = max(pos_l);
        end_r = max(pos_r);
        end_gesamt = max(end_l, end_r);
        if end_gesamt > limit
            limit = end_gesamt;
        end
    end
end
    
hrirshort_l = hrir_l(:, :, 1:max(hrir_length,limit)); 
hrirshort_r = hrir_r(:, :, 1:max(hrir_length,limit));

signallength_short = size(hrirshort_l,3);
frequencies_short = (0:(signallength_short-1))*samplingrate/signallength_short; % frequencies from 0 to ~ 100000 Hz with a new step size
time_short = (0:(signallength_short-1))*1/samplingrate;

% Normalization of the HRIRs (for the EPhys-Experiments)
idx_azi_0 = find(AZIMUTHS == 0);
idx_ele_0 = find(ELEVATIONS == 0);
left = squeeze(hrirshort_l(idx_azi_0,idx_ele_0,:)).';
right = squeeze(hrirshort_r(idx_azi_0,idx_ele_0,:)).';
energy_l = left * left.';
energy_r = right * right.';
energy = sqrt((energy_l + energy_r)/2);
hrirshort_l = 1/energy * hrirshort_l;
hrirshort_r = 1/energy * hrirshort_r;

%% Create and save matlabfile for HRIR use in electrophysiological experiments
hrir_l = hrirshort_l; % for saving the short HRIRs are renamed to hrir_l, hrir_r 
hrir_r = hrirshort_r;
azimuths = AZIMUTHS';
elevations = ELEVATIONS;
dataname = [owl '__' condition '__hrir_result'];
name = ['Owl: ' owl ';' ' Condition: ' condition ';' ' Date of experiment: ' num2str(date_of_experiment) '; ' 'Bandwidth: ' num2str(bandwidth)];
save(fullfile(doctorade_resultdir, dataname), 'hrir_l' , 'hrir_r' , 'azimuths' , 'elevations' , 'samplingrate' , 'name');

toc