% *******************************************************************
% Simulation of soft walls with a low-frequency approximation of the 
% vocal tract for the paper "Acoustic comparison of physical vocal 
% tract models with hard and soft walls".
% Copyright (C) by Peter Birkholz, TU Dresden, Germany, 2021
% *******************************************************************

global AIR_DENSITY;
global SOUND_VELOCITY;

% *******************************************************************
% The calculations here are based on the FDS in Peter Birkholz' paper:
% "Influence of temporal discretization schemes on formant frequencies 
% and bandwidths in time domain simulations of the vocal tract system"
% *******************************************************************

% Everything is in CGS units !!
AIR_DENSITY = 0.00114;      % g / cm^3
SOUND_VELOCITY = 35000.0;   % cm / s

numFreqPoints = 1000;   % Number of frequency points
numTubes = 8;           % Number of different tube models
frequencies_Hz = 1:1:numFreqPoints;
magnitudeHard = zeros(numTubes, numFreqPoints);
magnitudeSoft1 = zeros(numTubes, numFreqPoints);
magnitudeSoft2 = zeros(numTubes, numFreqPoints);

% *******************************************************************
% Define the geometry of the two sections of each model.
% *******************************************************************

% Both tube section lengths are 8 cm for all models.
length1_cm = 8 * ones(numTubes, 1);
length2_cm = 8 * ones(numTubes, 1);
area1_cm2 = zeros(numTubes, 1);
area2_cm2 = zeros(numTubes, 1);

% Area of the first tube section decreases, and that of the 2nd increases.
for n = 1:1:numTubes
  area1_cm2(n) = 5.657 * (1/sqrt(2))^(n-1);
  area2_cm2(n) = 0.5 * sqrt(2)^(n-1);
end

% *******************************************************************
% Determine and plot the volume velocity transfer functions.
% *******************************************************************

resonanceFreqHard_Hz = zeros(numTubes, 1);
resonanceBwHard_Hz = zeros(numTubes, 1);

resonanceFreqSoft1_Hz = zeros(numTubes, 1);
resonanceBwSoft1_Hz = zeros(numTubes, 1);

resonanceFreqSoft2_Hz = zeros(numTubes, 1);
resonanceBwSoft2_Hz = zeros(numTubes, 1);

for tubeIndex = 1:1:numTubes
    for n = 1:1:numFreqPoints
        magnitudeHard(tubeIndex, n) = getTransferFunction(frequencies_Hz(n), ...
            length1_cm(tubeIndex), area1_cm2(tubeIndex), length2_cm(tubeIndex), area2_cm2(tubeIndex), 0);
        
        magnitudeSoft1(tubeIndex, n) = getTransferFunction(frequencies_Hz(n), ...
            length1_cm(tubeIndex), area1_cm2(tubeIndex), length2_cm(tubeIndex), area2_cm2(tubeIndex), 1);

        magnitudeSoft2(tubeIndex, n) = getTransferFunction(frequencies_Hz(n), ...
            length1_cm(tubeIndex), area1_cm2(tubeIndex), length2_cm(tubeIndex), area2_cm2(tubeIndex), 2);
    end
    
    % Obtain the data about the main resonance.
    [resonanceFreqHard_Hz(tubeIndex), resonanceBwHard_Hz(tubeIndex)] = getResonanceData(magnitudeHard(tubeIndex, :));
    [resonanceFreqSoft1_Hz(tubeIndex), resonanceBwSoft1_Hz(tubeIndex)] = getResonanceData(magnitudeSoft1(tubeIndex, :));
    [resonanceFreqSoft2_Hz(tubeIndex), resonanceBwSoft2_Hz(tubeIndex)] = getResonanceData(magnitudeSoft2(tubeIndex, :));

    % Plot the transfer functions.
    subplot(3, 3, tubeIndex);
    plot(frequencies_Hz, 20*log10(magnitudeHard(tubeIndex, :)), 'k', ...
        frequencies_Hz, 20*log10(magnitudeSoft1(tubeIndex, :)), 'r', ...
        frequencies_Hz, 20*log10(magnitudeSoft2(tubeIndex, :)), 'g');
    hold on;
    min_dB = -20;
    max_dB = 50;
    plot([resonanceFreqHard_Hz(tubeIndex), resonanceFreqHard_Hz(tubeIndex)], [min_dB max_dB], 'k');
    plot([resonanceFreqSoft1_Hz(tubeIndex), resonanceFreqSoft1_Hz(tubeIndex)], [min_dB max_dB], 'r');
    plot([resonanceFreqSoft2_Hz(tubeIndex), resonanceFreqSoft2_Hz(tubeIndex)], [min_dB max_dB], 'g');
    hold off;
    ylim([min_dB, max_dB]);
    xlabel('Frequency in Hz');
    ylabel('Magnitude in dB');
    grid on;
end


% *******************************************************************
% Plot the resonance frequency with hard vs. soft walls.
% *******************************************************************

subplot(3, 3, 7);
plot(resonanceFreqHard_Hz, resonanceFreqSoft1_Hz - resonanceFreqHard_Hz, 'r-o', ...
    resonanceFreqHard_Hz, resonanceFreqSoft2_Hz - resonanceFreqHard_Hz, 'g-o');
minFreq_Hz = 150;
maxFreq_Hz = 850;
xlim([minFreq_Hz, maxFreq_Hz]);
ylim([0, 80]);
xlabel('f_{R1,hard} in Hz');
ylabel('f_{R1,soft} - f_{R1,hard} in Hz');
grid on;

