function NormLogFormants = CrossLangNormalise(GroupingIndices, RawLogFormants)
%
% function NormLogFormants = CrossLangNormalise(GroupingIndices, RawLogFormants)
%
% Geoffrey Stewart Morrison
% Release 2008-12-29 (almost identical to Release 2007-04-20)
%
% Cross-language / cross-dialect version of constant log interval normalisation of vowel formant values.
% Also normalises speakers' L2 vowels on the basis of their L1 data.
%
% Procedure described in:
%       Morrison, G. S., & Nearey, T. M.  (in prep). A simple pattern recognition model and cross-language vowel normalization 
%           for prediction of cross-language vowel perception. 
% Variant of procedure described in:
%       Morrison, G. S., & Nearey, T. M.  (2006). A cross-language vowel normalisation procedure. 
%           Canadian Acoustics, 34(3), 94–95.
%
%
% INPUT ARGUMENTSS:
%   GroupingIndices:    A matrix of numeric grouping indices. Rows are tokens. 
%                       Columns are [lang1, lang_spoken, gender, speaker, vowel_cat]
%                                       1       2           3       4       5
%
%                       lang1       indicates speaker's first language or dialect
%
%                       lang_spoken indicates language being spoken (optional column, default value = lang1)
%                                       Tokens with lang_spoken ~= lang1 are not used to calculate speaker and language normalisation
%                                       factors (Gspeaker & Glang), but are normalised using the Gspeaker & Glang values calculated using
%                                       the lang_spoken = lang1 data. Use lang_spoken to specify test tokens such as the speakers' L2 vowels.
%
%                       speaker     ids must be unique, i.e., use combinations such as
%                                       [lang1=1 speaker=1; lang1=1 speaker=2; ... lang1=2 speaker=11; lang1=2 speaker=12]
%                                       not [lang1=1 speaker=1; lang1=1 speaker=2; ... lang1=2 speaker=1; lang1=2 speaker=2]
%
%                       vowel_cat   ids must be unique
%
%   RawLogFormants:     A matrix of formant values (in log Hertz). Rows are tokens. 
%                           Columns are measures taken at points on formant tracks, e.g., [F1@25 F1@75 F2@25 F2@75]
%                           (use same number of points from each formant)
%                       OR
%                       A unidimensional cell array. Cells are tokens.
%                           Each cell contains a matrix of formant tracks, e.g., [F1_track F2_track F3_track]
%                           (time can be along either dimension)
%
% OUTPUT ARGUMENT:
%   NormLogFormants:    Normalised formant values (in log Hertz) arranged in same format as RawLogFormants.


if nargin == 0 
    % sample input data (only one token per vowel cat per speaker, only two male and two female speakers per lang)
    % speaker: 1=Spanish 2=English      gender: 0=Female 1=Male
    % vowel: 1=Sp/i/ 2=Sp/ei/ 3=Sp/e/ 4=Eng/i/ 5=Eng/I/ 6=Eng/e/ 7=Eng/E/
    GroupingIndices =...
    [1     0   206     1
     1     0   206     2
     1     0   206     3
     1     0   209     1
     1     0   209     2
     1     0   209     3
     1     1   200     1
     1     1   200     2
     1     1   200     3
     1     1   201     1
     1     1   201     2
     1     1   201     3
     2     0   165     4
     2     0   165     5
     2     0   165     6
     2     0   165     7
     2     0   167     4
     2     0   167     5
     2     0   167     6
     2     0   167     7
     2     1   166     4
     2     1   166     5
     2     1   166     6
     2     1   166     7
     2     1   175     4
     2     1   175     5
     2     1   175     6
     2     1   175     7];
    RawLogFormants =...
   [5.8859    5.8547    7.9167    7.9151
    6.3490    6.0536    7.7448    7.8548
    6.2899    6.2698    7.7153    7.7413
    5.9041    5.8424    7.8288    7.8694
    6.3401    5.9541    7.7615    7.8931
    6.1871    6.2079    7.7229    7.7252
    5.8675    5.8075    7.5662    7.5898
    6.2724    5.9672    7.4233    7.5649
    6.2422    6.2359    7.3461    7.3526
    5.7206    5.6800    7.5838    7.5709
    6.0621    5.7661    7.4527    7.6007
    6.0514    6.0679    7.4051    7.4289
    5.6702    5.7970    7.8323    7.8269
    6.2297    6.4688    7.6098    7.5711
    6.1980    6.1327    7.7110    7.7680
    6.5402    6.6139    7.5188    7.4597
    5.9205    5.9394    7.9543    7.9476
    6.4793    6.5030    7.5792    7.5800
    6.3094    6.1912    7.8575    7.9044
    6.6902    6.7667    7.5510    7.5219
    5.9513    5.7251    7.6815    7.6926
    6.1013    6.1083    7.5013    7.4530
    6.1748    6.0765    7.5503    7.6127
    6.3397    6.4229    7.4091    7.3083
    5.6346    5.6906    7.6608    7.6606
    6.1442    6.1783    7.3817    7.3173
    6.0755    6.0236    7.5348    7.6898
    6.2842    6.3597    7.3283    7.3083];
