% *******************************************************************
% This script plots the acoustic transfer functions of the resonators
% made of silicone and PLA and determines the frequencies and 
% bandwidths of their first three resonances.
%
% Copyright (C) by Peter Birkholz, TU Dresden, Germany, 2021
% *******************************************************************

resonatorNames = {'a', 'e', 'i', 'o', 'u', 'ae', 'oe', 'y', 'schwa', 'cylinder'};
numResonators = length(resonatorNames);

f1Pla_Hz = zeros(numResonators, 1);
f2Pla_Hz = zeros(numResonators, 1);
f3Pla_Hz = zeros(numResonators, 1);
bw1Pla_Hz = zeros(numResonators, 1);
bw2Pla_Hz = zeros(numResonators, 1);
bw3Pla_Hz = zeros(numResonators, 1);
f1Silicone_Hz = zeros(numResonators, 1);
f2Silicone_Hz = zeros(numResonators, 1);
f3Silicone_Hz = zeros(numResonators, 1);
bw1Silicone_Hz = zeros(numResonators, 1);
bw2Silicone_Hz = zeros(numResonators, 1);
bw3Silicone_Hz = zeros(numResonators, 1);

disp(' ');

for n=1:1:numResonators
    [f1Pla_Hz(n), bw1Pla_Hz(n), f1Silicone_Hz(n), bw1Silicone_Hz(n), ...
     f2Pla_Hz(n), bw2Pla_Hz(n), f2Silicone_Hz(n), bw2Silicone_Hz(n), ...
     f3Pla_Hz(n), bw3Pla_Hz(n), f3Silicone_Hz(n), bw3Silicone_Hz(n)] = ...
        plotResonatorCurves(n, resonatorNames{n});    
    
    % Print the resonance frequencies and bandwidths.
    disp([resonatorNames{n} ':  ' num2str(round(f1Pla_Hz(n))), ' (', ...
         num2str(round(bw1Pla_Hz(n))), ')   ' ...
         num2str(round(f1Silicone_Hz(n))), ' (', ...
         num2str(round(bw1Silicone_Hz(n))), ')   ', ...
         num2str(round(f1Silicone_Hz(n) - f1Pla_Hz(n)))]);
end

disp(' ');

% Plot the shift of the resonances due to the soft walls.

allFreqPla = [ f1Pla_Hz; f2Pla_Hz; f3Pla_Hz ];
allBwPla = [ bw1Pla_Hz; bw2Pla_Hz; bw3Pla_Hz ];

allFreqSilicone = [ f1Silicone_Hz; f2Silicone_Hz; f3Silicone_Hz ];
allBwSilicone = [ bw1Silicone_Hz; bw2Silicone_Hz; bw3Silicone_Hz ];

subplot(4, 3, 11);
plot(allFreqPla, allFreqSilicone - allFreqPla, 'ko');
xlabel('f_{R1,hard} in Hz');
ylabel('f_{R1,soft} - f_{R1,hard} in Hz');
grid on;

% Plot the bandwidths of the resonances for the soft and the hard models
% as a function of the frequency.

subplot(4, 3, 12);
plot(allFreqPla, allBwPla, 'ko', allFreqSilicone, allBwSilicone, 'ro');
xlabel('f_{R1} in Hz');
ylabel('Bandwidth in Hz');
grid on;


% *******************************************************************
% *******************************************************************

