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

resonatorNames = {'tense-a-pla', 'tense-e-pla', 'tense-i-pla', ...
    'tense-o-pla', 'tense-u-pla', 'tense-ae-pla', 'tense-oe-pla', ...
    'tense-y-pla', 'lax-ae-pla', 'lax-i-pla', 'lax-o-pla', ...
    'lax-u-pla', 'lax-y-pla', 'lax-oe-pla', 'schwa-pla', 'tiefschwa-pla'};

numResonators = length(resonatorNames);

% *******************************************************************
% Initial estimates of the first 6 resonance frequencies to help the 
% peak-picking algorithm. Values of zero mean "don't use", for example,
% because they are above 6 kHz.
% *******************************************************************

freqEstimate = [ ...
    598, 1151, 2326, 2861, 4250, 5896; ...       % /tense-a/
    296, 2064, 2630, 3345, 4682, 0; ...       % /tense-e/
    189, 1914, 2617, 3287, 4300, 5241; ...       % /tense-i/
    293, 594, 3213, 4132, 5090, 6087; ...       % /tense-o/
    226, 743, 2823, 3997, 4498, 6196; ...       % /tense-u/
    472, 1796, 2414, 3055, 4532, 5558; ...       % /tense-ae/
    276, 1638, 1892, 2895, 3807, 4662; ...       % /tense-oe/
    236, 1642, 1983, 2847, 4226, 4883; ...       % /tense-y/
    475, 1672, 2337, 2978, 4529, 5548; ...       % /lax-ae/
    307, 1767, 2346, 2895, 4428, 5701; ...       % /lax-i/
    471, 895, 2281, 3069, 4072, 5189; ...       % /lax-o/
    319, 791, 0, 3104, 3817, 4836; ...          % /lax-u/
    354, 1303, 2109, 2662, 3907, 4744; ...       % /lax-y/
    458, 1352, 2128, 2832, 4114, 5145; ...       % /lax-oe/
    402, 1518, 2211, 2782, 4431, 5462; ...      % /schwa/
    575, 1332, 2209, 2908, 4274, 5465];      
    

% *******************************************************************
% Final lists of determined frequencies and bandwidths.
% *******************************************************************

vtlResonatorsFreq_Hz = [];
vtlResonatorsBw_Hz = [];
vtlResonatorsIndex = [];        % Resonator index: 1, 2, 3, 4, ...

disp(' ');

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

for n = 1:1:numResonators
    
    if ((n == 1) || (n == 9))
        figure;
    end
    subplot(3, 3, 1 + mod(n-1, 8));
    
    % The first index in the function calls is the subplot index of a 3 x 3 
    % plot matrix.
    fileName = [resonatorNames{n} '.txt'];
    plotSpectrum(fileName);

    [peakFreq, bwLeftFreq, bwRightFreq] = getResonanceFreqAndBw([resonatorNames{n} '.txt'], freqEstimate(n, :));
    plotFoundResonances(peakFreq, bwLeftFreq, bwRightFreq);
    vtlResonatorsFreq_Hz = [vtlResonatorsFreq_Hz peakFreq];
    vtlResonatorsBw_Hz = [vtlResonatorsBw_Hz bwRightFreq - bwLeftFreq];
    vtlResonatorsIndex = [vtlResonatorsIndex n*ones(1, length(peakFreq))];
end

% Final plot with the bandwidths...

figure;
plot(vtlResonatorsFreq_Hz, vtlResonatorsBw_Hz, 'ko','MarkerFaceColor','black');
grid on;
xlabel('Frequency in Hz');
ylabel('Bandwidth in Hz');

% Export the frequency and bandwidths data into a CSV file.

titles = {'resonator_index', 'frequency_Hz','bandwidth_Hz'};
data = [vtlResonatorsIndex', vtlResonatorsFreq_Hz', vtlResonatorsBw_Hz']
writecell([titles; num2cell(data)], '../tubes_frequencies_bandwidths.txt')

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

function plotSpectrum(primaryFileName)

    % Plot the primary spectrum.
    if (exist(primaryFileName, 'file'))
        data = dlmread(primaryFileName, ' ', 1, 0);    % Skip the first row.
        frequencies_Hz = data(:, 1);
        magnitudes_dB = 20.0*log10(data(:, 2));
        plot(frequencies_Hz, magnitudes_dB, 'k');
    end

    xlim([0, 6500]);
    ylim([-20, 50]);
    title(primaryFileName);
    xlabel('Frequency in Hz');
    ylabel('Magnitude in dB');
    grid on;
end


% *******************************************************************
% Mark the resonances and bandwiths in the given subplot.
% *******************************************************************

function plotFoundResonances(peakFreq, bwLeftFreq, bwRightFreq)

    hold on;
   
    ymin = -20;
    ymax = 50;
    
    for n = 1:1:length(peakFreq)
        plot([peakFreq(n) peakFreq(n)], [ymin, ymax], 'k', ...
            [bwLeftFreq(n) bwLeftFreq(n)], [ymin, ymax], 'r', ...
            [bwRightFreq(n) bwRightFreq(n)], [ymin, ymax], 'r');
    end
    hold off;
end

    
% *******************************************************************
% Determine the exact frequencies and bandwidths of the resonances of the
% given spectrum (fileName) based on the given rough frequency estimates.
% *******************************************************************

function [peakFreq, bwLeftFreq, bwRightFreq] = getResonanceFreqAndBw(fileName, freqEstimates)
    peakFreq = [];
    bwLeftFreq = [];
    bwRightFreq = [];

    if (~exist(fileName, 'file'))
        return;
    end

    data = dlmread(fileName, ' ', 1, 0);    % Skip the first row.
    frequencies_Hz = data(:, 1);
    magnitudes_dB = 20.0*log10(data(:, 2));
    
    for n = 1:1:length(freqEstimates)
        % When a frequency estimate = 0 this means that we do not use it.
        if (freqEstimates(n) ~= 0)
            
            % 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), freqEstimates(n) - deltaFreq_Hz, 'nearest');
            maxSearchIndex = interp1(frequencies_Hz, 1:length(frequencies_Hz), freqEstimates(n) + deltaFreq_Hz, 'nearest');

            [dummy, peakIndex] = max(magnitudes_dB(minSearchIndex:maxSearchIndex));

            peakIndex = peakIndex + minSearchIndex - 1;

            thisPeakFreq = 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
            thisBwLeftFreq = frequencies_Hz(index);

            index = peakIndex;
            while ((index < length(frequencies_Hz)) && (magnitudes_dB(index) > peakAmp_dB - 3.0))
                index = index + 1;
            end
            thisBwRightFreq = frequencies_Hz(index);

            % Append the data of the new resonance to the lists of return
            % values.
            peakFreq = [peakFreq thisPeakFreq];
            bwLeftFreq = [bwLeftFreq thisBwLeftFreq];
            bwRightFreq = [bwRightFreq thisBwRightFreq];
        end
    end
end