end


% Gtoken: mean of formant values for each token
if iscell(RawLogFormants)
    Gtoken = cellfun(@mean,RawLogFormants,'UniformOutput',0);
    Gtoken = cellfun(@mean,Gtoken,'UniformOutput',1);
else
    Gtoken = mean(RawLogFormants,2); % Gtoken
end

% if lang_spoken is not provided, assume pure cross-language situation (assume missing column is lang_spoken)
if size(GroupingIndices,2) == 4 
    GroupingIndices = [GroupingIndices(:,1), GroupingIndices(:,1), GroupingIndices(:,2:4)];
end

% CALCULATE NORMALISATION FACTORS
%   Gspeaker:   mean of formant values for each speaker
%   Glang:      mean of formant values for each language 

langs = unique(GroupingIndices(:,1));     num_langs = length(langs);
genders = unique(GroupingIndices(:,3));   num_genders = length(genders);
speakers = unique(GroupingIndices(:,4));  num_speakers = length(speakers);

II_gender = cell(num_genders,1); % using cells because ther may be unequal numbers of speakers per gender etc.
for Igender = 1:num_genders
    II_gender{Igender} = find(GroupingIndices(:,3) == genders(Igender));
end

II_speaker = cell(num_speakers,1);
for Ispeaker=1:num_speakers
    II_speaker{Ispeaker}=  find(GroupingIndices(:,4) == speakers(Ispeaker));
end

Glang_gender = zeros(num_langs,num_genders);
Gspeaker = zeros(num_speakers,1);

II_lang = cell(num_langs,1);
for Ilang = 1:num_langs
    II_lang{Ilang} = find(GroupingIndices(:,1) == langs(Ilang) & GroupingIndices(:,2) == langs(Ilang)); % use L1 data only
    
    for Igender = 1:num_genders % gender within language
        II_lang_gender = intersect(II_lang{Ilang}, II_gender{Igender}); 
        speakers_sub = unique(GroupingIndices(II_lang_gender,4));
        num_speakers_sub = length(speakers_sub);
        Gspeaker_sub = zeros(num_speakers_sub,1);
        
        for Ispeaker_sub = 1:num_speakers_sub % speaker within gender within language
            Ispeaker = find(speakers == speakers_sub(Ispeaker_sub));
            II_lang_gender_speaker = intersect(II_lang_gender,II_speaker{Ispeaker});
            vowels_sub = unique(GroupingIndices(II_lang_gender_speaker,5));
            num_vowels_sub = length(vowels_sub);
            Gvowel_sub = zeros(num_vowels_sub,1);
            
            for Ivowel_sub = 1:num_vowels_sub % vowel cat within speaker within gender within language
                II_vowel_sub = GroupingIndices(II_lang_gender_speaker,5) == vowels_sub(Ivowel_sub);
                Gvowel_sub(Ivowel_sub) = mean(Gtoken(II_lang_gender_speaker(II_vowel_sub)));
            end
            
            Gspeaker(Ispeaker) = mean(Gvowel_sub); % Gspeaker
            Gspeaker_sub(Ispeaker_sub) = Gspeaker(Ispeaker);
        end
        
        Glang_gender(Ilang,Igender)=mean(Gspeaker_sub); % Glang_gender (intermediate step to Glang)
    end 
end

Glang = mean(Glang_gender,2); % Glang


% NORMALISE FORMANT DATA
if iscell(RawLogFormants)
    NormLogFormants=cell(size(RawLogFormants));
else
    NormLogFormants=zeros(size(RawLogFormants));
end

for Ispeaker=1:num_speakers
    Ilang = langs == GroupingIndices(II_speaker{Ispeaker}(1),1);
    % Ntoken = Gtoken - Gspeaker + Glang
    if iscell(RawLogFormants)
        NormLogFormants(II_speaker{Ispeaker}) = cellfun(@(x) x - Gspeaker(Ispeaker) + Glang(Ilang), RawLogFormants(II_speaker{Ispeaker}), 'UniformOutput', 0);
    else
        NormLogFormants(II_speaker{Ispeaker},:) = RawLogFormants(II_speaker{Ispeaker},:) - Gspeaker(Ispeaker) + Glang(Ilang);
    end
end

return