% *************************************************************************
% This script transforms the piece-wise linear area functions exported from
% VTL 2.3 to the volume-velocity transfer functions with a closed glottis
% and hard walls based on the FDS method. 
% The only losses considered here are the viscous losses at the 
% tube walls, heat conduction losses, and the radiation loss.
% The frequencies and bandwidths of the transfer functions are then 
% automatically determined and saved.
%
% This script must be run TWICE: once with 
%   useCoarseLosses = true
% and once with 
%   useCoarseLosses = false
%
% This will generate the files "fds_approx_bandwidths.m" and 
% "fds_exact_bandwidths.m" in the parent directory for the two different
% settings for the viscous, heat conduction, and radiation losses.
% *************************************************************************

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

numResonators = length(resonatorNames);

% If true, coarse (approximate) losses are used for viscous friction 
% and radiation, and the heat conduction loss is completely omitted.

global useCoarseLosses;
useCoarseLosses = false;

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

fdsResonatorsFreq_Hz = [];
fdsResonatorsBw_Hz = [];
fdsResonatorsIndex = [];        % Resonator index: 1, 2, 3, 4, ...

disp(' ');


% *******************************************************************
% Loop through all resonators.
% *******************************************************************

figure;

for n = 1:1:numResonators
    % Create a second figure for the lax vowels.
    if (n == 9)
        figure;
    end
    
    if (n <= 8)
        subplot(3, 3, n);
    else
        subplot(3, 3, (n-8));
    end
    
    % Calculate the transfer function for the current area function.
    
    [frequencies_Hz, magnitudes] = calcTransferFunction(resonatorNames{n});
    
    % Plot the magnitude spectrum in the current subplot.
    
    magnitudes_dB = 20.0*log10(abs(magnitudes));
    plot(frequencies_Hz, magnitudes_dB, 'k');

    xlim([0, 7000]);
    ylim([-20, 50]);
    title(resonatorNames{n});
    xlabel('Frequency in Hz');
    ylabel('Magnitude in dB');
    grid on;
    
    % Get the resonance frequencies and bandwidths.
    
    [peakFreq, bwLeftFreq, bwRightFreq] = getResonanceFreqAndBw(frequencies_Hz, magnitudes_dB);

    fdsResonatorsFreq_Hz = [fdsResonatorsFreq_Hz peakFreq];
    fdsResonatorsBw_Hz = [fdsResonatorsBw_Hz bwRightFreq - bwLeftFreq];
    fdsResonatorsIndex = [fdsResonatorsIndex n*ones(1, length(peakFreq))];
    
    % Output the freq. and BW on the console.
    disp([resonatorNames{n} ' :']);
    peakFreq
    bandwidths = bwRightFreq - bwLeftFreq
    
    % Mark the resonances and bandwiths in the given subplot.

    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

% *********************************************************
% End of the loop.
% *********************************************************

% Final plot with the bandwidths...

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

% Save the frequencies and bandwidths to files.

