excesslevel.m

计算声质中的Aures Tonality

function [RelIndex,err]= ExcessLevel(TonalIndex,TONAL,NOISE,f,document,showplot)
fullscreen=get(0,'ScreenSize');
screen_stopx=1024;
screen_stopy=728;
fullscreen=[1 1+fullscreen(4)-screen_stopy-93 screen_stopx screen_stopy];
err=1;
z=13*atan(0.76*f./1000)+3.5*atan((f/7500).^2);

L=zeros(1,size(TONAL,2));
L(TonalIndex)=TONAL(TonalIndex);
Lth=zeros(1,size(TONAL,2));
Lth(TonalIndex)=3.64*(f(TonalIndex)/1000).^(-0.8)-6.5*exp(-0.6*((f(TonalIndex)/1000)-3.3).^2)+(10^(-3))*(f(TonalIndex)/1000).^4;
Ehs=zeros(1,size(TONAL,2));
Ehs(TonalIndex)=10.^(Lth(TonalIndex)/10);

if showplot=='yes'
    set(figure,'Position',fullscreen)
    subplot(211);
    plot(f,L)
    set(get(gcf,'CurrentAxes'),'FontSize',18)
    title(['From ExcessLevel',10,'Main Excitation Level vs. Frequency'])
    xlabel('Frequency, Hz')
    ylabel('Level,dB')
    subplot(212)
    plot(f,Lth)
   set(get(gcf,'CurrentAxes'),'FontSize',18)
   title(['From Excess Level',10,'Threshold of Hearing, dB vs. Frequency'])
   xlabel('Frequency,Hz')
   ylabel('Level,dB')
end

