findband4.m

计算声质中的Aures Tonality

function [Ynotone,Ytone,band_index,ToneCount,NdBperOctave,err]=FindBand4(f,YdB,document,showplot,MaxToneCount,Fraction,PercentBelow)
fullscreen=get(0,'ScreenSize');
screen_stopx=1024;
screen_stopy=728;
fullscreen=[1 1+fullscreen(4)-screen_stopy-93 screen_stopx screen_stopy];
err=1;
errCount=0;

if ~exist('document')
    document='no';
end
if ~exist('showplot')
    showplot='no';
end
if ~exist('MaxToneCount')
    MaxToneCount=1;
    disp(['MaxTone has been set to' num2str(MaxTone)])
end
if ~exist('Fraction')
    Fraction=1/2;
end
if ~exist('PercentBelow')
    PercentBelow=[0.1 0.1];
end

band_index=[-1 -1 -1 -1]';
df=f(2)-f(1);
Ytone(1:length(YdB))=-100;
Yoriginal=YdB;
AX=[0 20000 -110 60];
LToneFlag=0;
RToneFlag=0;
ToneCount=0;
StepOne=3;
StepTwo=6;
BWt=3;
side=(BWt-1)/2;
if round(side)~=side
    disp('The averaging region was not selected correctly, the region has been set to:')
    side=round(side);
    BWt=2*side+1;
    disp(num2str(BWt))
end 
for ink=2:length(YdB)-side
    Yaveraged(ink)=sum(YdB(ink-side:ink+side))-10*log10(BWt);
end

for ink=1:side
    Yaveraged(ink)=YdB(ink);
end
for ink=0:side
    Yaveraged(length(YdB)-ink)=YdB(length(YdB)-ink);
end

