% ************************************************************************
% Load the precise transfer functions (radiation characteristics) with a 
% high frequency resolution (25 Hz) for the simple radiation model, for 
% the piston-in-sphere model, and for the head-torso model, and plot their 
% magnitude responses.
% ************************************************************************

T = readtable('radiation_characteristics.txt');

frequencies = T{:, 1};
H_simple_dB = 20*log10(T{:, 2});
H_piston_sphere_dB = 20*log10(T{:, 4});
H_head_torso_dB = 20*log10(T{:, 6});

semilogx(...
    frequencies,    H_simple_dB, 'k', ...
    frequencies,    H_piston_sphere_dB, 'b', ...
    frequencies,    H_head_torso_dB, 'r');

title('Magnitude of the radiation characteristics: simple (black), piston in sphere (blue), head-torso (red)');
xlabel('Frequency in Hz');
ylabel('Magnitude in dB');
xlim([100, 24000]);
grid on;

% ************************************************************************
% Plot the frequency response of one of the FIR filter approximations.
% ************************************************************************

figure;
coeff = readmatrix('FIR_filter_head_torso_16000Hz.txt');
freqz(coeff, [1]);
% Make the frequency axis logarithmic for a better comparison with the
% previous plot.
ax = findall(gcf, 'Type', 'axes');
set(ax, 'XScale', 'log');