function [f1Pla_Hz, bw1Pla_Hz, f1Silicone_Hz, bw1Silicone_Hz, ...
          f2Pla_Hz, bw2Pla_Hz, f2Silicone_Hz, bw2Silicone_Hz, ...
          f3Pla_Hz, bw3Pla_Hz, f3Silicone_Hz, bw3Silicone_Hz] = ...
    plotResonatorCurves(index, resonatorName)

    subplot(4, 3, index);

    fileNamePla = [resonatorName '-pla.txt'];
    fileNameSilicone = [resonatorName '-silicone.txt'];

    if (~exist(fileNamePla, 'file')) || (~exist(fileNameSilicone, 'file'))
        plot([0, 1], [0, 1], 'r');
        return;
    end

    % This offset has to be subtracted from the magnitude spectra to
    % compensate for the differences in the microphone sensitivities.
    offset_dB = 13.5;
    
    data = dlmread(fileNamePla, ' ', 1, 0);    % Skip the first row.
    frequencies_Hz = data(:, 1);
    magnitudesPla_dB = 20.0*log10(data(:, 2)) - offset_dB;

    data = dlmread(fileNameSilicone, ' ', 1, 0);    % Skip the first row.
    frequencies_Hz = data(:, 1);
    magnitudesSilicone_dB = 20.0*log10(data(:, 2)) - offset_dB;

    plot(frequencies_Hz, magnitudesPla_dB, 'k', ...
         frequencies_Hz, magnitudesSilicone_dB, 'r');

    xlim([0, 5000]);
    ylim([-34, 37]);
    title(['Resonator ' resonatorName]);
    xlabel('Frequency in Hz');
    ylabel('Magnitude in dB');
    grid on;
    
    % Rough visual estimates of the resonance frequencies as initial values 
    % for the fine measurements. Values of zeros are resonances to be ignored.
    
    freqEstimatePla = [ ...
        603, 1150, 2320;
        297, 2055, 2640;
        190, 1925, 2605;
        295, 598, 0;
        222, 740, 0;
        473, 1790, 2415;
        280, 1630, 1903;
        237, 1645, 1980;
        406, 1520, 2211;
        485, 1470, 2444 ];

    freqEstimateSilicone = [ ...
        713, 1150, 2320;
        437, 2055, 2640;
        327, 1925, 2605;
        440, 598, 0;
        390, 740, 0;
        590, 1790, 2415;
        426, 1630, 1903;
        393, 1645, 1980;
        505, 1520, 2211;
        605, 1470, 2444 ];
    
    [f1Pla_Hz, bw1Pla_Hz] = getResonanceFrequencyAndBandwidth(frequencies_Hz, magnitudesPla_dB, freqEstimatePla(index, 1));
    [f2Pla_Hz, bw2Pla_Hz] = getResonanceFrequencyAndBandwidth(frequencies_Hz, magnitudesPla_dB, freqEstimatePla(index, 2));
    [f3Pla_Hz, bw3Pla_Hz] = getResonanceFrequencyAndBandwidth(frequencies_Hz, magnitudesPla_dB, freqEstimatePla(index, 3));

    [f1Silicone_Hz, bw1Silicone_Hz] = getResonanceFrequencyAndBandwidth(frequencies_Hz, magnitudesSilicone_dB, freqEstimateSilicone(index, 1));
    [f2Silicone_Hz, bw2Silicone_Hz] = getResonanceFrequencyAndBandwidth(frequencies_Hz, magnitudesSilicone_dB, freqEstimateSilicone(index, 2));
    [f3Silicone_Hz, bw3Silicone_Hz] = getResonanceFrequencyAndBandwidth(frequencies_Hz, magnitudesSilicone_dB, freqEstimateSilicone(index, 3));
    
    hold on

    minY = -34;
    maxY = 37;
    
    plot([f1Pla_Hz, f1Pla_Hz], [minY, maxY], 'k', ...
         [f2Pla_Hz, f2Pla_Hz], [minY, maxY], 'k', ...
         [f3Pla_Hz, f3Pla_Hz], [minY, maxY], 'k', ...
         [f1Silicone_Hz, f1Silicone_Hz], [minY, maxY], 'r', ...
         [f2Silicone_Hz, f2Silicone_Hz], [minY, maxY], 'r', ...
         [f3Silicone_Hz, f3Silicone_Hz], [minY, maxY], 'r' ...
    );
    hold off;

end


% *******************************************************************
% Determine the resonance frequency and bandwidth around the given
% initialFreqEstimate_Hz in the given spectrum.
% *******************************************************************

function [freq_Hz, bw_Hz] = getResonanceFrequencyAndBandwidth(frequencies_Hz, magnitudes_dB, initialFreqEstimate_Hz)
    
    if (initialFreqEstimate_Hz == 0)
        freq_Hz = 0;
        bw_Hz = 0;
        return;
    end
    
    % Get the indices of the search range.
    deltaFreq_Hz = 100;  % Search between 100 Hz left and right from the initial estimate.
    minSearchIndex = interp1(frequencies_Hz, 1:length(frequencies_Hz), initialFreqEstimate_Hz - deltaFreq_Hz, 'nearest');
    maxSearchIndex = interp1(frequencies_Hz, 1:length(frequencies_Hz), initialFreqEstimate_Hz + deltaFreq_Hz, 'nearest');
    
    [dummy, peakIndex] = max(magnitudes_dB(minSearchIndex:maxSearchIndex));
    
    peakIndex = peakIndex + minSearchIndex - 1;
    
    freq_Hz = frequencies_Hz(peakIndex);
    peakAmp_dB = magnitudes_dB(peakIndex);
    
    % Go left and right from the peak frequency until the magnitude dropped
    % by 3 dB.
    index = peakIndex;
    while ((index > 1) && (magnitudes_dB(index) > peakAmp_dB - 3.0))
        index = index - 1;
    end
    leftFreq_Hz = frequencies_Hz(index);
    
    index = peakIndex;
    while ((index < length(frequencies_Hz)) && (magnitudes_dB(index) > peakAmp_dB - 3.0))
        index = index + 1;
    end
    rightFreq_Hz = frequencies_Hz(index);
    
    bw_Hz = rightFreq_Hz - leftFreq_Hz;
end


   