    %Measuring directional characteristics of sound sources
%Made by Benjamin Nimmerrichter 2022/23 for BUT

%analyze measured data from INGEST folder
% make sure to fill out the mic_pos file
clear
%% SETTINGS
option = 1; %1 = polar plot, 2 = spherical plot, 3 = fractional octave char - single band, 
            %4 = fractional octave single mic - multiband, 5 = single mic FFT
mic_ID = 1;                 % for single mic, otherwise ignored
band_ID = 10;                % for multi mic songle band, otherwise ignored
octave_frac = '1/3 octave'; % for setting the octave fractions try '1/3 octave', when empty it defaults to octave;
file_ID = 1;                % if you want to analyze a single file, set to 0 if you want to analyze all files
SPL = true;                 % if true, measurements will load reference

%% LOAD FILES
if isfile("INGEST\mic-pos.txt")
mic_pos = readmatrix("INGEST\mic-pos.txt");
disp('Soubor umístění mikrofonů byl načten.')


mics = max(mic_pos(:,1));
passes = max(mic_pos(:,2));


if ~(option == 2 || option == 1)
    passes = 1; %prevent problems with non compatible graphs
end
files = dir('INGEST\*.mat');
fix = 0;
for i=1:length(files)
    if (files(i).name == "config.mat")
    load("INGEST\config.mat",'conf');
    fix = -1;
    else
    load(['INGEST\' files(i).name]);
    
    for j=1:mics*passes
    data_dB(i+fix,j) = meas.dB(j);
    data_raw(i+fix,j,:) = zeros(length(meas.raw(:,j)),1);
    data_raw(i+fix,j,:) = meas.raw(:,j);
    end
    end
end
if ~exist("conf","var") %if you don't input a config file this prevents errors
    SPL = false;
    Fs = 48000;
else
    Fs = conf.Fs;
end

%% DATA PROCESSING

%INTERPOLATION
if (mics*passes > 1 )

    
azi = mic_pos(:,3); %azi
elev = mic_pos(:,4); %elev
if (file_ID == 0)
spl = data_dB(:,:); %spl
file_ID = 1;
else
spl = data_dB(file_ID,:); %spl   
end

if (option == 4 && option == 5) % for single microphone measurements
data_dB = data_dB(file_ID,mic_ID);
data_raw = data_raw(file_ID,mic_ID,:);
end

for i = 1:passes
start = mics*(i-1)+1;
stop = mics*i;
if (i == 1)
Iazi = interpn(start:stop,azi(start:stop),start:0.1:stop,'cubic');
Ielev = interpn(start:stop,elev(start:stop),start:0.1:stop,'cubic');
Ispl = round(interpn(start:stop,spl(start:stop),start:0.1:stop,'cubic'),1);
else
Iazi = cat(2,Iazi,interpn(start:stop,azi(start:stop),start:0.1:stop,'cubic'));
Ielev = cat(2,Ielev,interpn(start:stop,elev(start:stop),start:0.1:stop,'cubic'));
Ispl = cat(2,Ispl,round(interpn(start:stop,spl(start:stop),start:0.1:stop,'cubic'),1));
end
end



%% SPHERICAL AND POLAR PLOT

offset = min(Ispl);
divider = max(Ispl)-offset;
minSpl = min(Ispl);
maxSpl = max(Ispl);
for k = 1:size(Iazi,2)
in(1,k) = Ielev(k);
in(2,k) = Iazi(k);
in(3,k) = ((Ispl(k)-(offset))/divider);
end
end
clear i j k;
if (option == 1)
    polarp;
elseif(option == 2)
    spherical;
elseif(option == 3)
    sz = mics*passes;
    for i = 1:size(data_raw,1)
        for j = 1:size(data_raw,2)
            signal = reshape(data_raw(i,j,:),[],1);
            
            if (SPL == true)
                [oct_out(i,j,:),cf_out(i,j,:)] = frac_oct(signal,octave_frac,Fs,false,conf.calibration(j));
            else 
                [oct_out(i,j,:),cf_out(i,j,:)] = frac_oct(signal,octave_frac,Fs,false);
            end
        end
    end
    IdB = round(interpn(1:sz,oct_out(1,1:sz,band_ID),1:0.1:sz,'cubic'),1);
    dB_off = min(IdB);
    dB_div = max(IdB)-dB_off;
    maxSpl = max(IdB);
    for k = 1:size(Iazi,2)
        in(3,k) = ((IdB(k)-(dB_off))/dB_div);
    end
    polarp;
elseif(option == 4)
    signal = reshape(data_raw(file_ID,mic_ID,:),[],1);
    if (SPL == true)
        [oct_out,cf_out] = frac_oct(signal,octave_frac,Fs,true,conf.calibration(mic_ID));
    else 
        [oct_out,cf_out] = frac_oct(signal,octave_frac,Fs,true);
    end
elseif(option == 5)
    signal = reshape(data_raw(file_ID,mic_ID,:),[],1);
    if (SPL == true)
        dft(signal,Fs,conf.calibration(mic_ID))
    else
        dft(signal,Fs)
    end
end

else
disp('Nebylo možné načíst soubor umístění mikrofonů.');
end