% *******************************************************************
% Plot the bandwidth for hard vs. soft walls.
% *******************************************************************

subplot(3, 3, 9);
plot(resonanceFreqHard_Hz, resonanceBwHard_Hz, 'k-o', ...
    resonanceFreqSoft1_Hz, resonanceBwSoft1_Hz, 'r-o', ...
    resonanceFreqSoft2_Hz, resonanceBwSoft2_Hz, 'g-o');
xlim([minFreq_Hz, maxFreq_Hz]);
xlabel('f_{R1} in Hz');
ylabel('Bandwidth in Hz');
grid on;


% *******************************************************************
% Determine the frequency and the bandwidth of the first (and only)
% resonance. The spectrum has a frequency resolution of 1 Hz.
% *******************************************************************

function [resonanceFreq_Hz, resonanceBw_Hz] = getResonanceData(spectrum)
    [maxValue, resonanceFreq_Hz] = max(spectrum(1:1000));
    threshold = maxValue / sqrt(2.0);       % = 3 dB below the maximum
    
    leftIndex = resonanceFreq_Hz;
    while (spectrum(leftIndex) > threshold)
        leftIndex = leftIndex - 1;
    end
    
    rightIndex = resonanceFreq_Hz;
    while (spectrum(rightIndex) > threshold)
        rightIndex = rightIndex + 1;
    end
    
    resonanceBw_Hz = rightIndex - leftIndex;
end


% *******************************************************************
% wallSoftness = 0 -> hard walls
% wallSoftness = 1 -> first setting for soft walls.
% wallSoftness = 2 -> second setting for soft walls.
% *******************************************************************

function magnitude = getTransferFunction(freq_Hz, length1_cm, area1_cm2, length2_cm, area2_cm2, wallSoftness)

    global AIR_DENSITY;
    global SOUND_VELOCITY;
    
    omega = 2 * pi * freq_Hz;
    
    Z_L1 = i * omega * AIR_DENSITY * length1_cm / (2 * area1_cm2);
    Z_L2 = i * omega * AIR_DENSITY * length2_cm / (2 * area2_cm2);
    Z_C1 = 1 / (i * omega * length1_cm * area1_cm2 / (AIR_DENSITY * SOUND_VELOCITY * SOUND_VELOCITY));
    Z_C2 = 1 / (i * omega * length2_cm * area2_cm2 / (AIR_DENSITY * SOUND_VELOCITY * SOUND_VELOCITY));
    Z_rad = getRadiationImpedance(freq_Hz, area2_cm2);

    K_total = [1 Z_L1; 0 1];    % Left inductivity of section 1.

    K = [1 0; 1/Z_C1 1];        % Capacity of section 1.
    K_total = K_total * K;
    
    % Use soft walls for section 1?
    if (wallSoftness > 0)
        Z_wall1 = getWallImpedance(freq_Hz, length1_cm, area1_cm2, wallSoftness);
        K = [1 0; 1/Z_wall1 1];
        K_total = K_total * K;
    end
    
    K = [1 (Z_L1 + Z_L2); 0 1]; % Right inductivity of section 1 and left inductivity of section 2.
    K_total = K_total * K;
    
    K = [1 0; 1/Z_C2 1];        % Capacity of section 2.
    K_total = K_total * K;
    
    % Use soft walls for section 2?
    if (wallSoftness > 0)
        Z_wall2 = getWallImpedance(freq_Hz, length2_cm, area2_cm2, wallSoftness);
        K = [1 0; 1/Z_wall2 1];
        K_total = K_total * K;
    end

    K = [1 Z_L2; 0 1];          % Right inductivity for section 2.
    K_total = K_total * K;

    % Final value of the transfer function.
    magnitude = 1.0 / (K_total(2,1) * Z_rad + K_total(2,2));
    magnitude = abs(magnitude);
end


% *******************************************************************
% wallSoftness may be 1 or 2 for two different settings of material
% parameters.
% *******************************************************************

function Z_wall = getWallImpedance(freq_Hz, sectionLength_cm, sectionArea_cm2, wallSoftness)
    
    if (wallSoftness == 1)
        % Values from the paper of Flanagan (1975)
        % Curves plotted in RED.
        B = 1600; %1600   % g / (cm^2*s)
        M = 1.5; %1.5;    % g / cm^2
    else
        % Values from the Stevens book, p. 157.
        % Curves plotted in GREEN.
        B = 1000;    % g / (cm^2*s)
        M = 2.0;    % g / cm^2
    end
    
    K = 84500;  % g / (m^2 * s^2)
    
    perimeter_cm = 2*sqrt(sectionArea_cm2 * pi);
    innerSurface_cm2 = sectionLength_cm * perimeter_cm;
    R = B / innerSurface_cm2;
    L = M / innerSurface_cm2;
%    C = innerSurface_cm2 / K;
    omega = 2 * pi * freq_Hz;
    Z_wall = R + i*omega*L; % + 1.0 / (i*omega*C);
end

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

function Z_rad = getRadiationImpedance(freq_Hz, area_cm2)
    global AIR_DENSITY;
    global SOUND_VELOCITY;

    R_rad = (128 * AIR_DENSITY * SOUND_VELOCITY) / (9 * pi^2 * area_cm2);
    L_rad = (8 * AIR_DENSITY) / (3 * pi * sqrt(pi * area_cm2));

    omega = 2 * pi * freq_Hz;
    Z_rad = (R_rad * i*omega*L_rad) / (R_rad + i*omega*L_rad);
end

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


