% Tumbuka intonation
% (adapted from Louis Goldstein's script on Mambila downdrift)
% 
% Miran Oh
% 1 April 2020

% Model the results in (11b) and (12) of Downing's article

% Idea:
% In Tumbuka, intonational phrase edges are marked by a low tone boundary
% and phonological phrases by a high tone boundary (& internal-IP with H%)
% 
% A. 
% The range of High tone pitch raising decreases from beginning to the end
% -> decreasing attractive force for H tone (h)
%
% B.
% Depending on phrasing, the internal pitch patterns differ
% e.g., pitch register is lowered when followed by an H% tone 
% -> stronger attractive force for L boundary tone? (l)
% 
% C. 
% H% blocks the downdrift 
% (not incorporated in the model yet)

clear 

% Set the initial f0 value
T(1) = 90;

% Set PP & IP boundary tones
H = 10;         % the increase in pitch for H tone in phonological phrases
L_ = T(1) - 10; % L boundary "tone" target (goal for the declination system)
H_ = T(1) + 10; % H boundary tone target

k  = .3;    % stiffness of the declination systems (individual diff?)
h0 = 1;     % scaling of attractive force for H tone (H and H%)
l0 = 0.3;   % scaling of attractive force for L tone (for additional L% effect)

nsyll = 9;
PP = [3,9]; % syllable location of phonological phrases
IP=[];

for i = 2:nsyll  
    h = ((nsyll-i)./nsyll)*h0; % change the degree of attractive force over time
    l = (i./nsyll)*l0; % stronger force towards the end
    
    if ismember(i, PP-1) 
        dT = -k*T(i-1) + k*L_ + h.*(H); % penult High tone
    elseif ismember(i, IP)
        dT = -k*T(i-1) + k*L_ + h.*(H_-T(i-1)); % IP
    elseif isempty(IP)
        dT = -k*T(i-1) + k*L_; % dynamics for single IP
    else
        dT = -k*T(i-1) + k*L_ - l.*(T(i-1)-L_); % (multiple IP)
    end
    T(i) = T(i-1) + dT;
end
k1=k; T1=T

% add sentence-internal IP
IP = [3];     % syllable location of sentence-internal IP phrases

for i = 2:nsyll  
    h = ((nsyll-i)./nsyll)*h0; % change the degree of attractive force over time
    l = (i./nsyll)*l0; % stronger force towards the end
    
    if ismember(i, PP-1) 
        dT = -k*T(i-1) + k*L_ + h.*(H); % penult High tone
    elseif ismember(i, IP)
        dT = -k*T(i-1) + k*L_ + h.*(H_-T(i-1)); % IP
    elseif isempty(IP)
        dT = -k*T(i-1) + k*L_; % dynamics for single IP
    else
        dT = -k*T(i-1) + k*L_ - l.*(T(i-1)-L_); % (multiple IP)
    end
    T(i) = T(i-1) + dT;
end

syll = 1:nsyll;

k2=k; T2=T
subplot(1,2,1)
plot (syll, T1, '-o', 'LineWidth' ,2);
title(['k = ', num2str(k1)]);ylim([80 100])
ylabel({'Pitch'}); xlabel({'n. of syllables'});

subplot(1,2,2)
plot (syll, T2, '-o', 'LineWidth' ,2);
line([IP(1)+.5 IP(1)+.5], [80 100]);
title(['k = ', num2str(k2), ',   IP after syllable ',num2str(IP)]);ylim([80 100]);
ylabel({'Pitch'}); xlabel({'n. of syllables'});
