clear all
close all

load dirchar_men_s.mat

vowphonid = [1 6 11]; % phonation modes [breathy,modal,pressed]
numofphon = length(vowphonid);

singerid = [1:4]; % singer
numofpers = length(singerid);

load windowsize.mat
spath = '../../figs_new/';
%% order in data
% freq,vowel,singer,pitch

pitch_idd  = 7:13;

numofpitch = length(pitch_idd);

engv_h_all_sel = engv_h_all(:,:,:,pitch_idd);
engv_v_all_sel = engv_v_all(:,:,:,pitch_idd);
HDI3_sel = HDI3(:,:,:,pitch_idd);
VDI3_sel = VDI3(:,:,:,pitch_idd);

engv_h_all_sel(find(engv_h_all_sel==0)) = NaN;
engv_v_all_sel(find(engv_v_all_sel==0)) = NaN;

HDI3_sel(find(HDI3_sel==0)) = NaN;
VDI3_sel(find(VDI3_sel==0)) = NaN;

pitch_Ref_Freqs = [123 131 147 165 175 196 220 247 262 294 330 350 392 440];
pitche_Strs =    {'H','c','d','e','f','g','a','h','c1','d1','e1','f1','g1','a1'};

low_xlim = 200;

%% HORIZONTAL AND VERTICAL BEAMWIDTH
if 1
    if numofpitch == 1
        pdata_a = reshape(engv_h_all_sel(:,vowphonid,singerid),[size(engv_h_all_sel,1),numofphon*numofpers*numofpitch]);
        pdata_e = reshape(engv_h_all_sel(:,vowphonid+1,singerid),[size(engv_h_all_sel,1),numofphon*numofpers*numofpitch]);
        pdata_i = reshape(engv_h_all_sel(:,vowphonid+2,singerid),[size(engv_h_all_sel,1),numofphon*numofpers*numofpitch]);
        pdata_o = reshape(engv_h_all_sel(:,vowphonid+3,singerid),[size(engv_h_all_sel,1),numofphon*numofpers*numofpitch]);
        pdata_u = reshape(engv_h_all_sel(:,vowphonid+4,singerid),[size(engv_h_all_sel,1),numofphon*numofpers*numofpitch]);
        tname = sprintf('%d Hz %s%',pitch_Ref_Freqs(pitch_idd),pitche_Strs{pitch_idd});
    else
        pdata_a = reshape(engv_h_all_sel(:,vowphonid,singerid,1:length(pitch_idd)),[size(engv_h_all_sel,1),numofphon*numofpers*numofpitch]);
        pdata_e = reshape(engv_h_all_sel(:,vowphonid+1,singerid,1:length(pitch_idd)),[size(engv_h_all_sel,1),numofphon*numofpers*numofpitch]);
        pdata_i = reshape(engv_h_all_sel(:,vowphonid+2,singerid,1:length(pitch_idd)),[size(engv_h_all_sel,1),numofphon*numofpers*numofpitch]);
        pdata_o = reshape(engv_h_all_sel(:,vowphonid+3,singerid,1:length(pitch_idd)),[size(engv_h_all_sel,1),numofphon*numofpers*numofpitch]);
        pdata_u = reshape(engv_h_all_sel(:,vowphonid+4,singerid,1:length(pitch_idd)),[size(engv_h_all_sel,1),numofphon*numofpers*numofpitch]);
        tname = sprintf('%d Hz to %d Hz, and %s to %s',pitch_Ref_Freqs(pitch_idd(1)),pitch_Ref_Freqs(pitch_idd(end)),pitche_Strs{pitch_idd(1)},pitche_Strs{pitch_idd(end)})
    end
    pdata_a_mean=mean(pdata_a,2,'omitnan');
    SD_a=std(pdata_a,0,2,'omitnan');
    BSD_a=(pdata_a_mean-SD_a);
    TSD_a=(pdata_a_mean+SD_a);

    pdata_e_mean=mean(pdata_e,2,'omitnan');
    SD_e=std(pdata_e,0,2,'omitnan');
    BSD_e=(pdata_e_mean-SD_e);
    TSD_e=(pdata_e_mean+SD_e);

    pdata_i_mean=mean(pdata_i,2,'omitnan');
    SD_i=std(pdata_i,0,2,'omitnan');
    BSD_i=(pdata_i_mean-SD_i);
    TSD_i=(pdata_i_mean+SD_i);

    pdata_o_mean=mean(pdata_o,2,'omitnan');
    SD_o=std(pdata_o,0,2,'omitnan');
    BSD_o=(pdata_o_mean-SD_o);
    TSD_o=(pdata_o_mean+SD_o);

    pdata_u_mean=mean(pdata_u,2,'omitnan');
    SD_u=std(pdata_u,0,2,'omitnan');
    BSD_u=(pdata_u_mean-SD_u);
    TSD_u=(pdata_u_mean+SD_u);

    figure
    patch([F(2:end); flipud(F(2:end))],[BSD_u(2:end); flipud(TSD_u(2:end))],'r','Facealpha',0.1,'EdgeColor',[.8 .8 .8])
    hold on;

    patch([F(2:end); flipud(F(2:end))],[BSD_o(2:end); flipud(TSD_o(2:end))],'g','Facealpha',0.1,'EdgeColor',[.8 .8 .8])

    patch([F(2:end); flipud(F(2:end))],[BSD_i(2:end); flipud(TSD_i(2:end))],'m','Facealpha',0.1,'EdgeColor',[.8 .8 .8])

    patch([F(2:end); flipud(F(2:end))],[BSD_e(2:end); flipud(TSD_e(2:end))],'b','Facealpha',0.1,'EdgeColor',[.8 .8 .8])

    patch([F(2:end); flipud(F(2:end))],[BSD_a(2:end); flipud(TSD_a(2:end))],'k','Facealpha',0.1,'EdgeColor',[.8 .8 .8])
    semilogx(F,pdata_u_mean,'r','LineWidth',3); 
    semilogx(F,pdata_o_mean,'g','LineWidth',3); 
    semilogx(F,pdata_i_mean,'m','LineWidth',3); 
    semilogx(F,pdata_e_mean,'b','LineWidth',3); 
    semilogx(F,pdata_a_mean,'k','LineWidth',3); 

    grid on
    set(gca,'XScale','log','FontSize',14)
    set(gcf,'Position',cpos)
    axis tight
    xlabel('Frequency in Hz')
    ylabel('Hor. beam width in degree')
    xlim([low_xlim 10000])
    ylim([80 190])
    box on
    legend('u','o','i','e','a','u','o','i','e','a','Location','northwest','Orientation','horizontal')
    title(tname)

    print([spath 'male_Hbeam'],'-depsc')
    %% VERTICAL BEAMWIDTH
    if numofpitch == 1
        pdata_a = reshape(engv_v_all_sel(:,vowphonid,singerid),[size(engv_v_all_sel,1),numofphon*numofpers*numofpitch]);
        pdata_e = reshape(engv_v_all_sel(:,vowphonid+1,singerid),[size(engv_v_all_sel,1),numofphon*numofpers*numofpitch]);
        pdata_i = reshape(engv_v_all_sel(:,vowphonid+2,singerid),[size(engv_v_all_sel,1),numofphon*numofpers*numofpitch]);
        pdata_o = reshape(engv_v_all_sel(:,vowphonid+3,singerid),[size(engv_v_all_sel,1),numofphon*numofpers*numofpitch]);
        pdata_u = reshape(engv_v_all_sel(:,vowphonid+4,singerid),[size(engv_v_all_sel,1),numofphon*numofpers*numofpitch]);
        tname = sprintf('%d Hz %s%',pitch_Ref_Freqs(pitch_idd),pitche_Strs{pitch_idd});
    else
        pdata_a = reshape(engv_v_all_sel(:,vowphonid,singerid,1:length(pitch_idd)),[size(engv_v_all_sel,1),numofphon*numofpers*numofpitch]);
        pdata_e = reshape(engv_v_all_sel(:,vowphonid+1,singerid,1:length(pitch_idd)),[size(engv_v_all_sel,1),numofphon*numofpers*numofpitch]);
        pdata_i = reshape(engv_v_all_sel(:,vowphonid+2,singerid,1:length(pitch_idd)),[size(engv_v_all_sel,1),numofphon*numofpers*numofpitch]);
        pdata_o = reshape(engv_v_all_sel(:,vowphonid+3,singerid,1:length(pitch_idd)),[size(engv_v_all_sel,1),numofphon*numofpers*numofpitch]);
        pdata_u = reshape(engv_v_all_sel(:,vowphonid+4,singerid,1:length(pitch_idd)),[size(engv_v_all_sel,1),numofphon*numofpers*numofpitch]);
        tname = sprintf('%d Hz to %d Hz, and %s to %s',pitch_Ref_Freqs(pitch_idd(1)),pitch_Ref_Freqs(pitch_idd(end)),pitche_Strs{pitch_idd(1)},pitche_Strs{pitch_idd(end)})
    end
    pdata_a_mean=mean(pdata_a,2,'omitnan');
    SD_a=std(pdata_a,0,2,'omitnan');
    BSD_a=(pdata_a_mean-SD_a);
    TSD_a=(pdata_a_mean+SD_a);

    pdata_e_mean=mean(pdata_e,2,'omitnan');
    SD_e=std(pdata_e,0,2,'omitnan');
    BSD_e=(pdata_e_mean-SD_e);
    TSD_e=(pdata_e_mean+SD_e);

    pdata_i_mean=mean(pdata_i,2,'omitnan');
    SD_i=std(pdata_i,0,2,'omitnan');
    BSD_i=(pdata_i_mean-SD_i);
    TSD_i=(pdata_i_mean+SD_i);

    pdata_o_mean=mean(pdata_o,2,'omitnan');
    SD_o=std(pdata_o,0,2,'omitnan');
    BSD_o=(pdata_o_mean-SD_o);
    TSD_o=(pdata_o_mean+SD_o);

    pdata_u_mean=mean(pdata_u,2,'omitnan');
    SD_u=std(pdata_u,0,2,'omitnan');
    BSD_u=(pdata_u_mean-SD_u);
    TSD_u=(pdata_u_mean+SD_u);
    %
    figure
    patch([F(2:end); flipud(F(2:end))],[BSD_u(2:end); flipud(TSD_u(2:end))],'r','Facealpha',0.1,'EdgeColor',[.8 .8 .8])
    hold on;

    patch([F(2:end); flipud(F(2:end))],[BSD_o(2:end); flipud(TSD_o(2:end))],'g','Facealpha',0.1,'EdgeColor',[.8 .8 .8])

    patch([F(2:end); flipud(F(2:end))],[BSD_i(2:end); flipud(TSD_i(2:end))],'m','Facealpha',0.1,'EdgeColor',[.8 .8 .8])

    patch([F(2:end); flipud(F(2:end))],[BSD_e(2:end); flipud(TSD_e(2:end))],'b','Facealpha',0.1,'EdgeColor',[.8 .8 .8])

    patch([F(2:end); flipud(F(2:end))],[BSD_a(2:end); flipud(TSD_a(2:end))],'k','Facealpha',0.1,'EdgeColor',[.8 .8 .8])
    semilogx(F,pdata_u_mean,'r','LineWidth',3); 
    semilogx(F,pdata_o_mean,'g','LineWidth',3); 
    semilogx(F,pdata_i_mean,'m','LineWidth',3); 
    semilogx(F,pdata_e_mean,'b','LineWidth',3); 
    semilogx(F,pdata_a_mean,'k','LineWidth',3); 

    grid on
    set(gca,'XScale','log','FontSize',14)
    set(gcf,'Position',cpos)
    axis tight
    xlabel('Frequency in Hz')
    ylabel('Vert. beam width in degree')
    xlim([low_xlim 10000])
    ylim([80 190])
    box on
    legend('u','o','i','e','a','u','o','i','e','a','Location','northwest','Orientation','horizontal')
    title(tname)
    
    print([spath 'male_Vbeam'],'-depsc')