Egr=zeros(1,size(TONAL,2));
for i=1:size(TonalIndex,2)
    Zcenter=z(TonalIndex(i));
    j=1;
    while 1
        if TonalIndex(i)-j<=0
            j=j-1;
            Zlower=z(TonalIndex(i)-j);
            break
        end
        Zlower=z(TonalIndex(i)-j);
        if Zlower>Zcenter-0.5
            j=j+1;
        else 
            break
        end
    end
    k=1;
    while 1
        if TonalIndex(i)+k>length(z)
            k=k-1;
            Zupper=z(TonalIndex(i)+k);
            break
        end
        Zupper=z(TonalIndex(i)+k);
        if Zupper<Zcenter+0.5
            k=k+1;
        else 
            Zdiff(i)=Zupper-Zlower;
            if Zdiff(i)>1.2
                disp(['Warning, the sum of noise for term' num2str(i) 'was conducted over' num2str(Zdiff(i)) 'Bark'])
                disp('This is most likely due to a component that is close to the beginning or end of the frequency spectrum')
                disp('This warning was generated in the ExcessLevel function')
            end
            if Zdiff(i)<0.8
                     
                disp(['Warning, the sum of noise for term' num2str(i) 'was conducted over' num2str(Zdiff(i)) 'Bark'])
                disp('This is most likely due to a component that is close to the beginning or end of the frequency spectrum')
                disp('This warning was generated in the ExcessLevel function')
            end
                break
        end 
    end
    
    sortTonalIndex=sort(TonalIndex);
    for i=1;length(TonalIndex)
        for j=1:i-1
            s=-24-230/f(sortTonalIndex(j))+0.2*L(sortTonalIndex(j));
            Lek(i,j)=L(sortTonalIndex(j))-s*(z(sortTonalIndex(j))-z(sortTonalIndex(i)));
        end
        
        for j=i+1:size(TonalIndex,2)
            s=27;
             Lek(i,j)=L(sortTonalIndex(j))-s*(z(sortTonalIndex(j))-z(sortTonalIndex(i)));
        end
    end
    if length(TonalIndex)~=1
        Aek=10.^(Lek/20);
    else
        Lek=0;
        Aek=0;
    end
    
    AekSUM=zeros(1,size(TONAL,2));
    AekSUM(sortTonalIndex)=sum(Aek');
    LX=zeros(1,size(TONAL,2));
    LX(sortTonalIndex)=L(sortTonalIndex)-10*log10(AekSUM(sortTonalIndex).^2+Egr(sortTonalIndex)+Ehs(sortTonalIndex));
    
    
    if 1
        set(figure,'Position',fullscreen)
        MS=6;
        LW=0.5;
        FontName='TimesNewRoman';
        FontSize=16;
        subplot(3,1,1:2)
        plot(f(sortTonalIndex),L(sortTonalIndex),'go')
        hold on
        warning off
        plot(f(sortTonalIndex),10*log10(Ehs(sortTonalIndex)),'bd')
          plot(f(sortTonalIndex),10*log10(Egr(sortTonalIndex)),'k')
      plot(f(sortTonalIndex),10*log10(AekSUM(sortTonalIndex).^2),'mx')
      warning on
      set(gca,'FontName',FontName)
      set(gca,'FontSize',FontSize)
      title(['Components Used to Generate Excess Level' 10 , 'This graph produced in Function ExcessLevel circa line 200'])
      legend('Tones','Audidbility Threshold','Noise','Accesory Masking')
      subplot(313)
      plot(f(sortTonalIndex),LX(sortTonalIndex),'go')
      set(gca,'FontName',FontName)
      set(gca,'FontSize',FontSize)
      title('Excess Level')
    end
    %%%%%%%%%%%%%
        
    if 1
        set(figure,'Position',fullscreen)
        MS=4;
        LW=0.5;
        FontName='TimesNewRoman';
        FontSize=12;
        subplot(4,4,[1 2 3 5 6 7])
        plot(f(sortTonalIndex),L(sortTonalIndex),'go')
        hold on
        warning off
        plot(f(sortTonalIndex),10*log10(Ehs(sortTonalIndex)),'bd')
          plot(f(sortTonalIndex),10*log10(Egr(sortTonalIndex)),'k')
      plot(f(sortTonalIndex),10*log10(AekSUM(sortTonalIndex).^2),'mx')
      warning on
      set(gca,'FontName',FontName)
      set(gca,'FontSize',FontSize)
      %title(['Components Used to Generate Excess Level' 10 , 'This graph produced in Function ExcessLevel circa line 200'])
      title('(a)')
      
      ylabel('Level,dB')
      legend('Tones','Audidbility Threshold','Noise','Accesory Masking')
      subplot(3,4,[9 10 11])
      plot(f(sortTonalIndex),LX(sortTonalIndex),'go')
      set(gca,'FontName',FontName)
      set(gca,'FontSize',FontSize)
      title('(b)')
      ylabel('Level,dB')
      xlabel('Frequency,Hz')
      print(gcf,'-deps2','FromExcessLevel_Line248.eps')
    end

     if document=='yes'
         diary Test.txt
         temp=sort(L);
         temp=fliplr(temp);
         disp(['Main Excitation,L=' num2str(temp(1:1))])
         temp=sort(Lth)
         temp=fliplr(temp);
         disp(['Threshold Level, Lth=' num2str(temp(1:1))])
       temp=sort(Egr)
         temp=fliplr(temp);
         disp(['Masking by noise, Egr(Indensity)=' num2str(temp(1:1))])
          temp=sort(AekSUM)
         temp=fliplr(temp);
         disp(['Masking by Acessory Excitation,AekSUM(Intensity)' num2str(temp(1:1))])
         diary off
     end
     
     
     
     if showplot=='yes'
         set(figure,'Position',fullscreen)
         subplot(211)
plot(f(sortTonalIndex),L(sortTonalIndex),'b^')
hold on
if AekSUM(sortTonalIndex)~=0
    plot(f(sortTonalIndex),10*log10(AekSUM(sortTonalIndex)),'rv')
end
set(get(gcf,'CurentAxes'),'FontSize',18)
title(['From ExcessLevel',10,'Excitation Level vs. f'])
xlabel('Frequency,Hz')
ylabel('Level,dB')
legend('Main Excitation Level','Accessory Excitation Level')
axis([min(f) max(f) 0 150])
subplot(212)
plot(f, LX)
set(get(gcf,'CurentAxes'),'FontSize',18)
title(['Excess Level vs. Frequency'])
xlabel('Frequency, Hz')
ylabel('Level,dB')
hold off
     end
     
     RelIndex=[0 0]';
     for i =1:size(TonalIndex,2)
         if LX(sortTonalIndex(i)>0)
             RelIndex(:,i)=[LX(sortTonalIndex(i)) sortTonalIndex(i)]';
         else
             RelIndex(:,i)=[0 sortTonalIndex(i)]';
         end
     end
     err=0;
     return