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

resonatorNames = {...
    's1-01-bahn-tense-a', ...
    's1-02-beet-tense-e', ...
    's1-03-tiere-tense-i', ...
    's1-04-boote-tense-o', ...
    's1-05-bude-tense-u', ...
    's1-06-laehmung-tense-ae', ...
    's1-07-hoehle-tense-oe', ...
    's1-08-guete-tense-y', ...
    's1-15-bass-lax-a', ...
    's1-16-bett-lax-ae', ...
    's1-17-mit-lax-i', ...
    's1-18-offen-lax-o', ...
    's1-19-butter-lax-u', ...
    's1-20-muetter-lax-y', ...
    's1-21-goetter-lax-oe', ...
    's1-22-ehe-schwa', ...
    's2-01-bahn-tense-a', ...
    's2-02-beet-tense-e', ...
    's2-03-tiere-tense-i', ...
    's2-04-boote-tense-o', ...
    's2-05-bude-tense-u', ...
    's2-06-laehmung-tense-ae', ...
    's2-07-hoehle-tense-oe', ...
    's2-08-guete-tense-y', ...
    's2-15-bass-lax-a', ...
    's2-16-bett-lax-ae', ...
    's2-17-mit-lax-i', ...
    's2-18-offen-lax-o', ...
    's2-19-butter-lax-u', ...
    's2-20-muetter-lax-y', ...
    's2-21-goetter-lax-oe', ...
    's2-22-ehe-schwa' };

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. We also excluded all resonances that
% lie directly next to an anti-resonance, because that may corrupt the
% bandwidth.
% *******************************************************************

freqEstimate = [ ...
    516, 1116, 0, 3263, 4993, 0; ... % 's1-01-bahn-tense-a'
    252, 1797, 2361, 2972, 0, 0; ... % 's1-02-beet-tense-e'
    183, 1803, 2626, 3007, 5154, 5889; ... % 's1-03-tiere-tense-i'
    259, 683, 2098, 2750, 3656, 3983; ... % 's1-04-boote-tense-o'
    194, 651, 1921, 2702, 3361, 5865; ... % 's1-05-bude-tense-u'
    438, 1599, 2239, 2939, 5397, 0; ... % 's1-06-laehmung-tense-ae'
    230, 1276, 1804, 2773, 3340, 0; ... % 's1-07-hoehle-tense-oe'
    155, 1323, 1766, 2757, 3577, 5763; ... % 's1-08-guete-tense-y'
    
    529, 968, 2380, 3030, 0, 0; ... % 's1-15-bass-lax-a'
    432, 1531, 2243, 2955, 0, 0; ... % 's1-16-bett-lax-ae'
    292, 1457, 2056, 2948, 0, 0; ... % 's1-17-mit-lax-i'
    364, 773, 2335, 3021, 3875, 4191; ... % 's1-18-offen-lax-o'
    247, 665, 2141, 2817, 0, 0; ... % 's1-19-butter-lax-u'
    301, 1163, 2055, 2852, 4078, 0; ... % 's1-20-muetter-lax-y'
    392, 1208, 2072, 3170, 3626, 0; ... % 's1-21-goetter-lax-oe'
    340, 1568, 2261, 3118, 3503, 0; ... % 's1-22-ehe-schwa'
    
    749, 1416, 3152, 3995, 4322, 0; ... % 's2-01-bahn-tense-a'
    315, 2484, 2982, 4317, 0, 0; ... % 's2-02-beet-tense-e'
    223, 2510, 3513, 4183, 0, 0; ... % 's2-03-tiere-tense-i'
    407, 904, 3169, 4059, 5178, 0; ... % 's2-04-boote-tense-o'
    333, 869, 2725, 0, 0, 0; ... % 's2-05-bude-tense-u'
    677, 1897, 2831, 4066, 0, 0; ... % 's2-06-laehmung-tense-ae'
    382, 1597, 2317, 3646, 0, 0; ... % 's2-07-hoehle-tense-oe'
    298, 1704, 2262, 3796, 0, 0; ... % 's2-08-guete-tense-y'
    
    700, 1218, 3000, 3924, 4334, 0; ... % 's2-15-bass-lax-a'
    576, 1716, 2707, 3939, 0, 0; ... % 's2-16-bett-lax-ae'
    358, 1987, 2627, 3740, 0, 0; ... % 's2-17-mit-lax-i'
    460, 1082, 3076, 4037, 0, 0; ... % 's2-18-offen-lax-o'
    378, 1083, 3006, 3843, 0, 0; ... % 's2-19-butter-lax-u'
    411, 1734, 2456, 3783, 0, 0; ... % 's2-20-muetter-lax-y'
    532, 1477, 2429, 3804, 4066, 0; ... % 's2-21-goetter-lax-oe'
    510, 1511, 2574, 3335, 4055, 0 ]; % 's2-22-ehe-schwa'


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

dvtdResonatorsFreq_Hz = [];
dvtdResonatorsBw_Hz = [];
dvtdResonatorsIndex = [];

disp(' ');

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

for n = 1:1:numResonators
    
    if ((n == 1) || (n == 9) || (n == 17) || (n == 25))
        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} '-vvtf-measured.txt'];
    plotSpectrum(fileName);

    [peakFreq, bwLeftFreq, bwRightFreq] = getResonanceFreqAndBw(fileName, freqEstimate(n, :));
    plotFoundResonances(peakFreq, bwLeftFreq, bwRightFreq);
    dvtdResonatorsFreq_Hz = [dvtdResonatorsFreq_Hz peakFreq];
    dvtdResonatorsBw_Hz = [dvtdResonatorsBw_Hz bwRightFreq - bwLeftFreq];
    dvtdResonatorsIndex = [dvtdResonatorsIndex n*ones(1, length(peakFreq))];
end

% Final plot with the bandwidths...

figure;
plot(dvtdResonatorsFreq_Hz, dvtdResonatorsBw_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 = [dvtdResonatorsIndex', dvtdResonatorsFreq_Hz', dvtdResonatorsBw_Hz']
writecell([titles; num2cell(data)], '../dvtd_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([-30, 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 = -30;
    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.
            MAX_BW_HZ = 400;
            index = peakIndex;
            
            while ((index > 1) && (index > peakIndex - MAX_BW_HZ/2) && (magnitudes_dB(index) > peakAmp_dB - 3.0))
                index = index - 1;
            end
            thisBwLeftFreq = frequencies_Hz(index);

            leftSideValid = true;
            if (magnitudes_dB(index) > peakAmp_dB - 3.0)
                leftSideValid = false;
            end

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

            rightSideValid = true;
            if (magnitudes_dB(index) > peakAmp_dB - 3.0)
                rightSideValid = false;
            end

            % If the measurement at the left or right sid is invalid, then
            % rely on the measurement of the other side only.
            
            if (~leftSideValid)
                thisBwLeftFreq = thisPeakFreq - (thisBwRightFreq - thisPeakFreq);
            end
            
            if (~rightSideValid)
                thisBwRightFreq = thisPeakFreq + (thisPeakFreq - thisBwLeftFreq);
            end
            
            % 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

  