end
if 1
    %% HDI & VDI
    if numofpitch == 1
        pdata_a = reshape(HDI3_sel(:,vowphonid,singerid),[size(HDI3_sel,1),numofphon*numofpers*numofpitch]);
        pdata_e = reshape(HDI3_sel(:,vowphonid+1,singerid),[size(HDI3_sel,1),numofphon*numofpers*numofpitch]);
        pdata_i = reshape(HDI3_sel(:,vowphonid+2,singerid),[size(HDI3_sel,1),numofphon*numofpers*numofpitch]);
        pdata_o = reshape(HDI3_sel(:,vowphonid+3,singerid),[size(HDI3_sel,1),numofphon*numofpers*numofpitch]);
        pdata_u = reshape(HDI3_sel(:,vowphonid+4,singerid),[size(HDI3_sel,1),numofphon*numofpers*numofpitch]);
        tname = sprintf('%d Hz %s%',pitch_Ref_Freqs(pitch_idd),pitche_Strs{pitch_idd});
    else
        pdata_a = reshape(HDI3_sel(:,vowphonid,singerid,1:length(pitch_idd)),[size(HDI3_sel,1),numofphon*numofpers*numofpitch]);
        pdata_e = reshape(HDI3_sel(:,vowphonid+1,singerid,1:length(pitch_idd)),[size(HDI3_sel,1),numofphon*numofpers*numofpitch]);
        pdata_i = reshape(HDI3_sel(:,vowphonid+2,singerid,1:length(pitch_idd)),[size(HDI3_sel,1),numofphon*numofpers*numofpitch]);
        pdata_o = reshape(HDI3_sel(:,vowphonid+3,singerid,1:length(pitch_idd)),[size(HDI3_sel,1),numofphon*numofpers*numofpitch]);
        pdata_u = reshape(HDI3_sel(:,vowphonid+4,singerid,1:length(pitch_idd)),[size(HDI3_sel,1),numofphon*numofpers*numofpitch]);
        tname = sprintf('%d Hz to %d Hz, and %s to %s',pitch_Ref_Freqs(pitch_idd(1)),pitch_Ref_Freqs(pitch_idd(end)),pitche_Strs{pitch_idd(1)},pitche_Strs{pitch_idd(end)})
    end

    pdata_a_mean=mean(pdata_a,2,'omitnan');
    SD_a=std(pdata_a,0,2,'omitnan');
    BSD_a=(pdata_a_mean-SD_a);
    TSD_a=(pdata_a_mean+SD_a);

    pdata_e_mean=mean(pdata_e,2,'omitnan');
    SD_e=std(pdata_e,0,2,'omitnan');
    BSD_e=(pdata_e_mean-SD_e);
    TSD_e=(pdata_e_mean+SD_e);

    pdata_i_mean=mean(pdata_i,2,'omitnan');
    SD_i=std(pdata_i,0,2,'omitnan');
    BSD_i=(pdata_i_mean-SD_i);
    TSD_i=(pdata_i_mean+SD_i);

    pdata_o_mean=mean(pdata_o,2,'omitnan');
    SD_o=std(pdata_o,0,2,'omitnan');
    BSD_o=(pdata_o_mean-SD_o);
    TSD_o=(pdata_o_mean+SD_o);

    pdata_u_mean=mean(pdata_u,2,'omitnan');
    SD_u=std(pdata_u,0,2,'omitnan');
    BSD_u=(pdata_u_mean-SD_u);
    TSD_u=(pdata_u_mean+SD_u);

    figure
    patch([F(2:end); flipud(F(2:end))],[BSD_u(2:end); flipud(TSD_u(2:end))],'r','Facealpha',0.1,'EdgeColor',[.8 .8 .8])
    hold on;

    patch([F(2:end); flipud(F(2:end))],[BSD_o(2:end); flipud(TSD_o(2:end))],'g','Facealpha',0.1,'EdgeColor',[.8 .8 .8])

    patch([F(2:end); flipud(F(2:end))],[BSD_i(2:end); flipud(TSD_i(2:end))],'m','Facealpha',0.1,'EdgeColor',[.8 .8 .8])

    patch([F(2:end); flipud(F(2:end))],[BSD_e(2:end); flipud(TSD_e(2:end))],'b','Facealpha',0.1,'EdgeColor',[.8 .8 .8])

    patch([F(2:end); flipud(F(2:end))],[BSD_a(2:end); flipud(TSD_a(2:end))],'k','Facealpha',0.1,'EdgeColor',[.8 .8 .8])
    semilogx(F,pdata_u_mean,'r','LineWidth',3); 
    semilogx(F,pdata_o_mean,'g','LineWidth',3); 
    semilogx(F,pdata_i_mean,'m','LineWidth',3); 
    semilogx(F,pdata_e_mean,'b','LineWidth',3); 
    semilogx(F,pdata_a_mean,'k','LineWidth',3); 

    grid on
    set(gca,'XScale','log','FontSize',14)
    set(gcf,'Position',cpos)
    axis tight
    xlim([low_xlim 10000])
    ylim([-8 15])
    box on
    xlabel('Frequency in Hz')
    ylabel('HDI')
    legend('u','o','i','e','a','u','o','i','e','a','Location','southwest','Orientation','horizontal')
    title(tname)
    
    print([spath 'male_HDI'],'-depsc')
    %% VDI
    if numofpitch == 1
        pdata_a = reshape(VDI3_sel(:,vowphonid,singerid),[size(VDI3_sel,1),numofphon*numofpers*numofpitch]);
        pdata_e = reshape(VDI3_sel(:,vowphonid+1,singerid),[size(VDI3_sel,1),numofphon*numofpers*numofpitch]);
        pdata_i = reshape(VDI3_sel(:,vowphonid+2,singerid),[size(VDI3_sel,1),numofphon*numofpers*numofpitch]);
        pdata_o = reshape(VDI3_sel(:,vowphonid+3,singerid),[size(VDI3_sel,1),numofphon*numofpers*numofpitch]);
        pdata_u = reshape(VDI3_sel(:,vowphonid+4,singerid),[size(VDI3_sel,1),numofphon*numofpers*numofpitch]);
        tname = sprintf('%d Hz %s%',pitch_Ref_Freqs(pitch_idd),pitche_Strs{pitch_idd});
    else
        pdata_a = reshape(VDI3_sel(:,vowphonid,singerid,1:length(pitch_idd)),[size(VDI3_sel,1),numofphon*numofpers*numofpitch]);
        pdata_e = reshape(VDI3_sel(:,vowphonid+1,singerid,1:length(pitch_idd)),[size(VDI3_sel,1),numofphon*numofpers*numofpitch]);
        pdata_i = reshape(VDI3_sel(:,vowphonid+2,singerid,1:length(pitch_idd)),[size(VDI3_sel,1),numofphon*numofpers*numofpitch]);
        pdata_o = reshape(VDI3_sel(:,vowphonid+3,singerid,1:length(pitch_idd)),[size(VDI3_sel,1),numofphon*numofpers*numofpitch]);
        pdata_u = reshape(VDI3_sel(:,vowphonid+4,singerid,1:length(pitch_idd)),[size(VDI3_sel,1),numofphon*numofpers*numofpitch]);
        tname = sprintf('%d Hz to %d Hz, and %s to %s',pitch_Ref_Freqs(pitch_idd(1)),pitch_Ref_Freqs(pitch_idd(end)),pitche_Strs{pitch_idd(1)},pitche_Strs{pitch_idd(end)})
    end

    pdata_a_mean=mean(pdata_a,2,'omitnan');
    SD_a=std(pdata_a,0,2,'omitnan');
    BSD_a=(pdata_a_mean-SD_a);
    TSD_a=(pdata_a_mean+SD_a);

    pdata_e_mean=mean(pdata_e,2,'omitnan');
    SD_e=std(pdata_e,0,2,'omitnan');
    BSD_e=(pdata_e_mean-SD_e);
    TSD_e=(pdata_e_mean+SD_e);

    pdata_i_mean=mean(pdata_i,2,'omitnan');
    SD_i=std(pdata_i,0,2,'omitnan');
    BSD_i=(pdata_i_mean-SD_i);
    TSD_i=(pdata_i_mean+SD_i);

    pdata_o_mean=mean(pdata_o,2,'omitnan');
    SD_o=std(pdata_o,0,2,'omitnan');
    BSD_o=(pdata_o_mean-SD_o);
    TSD_o=(pdata_o_mean+SD_o);

    pdata_u_mean=mean(pdata_u,2,'omitnan');
    SD_u=std(pdata_u,0,2,'omitnan');
    BSD_u=(pdata_u_mean-SD_u);
    TSD_u=(pdata_u_mean+SD_u);

    figure
    patch([F(2:end); flipud(F(2:end))],[BSD_u(2:end); flipud(TSD_u(2:end))],'r','Facealpha',0.1,'EdgeColor',[.8 .8 .8])
    hold on;

    patch([F(2:end); flipud(F(2:end))],[BSD_o(2:end); flipud(TSD_o(2:end))],'g','Facealpha',0.1,'EdgeColor',[.8 .8 .8])

    patch([F(2:end); flipud(F(2:end))],[BSD_i(2:end); flipud(TSD_i(2:end))],'m','Facealpha',0.1,'EdgeColor',[.8 .8 .8])

    patch([F(2:end); flipud(F(2:end))],[BSD_e(2:end); flipud(TSD_e(2:end))],'b','Facealpha',0.1,'EdgeColor',[.8 .8 .8])

    patch([F(2:end); flipud(F(2:end))],[BSD_a(2:end); flipud(TSD_a(2:end))],'k','Facealpha',0.1,'EdgeColor',[.8 .8 .8])
    semilogx(F,pdata_u_mean,'r','LineWidth',3); 
    semilogx(F,pdata_o_mean,'g','LineWidth',3); 
    semilogx(F,pdata_i_mean,'m','LineWidth',3); 
    semilogx(F,pdata_e_mean,'b','LineWidth',3); 
    semilogx(F,pdata_a_mean,'k','LineWidth',3); 

    grid on
    set(gca,'XScale','log','FontSize',14)
    set(gcf,'Position',cpos)
    axis tight
    xlim([low_xlim 10000])
    legend('u','o','i','e','a','u','o','i','e','a','Location','northwest','Orientation','horizontal')
    ylim([-8 15])
    box on
    xlabel('Frequency in Hz')
    ylabel('VDI')
    title(tname)
    
    print([spath 'male_VDI'],'-depsc')
end