%% Matthew Neal
% 7/25/2020
% Code examples for JAES Publication: A compact spherical loudspeaker array
% for efficiently recreating instrument directivities.
%
% Plotting array simulation results vs. ideal directivity patterns

clear;
clc;
close all;

addpath(genpath('C:\MATLAB Toolboxes\SOFA')); %SOFA API for MATLAB. See https://www.sofaconventions.org/.
addpath('Matlab Directivity Functions')
addpath('Simulated Directivities - Point Source');
addpath('Simulated Directivities - Spherical Cap');

%% Load ideal instrument directivity from SH Database
addpath('C:\MATLAB Toolboxes\Datasets\Berlin Instrument Directivity'); %add path to TU Berlin Instrument directivity data: http://dx.doi.org/10.14279/depositonce-5861.2

% instrument = 'bassoon';
% instrument = 'cello';
% instrument = 'clarinet';
% instrument = 'double_bass';
% instrument = 'french_horn';
instrument = 'oboe';
% instrument = 'tenor_trombone';
% instrument = 'transverse_flute';
% instrument = 'trumpet';
% instrument = 'tuba';
% instrument = 'viola';
% instrument = 'violin';

era = 'modern';
dynamic = 'ff';
try
    filename = [instrument '_' era '_et_' dynamic '.mat'];
    load(filename);
catch
    filename = [instrument '_et_' dynamic '.mat'];
    load(filename);
end

%% Generate the spatial sampling grid for the simulations
azs = 0:10:360;
els = 90:-5:-90;
[azg, elg] = meshgrid(azs, els);

%% Load the array's radiation control filters for the selected instrument
SOFAstart();

addpath('CSLA Driver Filters');
filt = SOFAload(['ArrayCtrlFilt_' instrument '.sofa']);
fs = filt.Data.SamplingRate;
f = (0:size(filt.Data.IR, 4)-1)*fs/size(filt.Data.IR, 4);
inst_filt = permute(filt.Data.IR(:, :, :, 1:length(f)), [3, 1, 4, 2]);
inst_emitters = filt.EmitterPosition;
clear filt;
inst_filt_tf = fft(inst_filt, [], 3);

%% Filter the point-source simulation with the instrument radiation
%  control filters
driver_sim_tf_PS = zeros(size(azg,1), size(azg,2), size(inst_filt_tf,3));
for nn=1:20
    load(['Simulated_CSLA_pointsource_76mm_baffle_driver_' num2str(nn) '.mat']);
    driver_sim_tf_PS = driver_sim_tf_PS + fft(driver_irs, [], 3).*repmat(inst_filt_tf(nn, 1, :), [size(azg,1), size(azg,2), 1]);
    disp(['Driver ' num2str(nn) ' point-source simulation loaded...']);
end

%% Filter the spherical cap simulation with the instrument radiation
%  control filters
driver_sim_tf_Cap = zeros(size(azg,1), size(azg,2), size(inst_filt_tf,3));
for nn=1:20
    load(['Simulated_CSLA_cap_76mm_baffle_40mm_driver_' num2str(nn) '.mat']);
    driver_sim_tf_Cap = driver_sim_tf_Cap + fft(driver_irs, [], 3).*repmat(inst_filt_tf(nn, 1, :), [size(azg,1), size(azg,2), 1]);
    disp(['Driver ' num2str(nn) ' spherical cap simulation loaded...']);
end

%% Plot comparisons of ideal dirctivity, spherical cap model, and point-
%  source model for the arrays performance for a specific instrument at a 
%  specific frequency
band_index = 18; % change this index to change frequencies!!!

if band_index>=16
    order = 3;
elseif band_index >=14
    order = 2;
elseif band_index>=13
    order = 1;
else
    order = 0;
end

[Ynm,n_ord,m_deg,N_norm] = SphericalHarmonics(azg, 90 - elg, order, 'Format', 'Ortho-Cplx');
bands = radiation.bands.center_frequencies;
Pnm = radiation.pnm(1:(order+1)^2,band_index);

ideal_directivity = db(abs(sum(Ynm.*repmat(permute(Pnm,[2,3,1]),[size(Ynm,1),size(Ynm,2),1]),3)));

% Plot ideal directivity
figure(1);
Balloon3D(ideal_directivity, azg, 90 - elg, 'PlotType', 'color_dB');
title(['Ideal Direcitivty: ' num2str(round(bands(band_index), 0)) ' Hz']);


[~, f_idx] = min(abs(f - bands(band_index)));

% Plot Simulated point-source directivity
plot_data_sim_PS = db(abs(driver_sim_tf_PS(:,:,f_idx)));

figure(2);
Balloon3D(plot_data_sim_PS, azg, 90 - elg, 'PlotType', 'color_dB');
title(['Point-source Simulation: ' num2str(round(bands(band_index), 0)) ' Hz']);

% Plot Simulated spherical cap directivity
plot_data_sim_Cap = db(abs(driver_sim_tf_Cap(:,:,f_idx)));

figure(3);
Balloon3D(plot_data_sim_Cap, azg, 90 - elg, 'PlotType', 'color_dB');
title(['Spherical Cap Simulation: ' num2str(round(bands(band_index), 0)) ' Hz']);