titles = {'resonator_index', 'frequency_Hz', 'bandwidth_Hz'};
data = [fdsResonatorsIndex', fdsResonatorsFreq_Hz', fdsResonatorsBw_Hz']

if (useCoarseLosses)
    writecell([titles; num2cell(data)], 'frequencies_bandwidths_coarse_losses.txt');
else
    writecell([titles; num2cell(data)], 'frequencies_bandwidths_detailed_losses.txt');
end
    
 
% *************************************************************************
% Calc. the VT transfer function based on the area function of the given
% resonator.
% *************************************************************************

function [frequencies_Hz, magnitudes] = calcTransferFunction(resonatorName)
  
    global useCoarseLosses;
    
    areasFileName = ['areas-' resonatorName '.txt'];    % Source file with the area function

    % *****************************************************
    % Load the discrete tube area function.
    % *****************************************************
    
    [x_cm A_cm2 perimeter articulator] = textread(areasFileName, '%f %f %f %s', 'headerlines', 2);

    N = length(x_cm);
    numTubeSections = N / 2;
   
    tubeLength_cm = zeros(numTubeSections, 1);
    tubeArea_cm2 = zeros(numTubeSections, 1);
    
    for n = 1:numTubeSections
        tubeLength_cm(n) = x_cm(2*n) - x_cm(2*n - 1);
        tubeArea_cm2(n) = A_cm2(2*n - 1);
    end

    % ***************************************************************
    % Calculate the transfer function using the chain matrix method.
    % The order of multiplication of the matrices is from glottis to 
    % mouth, i.e. in the order as described in Birkholz' PhD thesis,
    % and as in the paper by Wakita & Fant (1978): "Toward a better 
    % vocal tract model".
    % ***************************************************************

    % Define some physical constants equal to those in VTL.
    AMBIENT_DENSITY_CGS = 1.14e-3;   % g/cm^3
    SOUND_VELOCITY_CGS = 3.5e4;      % cm/s
    AIR_VISCOSITY_CGS = 1.86e-4;     % dyne-s/cm^2
    LAMBDA_DIV_CP_CGS = 0.229e-3;    % g/(cm*s)
    
    % Run through all frequencies from 1 to 7000 in 1-Hz steps.
    numFreqPoints = 7000;
    transferFunction = zeros(numFreqPoints);
    frequencies_Hz = 1:1:numFreqPoints;
    
    for freqIndex = frequencies_Hz
        
        omega = 2 * pi * freqIndex;
        
        % Calc. the total chain matrix for the current frequency.
        M_total = [1 0; 0 1];       % Unit matrix
    
        for n = 1:1:numTubeSections
            length_cm = tubeLength_cm(n);
            A = tubeArea_cm2(n);
            S = 2*sqrt(A * pi);     % Perimeter

            % Acoustic elements for the current tube section
            
            if (useCoarseLosses)
                % Simple resistance for laminar flow as in Maeda (1982)
                R = (length_cm/2) * (8*pi*AIR_VISCOSITY_CGS) / (A*A);
            else
                % Frequency-dependent resistance as in Flanagans book.
                R = (length_cm/2) * (S/(A*A)) * sqrt(AMBIENT_DENSITY_CGS * omega * AIR_VISCOSITY_CGS / 2);
            end
            L = AMBIENT_DENSITY_CGS * length_cm / (2*A);
            C = A*length_cm / (AMBIENT_DENSITY_CGS * SOUND_VELOCITY_CGS * SOUND_VELOCITY_CGS);
            G = length_cm * (S * 0.4 / (AMBIENT_DENSITY_CGS * (SOUND_VELOCITY_CGS^2))) * sqrt(omega * LAMBDA_DIV_CP_CGS / (2 * AMBIENT_DENSITY_CGS));

            % Z_a and Z_b of the T-type two-port-network
            Z_a = R + i*omega*L;

            if (useCoarseLosses)
                Z_b = 1 / (i*omega*C);
            else
                Z_b = 1 / (i*omega*C + G);      % Include the conductance for the heat-conduction loss.
            end

            % Chain matrix for the tube section
            M = [1 + Z_a/Z_b, 2*Z_a + Z_a*Z_a/Z_b; 1/Z_b, 1 + Z_a/Z_b];
            M_total = M_total * M;
        end
        
        % Radiation impedance at the current frequency
        Z_rad = 0;
        mouthArea_cm2 = tubeArea_cm2(numTubeSections);

        if (useCoarseLosses)
            % Use parallel circuit of R and L for the radiation (Flanagan's book).

            R_rad = (128 * AMBIENT_DENSITY_CGS * SOUND_VELOCITY_CGS) / (9 * pi * pi * mouthArea_cm2);
            L_rad = (8 * AMBIENT_DENSITY_CGS) / (3 * pi * sqrt(pi * mouthArea_cm2));
            Z_rad = (R_rad * i*omega*L_rad) / (R_rad + i*omega*L_rad);
        else
            % Use the exact radiation impedance for a piston in an infinite wall,
            % according to Birkholz' thesis, page 62.
            % There is no implementation of the Struve function in Matlab, so we
            % use a good approximation below from a JASA paper.

            ka = sqrt(mouthArea_cm2 / pi) * omega / SOUND_VELOCITY_CGS;     % wave number times the mouth radius
            Z_rad = (AMBIENT_DENSITY_CGS * SOUND_VELOCITY_CGS / mouthArea_cm2) * ...
                (1 - besselj(1, 2*ka)/ka + i*struve(2*ka)/ka);
        end
        
        transferFunction(freqIndex) = 1 / (M_total(2,1)*Z_rad + M_total(2,2));
    end
    
    magnitudes = transferFunction;
end


% *******************************************************************
% Determine the exact frequencies and bandwidths of the resonances of the
% given spectrum.
% *******************************************************************

function [peakFreq, bwLeftFreq, bwRightFreq] = getResonanceFreqAndBw(frequencies_Hz, magnitudes_dB)
    peakFreq = [];
    bwLeftFreq = [];
    bwRightFreq = [];

    % *****************************************************
    % Find the first six local maxima in the amp. spectrum.
    % *****************************************************

    numPoints = length(frequencies_Hz);
    MAX_PEAKS = 6;
    numPeaksFound = 0;
    freqIndexOfPeak = zeros(MAX_PEAKS, 1);
    
    for n = 2:1:(numPoints-1)
        if (magnitudes_dB(n) > magnitudes_dB(n-1)) && (magnitudes_dB(n) >= magnitudes_dB(n+1)) && (numPeaksFound < MAX_PEAKS)
            numPeaksFound = numPeaksFound + 1;
            freqIndexOfPeak(numPeaksFound) = n;
        end
    end
    
    % *****************************************************
    % *****************************************************
    
    for n = 1:1:numPeaksFound
        
        peakIndex = freqIndexOfPeak(n);
        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


% *******************************************************************
% Approximation equation according to the paper:
% Aarts, R. M. & Janssen, A. J. (2003). Approximation of the Struve 
% function H 1 occurring in impedance calculations. The Journal of the 
% Acoustical Society of America, 113(5), 2635-2637.
% *******************************************************************

function value = struve(x)
    value = 2/pi - besselj(0, x) + (16/pi - 5)*sin(x)/x + (12 - 36/pi)*(1 - cos(x))/(x^2);
end