Ysearch=Yaveraged;
Ysearcold=Ysearch;
chalk=0;
while 1
    chalk=chalk+1;
    [Ymax,PeakIndex]=max(Ysearch);
    if 0
        figure(2)
        hold off
        plot(f,Ysearch)
        hold on
        plot(f(PeakIndex),Ymax,'ro')
    end
    
    MinPeakLevel=0;
    if Ymax<MinPeakLevel
        break
    end
    ink=0;
    while 1
        ink=ink+1;
        if PeakIndex-ink<=0;
            LeftPowerIndex=1;
            break
        end
        if Yaveraged(PeakIndex-ink)+3<Ymax
            LeftPowerIndex=PeakIndex-ink;
            break
        end
        if Ysearch(PeakIndex-ink)==-100
            LeftPowerIndex=PeakIndex-ink;
            LToneFlag=1;
            break
        end
    end
    
    
   ink=0;
    while 1
        ink=ink+1;
        if PeakIndex+ink>length(Yaveraged);
           RightPowerIndex=length(Yaveraged);
            break
        end
        if Yaveraged(PeakIndex+ink)+3<Ymax
            RightPowerIndex=PeakIndex+ink;
            break
        end
        if Ysearch(PeakIndex+ink)==-100
            RightPowerIndex=PeakIndex+ink;
            RToneFlag=1;
            break
        end
    end
    if LeftPowerIndex*df<500
        StartLeftCBWIndex=LeftPowerIndex-round(Fraction*100/df);
    else
        StartLeftCBWIndex=LeftPowerIndex-round(Fraction*0.2*LeftPowerIndex);
    end
    
  if StartLeftCBWIndex<1
  StartLeftCBWIndex=1;
  end
  
  if LToneFlag==1
      LeftTonalRegionStartIndex=LeftPowerIndex;
      LTonalFlag=0;
  else
      tmp=sort(Yaveraged(StartLeftCBWIndex:LeftPowerIndex));
      ThreshIndex=floor(PercentBelow(1)*length(tmp))-1;
      if ThreshIndex<1
          ThreshIndex=1;
      end
      LeftNoiseFloor=tmp(ThreshIndex);
      
      for ink=LeftPowerIndex:-1:StartLeftCBWIndex
          if Yaveraged(ink)<=LeftNoiseFloor
              LeftTonalRegionStartIndex=ink;
              break
          end
      end
      
      
      for ink=LeftPowerIndex-2:-1:StartLeftCBWIndex
          if Yaveraged(ink)>Yaveraged(ink+1)+StepOne
              if Yaveraged(ink)>Yaveraged(ink+2)+StepTwo
                  if ink>LeftTonalRegionStartIndex
                      LeftTonalRegionStartIndex=ink;
                      [Ymin,ValleyIndex]=min(Yaveraged(LeftTonalRegionStartIndex:LeftPowerIndex));
                      LeftTonalRegionStartIndex=LeftTonalRegionStartIndex+ValleyIndex;
                      disp('A nearby tone was identified. See FindBnad_4 circa line350.')
                  end
                  break
              end
          end
      end
  end
  
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     if RightPowerIndex*df<500
        StopRightCBWIndex=RightPowerIndex+round(Fraction*100/df);
    else
        StopRightCBWIndex=RightPowerIndex+round(Fraction*0.2*RightPowerIndex);
    end
    
  if StopRightCBWIndex>length(Yaveraged)
  StopRightCBWIndex=length(Yaveraged);
  end
  
  if RToneFlag==1
      RightTonalRegionStartIndex=RightPowerIndex;
      RTonalFlag=0;
  else
      tmp=sort(YdB(RightPowerIndex:StopRightCBWIndex));
      ThreshIndex=floor(PercentBelow(2)*length(tmp))-1;
      if ThreshIndex<1
          ThreshIndex=1;
      end
      RightNoiseFloor=tmp(ThreshIndex);
      
      for ink=RightPowerIndex:+1:StopRightCBWIndex
          if YdB(ink)<=RightNoiseFloor
              RightTonalRegionStopIndex=ink;
              break
          end
      end
      
      
      for ink=RightPowerIndex+1:StopRightCBWIndex-1
          if Yaveraged(ink)<Yaveraged(ink+1)-StepOne
              if Yaveraged(ink)<Yaveraged(ink+2)-StepTwo
                  if ink<RightTonalRegionStopIndex
                      RightTonalRegionStopIndex=ink;
                      [Ymin,ValleyIndex]=min(Yaveraged(RightPowerIndex:RightTonalRegionStartIndex));
                      RightTonalRegionStartIndex=RightTonalRegionStartIndex-length(RightPowerIndex:RightTonalRegionStopIndex-ValleyIndex);
                      disp('A nearby tone was identified. See FindBnad_4 circa line420.')
                  end
                  break
              end
          end
      end
  end
    LeftSlopeIndex=LeftTonalRegionStartIndex;
    RightSlopeIndex=RightTonalRegionStopIndex;
    Ysearchold=Ysearch;
    Ysearch(LeftSlopeIndex:RightSlopeIndex)=-100;
    yLower=YdB(LeftSlopeIndex);
    yUpper=YdB(RightSlopeIndex);
    slope=(yUpper-yLower)/(RightSlopeIndex-LeftSlopeIndex);
    b=yLower;
    YdBold=YdB;
    for lead=LeftSlopeIndex:RightSlopeIndex
        YdB(lead)=slope*(lead-LeftSlopeIndex)+b;
    end
    
    fc_index=round(sqrt(LeftPowerIndex*RightPowerIndex));
    CBW_f=25+75*(1+1.4*(fc_index*df/1000)^2)^0.69;
    
    if CBW_f>(LeftPowerIndex-RightPowerIndex)*df
        [TMP7,err]=BandSumDiff(YdBold(LeftPowerIndex:RightPowerIndex),YdB(LeftPowerIndex:RightPowerIndex));
        if err==7 & 0
