%% Matthew Neal
% 7/25/2020
% Code examples for JAES Publication: A compact spherical loudspeaker array
% for efficiently recreating instrument directivities.
%
% Plotting turntable measurement results

clear;
clc;
close all;

addpath(genpath('C:\MATLAB Toolboxes\SOFA')); %SOFA API for MATLAB. See https://www.sofaconventions.org/.
addpath('Matlab Directivity Functions')

%% 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

%% Load Measured Directivity
addpath('Measured Directivities - Turntable');
SOFAstart();

meas = SOFAload(['Measured_SPRAL_CSLA_TFs_' instrument '.sofa']);

[meas_data, azg, elg] = Sofa2GridTF(meas, 10, 5); % convert to grid-like data for visualization

%% Calculate ideal directivity pattern from spherical harmonic weights
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)));

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

% Plot turntable measurements at same frequency (bin)
[~, f_idx] = min(abs(meas.N - bands(band_index)));

figure(2);
Balloon3D(db(abs(meas_data(:,:,f_idx))), azg, 90 - elg, 'PlotType', 'color_dB');
title(['Measured Directivity: ' num2str(round(bands(band_index), 0)) ' Hz']);