figure(77)
plot(YdBold(LeftPowerIndex:RightPowerIndex),'b')
    hold on
    plot(YdB(LeftPowerIndex:RightPowerIndex),'k')
    legend('Upper','Lower')
    errCount=errCount+1;
    disp(['errCount=' num2str(errCount)])
    disp(['ToneCount =' num2str(ToneCount)])
    pause
    clc
    close(77)
        end
        if err==7
            chalk=chalk-1;
        end
        if err ~=7
            Ytone(fc_index)=TMP7;
            fupper=RightPowerIndex*df;
            flower=LeftPowerIndex*df;
            dz_bark(fc_index)=13*atan(0.76*fupper/1000)+3.5*atan((fupper/7500).^2)-13*stan(0.76*flower/1000)-3.5*atan((flower/7500).^2);
            
            
            
          LeftNdBperOctaveStartIndex=LeftPowerIndex-round((LeftPowerIndex-LeftTonalRegionStartIndex)*SlopeFraction(1));
          Fl_Fu_LeftSlope=[LeftNdBOctaveStartIndex LeftPowerIndex]*df;
          RightNdBperOctaveStartIndex=round((RightTonalRegionStopIndex-RightPowerIndex)*SlopeFraction(2))+RightPowerIndex;
          Fl_Fu_RightSlope=[RightPowerIndex RightNdBperOctaveStopIndex]*df;
          Left_Octaves=(log(Fl_Fu_LeftSlope(2))-log(Fl_Fu_LeftSlope(1)))/log(2);
          Right_Octaves=(log(Fl_Fu_RightSlope(2))-log(Fl_Fu_RightSlope(1)))/log(2);
          
          LeftNdBChange=Yaveraged(LeftPowerIndex)-Yaveraged(LeftNdBperOctaveStartIndex);
           RightNdBChange=Yaveraged(RightPowerIndex)-Yaveraged(RightNdBperOctaveStartIndex);
           
            NdBperOctave=(LeftNdBChange/Left_Octaves+RightNdBChange/Right_Octaves)/2;
            
            band_index(:,chalk)=[fc_index RightPowerIndex-LeftPowerIndex+1 dz_bark(fc_index) NdBperOctave]';
            ToneCount=ToneCount+1;
        end
    end
    
    LW=1;
    MS=8;
    showplot='yes';
    disp('Show plot set to yes circa line 517 in FindBAND_V')
   if showplot=='yes'
       figure(37)
       subplot(311)
       hold off
       
      set(gcf, 'Position',fullscreen)
plot(f,Yoriginal,'g-')
set(get(gcf, 'CurrentAxes'), 'FontSize',12)
title(['From FindBband4.m',10,'Stage of Tonal Component Parameter Identification',10,'Finding Peak and 1/2 Power Points of Tone'])
ylabel('Amplitude, dB')
hold on
plot(f(LeftPowerIndex),Yoriginal(LeftPowerIndex),'b<')
           plot(f(RightPowerIndex),Yoriginal(RightPowerIndex),'b<') 
            plot(f(PeakIndex),Yoriginal(PeakIndex),'bv','LineWidth')
            
            plot(f,YdBold,'k-.')
            plot(f,Ysearchold,'r')
            zoom on
subplot(312)
hold off
plot(f,Yoriginal,'g-')
set(get(gcf,'CurrentAxes'),'FontSize',12)
title(['Removing Tone from Search'])
ylabel('Amplitude,dB')
hold on
plot(f(fc_index),Yoriginal(fc_index),'rv','LineWidth')
plot(f(StartLeftCBWIndex:LeftPowerIndex),Yoriginal(StartLeftCBWIndex:LeftPowerIndex),'k-.')
plot(f(StartLeftCBWIndex:LeftPowerIndex),LeftNoiseFloor*ones(size(StartLeftCBWIndex:LeftPowerIndex)),'m')
plot(f(LeftSlopeIndex),Yoriginal(LeftSlopeIndex),'r<')
plot(f(RightPowerIndex:StopRightCBWIndex),Yoriginal(RightPowerIndex:StopRightCBWIndex),'k-.')
plot(f(RightPowerIndex:StopRightCBWIndex),RightNoiseFloor*Ones(size(RightPowerIndex:StopRightCBWIndex)),'m')
plot(f(RightSlopeIndex),Yoriginal(RightSlopeIndex),'r>')
plot(f,Ysearch)
zoom on
subplot(312)
plot(f,Yaveraged)
   end
   Ynotone=YdB;
   
   if showplot=='yes'
       subplot(313)
hold off
plot(f,Yoriginal,'g-')
set(get(gcf,'CurrentAxes'),'FoneSize',12)
title(['Original Spectru, and Noise Spectrum'])
xlabel('Frequency, Hz')
ylabel('Amplitude,dB')
hold on
plot(f,Ynotone,'k-.')
zoom on
   end
   if ToneCount==MaxToneCount
       break
   end
end

if showplot=='yes'
    figure
    FontName='TimesNewRoman';
    FontSize=12;
    subplot(211)
    plot(f,Yaveraged,'g-')
    set(gca,'FontName',FoneName)
    set(gca,'FontSize',FontSize)
    title(['From FindBand4.m',10,'Original Spectrum and Noise Floor'])
    set(gcf,'Position',fullscreen)
    set(get(gcf,'CurrentAxes'),'FontSize',12)
hold on
plot(f,Ynotone,'k-.')
plot(f(LeftPowerIndex),Yaveraged(LeftPowerIndex),'r<')
plot(f(LeftNdNperOctaveStartIndex),Yaveraged(LeftNdBperOctaveStartIndex),'r>')
plot(f(RightNdBperOctaveStartIndex),Yaveraged(LeftNdBperOctaveStartIndex),'b<')
plot(f(RightPowerIndex),Yaveraged(RightPowerIndex),'b>')

subplot(212)
plot(f,Yaveraged,'g-')
set(gca,'FontName',FontName)
set(gca,'FontSize',FontSize)
title('Original Spectrum and Identified Tone Spectrum')
hold on
plot(f,Ytone,'b')
zoom on



if 1
    figure
    FontName='TimesNewRoman';
    FontSize=12;
    subplot(2,4,[1 2 3])
hcc=plot(f,Yaveraged);
set(hcc,'Color',[0.5 0.5 0.5])
    set(gca,'FontName',FontName)
    set(gca,'FontSize',FontSize)
    title('(a)')
    set(gcf,'Position',fullscreen)
    set(get(gcf,'CurrentAxes'),'FontSize',12)
    axis([0 7000 -40 60])
    hold on
    plot(f,Ynotone,'k:','LineWidth',2)
    
    subplot(2,4,[5 6 7])
    hcc=plot(f, Yaveraged);
    set(hcc,'Color',[0.5 0.5 0.5])
    set(gca,'FontName',FontName)
    set(gca,'FontSize',FontSize)
    title('(b)')
    axis([0 7000 -40 60])
    hold on
    plot(f,Ytone,'k-')
    print(gcf,'-depsc2','FromFindBand_V635.eps')
end
figure
subplot(3,3,[1 2 4 5])
plot(f,Yaveraged,'g-')
title(['From FindBand4.m',10,'Original Spectrum and Noise Floor'])
set(gcf,'Position',fullscreen)
set(get(gcf,'CurrentAxes'),'FontSize',12)
set(gca,'FontName','TimesNewRoman')
hold on
plot(f,Ytone,'k-.')

xlabel('Freuqncy,Hz')
ylabel('Sound Pressure Lvel,dB')

text(10,10,['Intersection with ',10,'noise floor'])
text(10,20,['Intersection with ',10,'noise floor'])
text(10,30,['Region to estimate ',10,'roll-off rate'])
text(10,40,['Region to estimate ',10,'roll-off rate'])
text(10,50,['Transition ',10,'Region'])
plot(f(LeftPowerIndex),Yaveraged(LeftPowerIndex),'r<')
plot(f(LeftNdBperOctaveStartIndex),Yaveraged(LeftNdBperOctaveStartIndex),'r>')
plot(f(RightNdBperOctaveStartIndex),Yaveraged(RightNdBperOctaveStartIndex),'b<')
plot(f(RightPowerIndex),Yaveraged(RightPowerIndex),'b>')
end
err=0;
if band_index(1,1)==-1
    err=6;
end
save Tmpfb4.mat dz_bark band_index Ytone ToneCount chalk fc_index
return