tone_iii.m

计算声质中的Aures Tonality

function [K,K_mod,w1,w1_mod,w2,w3,Wt,Wgr,TLX,Tdz,NoiseN,YdBN,NoiseNs,YdBNs]=Tone_III(fname,field,SPLwant,document,showplot,MaxToneCount,fname_path,Fraction,PercentBelow)
close all
if nargin<1
    clc
    disp(['Error,not enough input arguments were given.'])
    return
end
if ~exist('field')
    field='f';
    disp(['Loudness was calculated with MS se to' field])
end
if ~exist ('SPLwant')
    SPLwant=0;
    dsip('No Calibration will be done')
end
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('fname_path')
    fname_path=pwd;
end
if ~exist('Fraction')
    Fraction=1/2;
end
if ~exist('PercentBelow')
    PercentBelow=0.1;
end

fullscreen=[1 1 1500 970]
TF = strcmp(document, 'yes')
if TF==1
    diary Test.txt
    disp(['File Nmae=' fname])
    diary off
end
[y,fs,nbots,err]=LoadTone(fname,fname_path);
if err~=0;
    disp('There was an error in loading the signal. The program will quit now')
    disp(['The error had a value of ',num2str(err)])
    help LoadTone.m
    return
end
if SPLwant~=0;
    [ycal,err]=Calibrate(y,SPLwant);
    y=ycal;
end
if 1
    if exist('fname_path')
        end_path=length(fname_path);
        if fname_path(end_path-1)=='G'
            disp('Clibration based on folder name. see Tone_II.m ~ line 161')
            GAIN=str2num(fname_path(end_path))
            y=y*GAIN;
        end
            if fname_path(end_path-2)=='G'
            disp('Clibration based on folder name. see Tone_II.m ~ line 161')
            GAIN=str2num(fname_path(end_path-1:end_path))
            y=y*GAIN;
            end
            if fname_path(end_path-3)=='G'
            disp('Clibration based on folder name. see Tone_II.m ~ line 161')
            GAIN=str2num(fname_path(end_path-2:end_path))
            y=y*GAIN;
            end
    end
end
df=1;
[Yxx,f]=PowSpecdF(y,fs,df);
df=f(2)-f(1);
[YdB,err]=Convert2dB(Yxx,1);
NoiseFloor=-30;
for ink=1:length(YdB)
    if YdB(ink)<NoiseFloor
        YdB(ink)=NoiseFloor;
    end
end
for ink =1:floor(20/df)
    YdB(ink)=NoiseFloor;
end
if err~=0;
    disp(['There was an error in calculating the Calibrated SPL''The program will quit now'])
    disp(['The error had a value of',num2str(err)])
help Convert2dB_ALH.m
return
end

[Ynotone,Ytone,band_index,ToneCount,NdBperOctave,err]=FindBand4(f,YdB,document,showplot,MaxToneCount,Fraction,PercentBelow);
if err~=0;
    if err~=2;
        disp(['There was an unknown error in finding narrow band components.''The program will quit now'])
        disp(['The error had a value of',num2str(err)])
        help FindBand_V.m
        return
    end
end

TONAL=Ytone;
NOISE=Ynotone;
TonalIndex=band_index(1,:);
if TonalIndex(1)==-1
    disp('No Tones were found')
    K=O;
    figure
    plot(d,YdB)
    title(['Your sound has no tones in it!',10,'From Tone.m'])
    xlabel('Frequency,Hz')
    ylabel('SPL,dB')
    return
end

[RelIndex,err]=ExcessLevel(TonalIndex,TONAL,NOISE,f,document,showplot);
if showplot=='yes'
    set(figure,'Position',fullscreen)
    subplot(211)
    plot(TONAL,'g')
    hold on
    plot(NOISE,'-.k')
    title(['Tones and Noise for' fname,10,'Plotted from Tone.m '])
    subplot(212)
    plot(RelIndex(1,:))
    title(['Excess Level Calculated for' fname])
saves(gcf,fname(1:length(fname)-4),'fig')
end

[Tf_index,TBW_index,Tdz,TLX,err]=ArrangeComponents(band_index,RelIndex);
if document=='yes'
    diary Test.txt
    disp(['Tonal Frequencies, Tf_index= ' num2str(Tf_index)])
    dsip(['Tonal Bnadwidths, TBW_index=' num2str(TBW_index)])
    dsip(['Tonal Roll-Off Rate=' num2str(band_index(4,:))])
     dsip(['OriginalTonal Delta Z, Tdz=' num2str(Tdz)])
      dsip(['Tonal Excess Level, TLX=' num2str(TLX)])
       dsip(['df=' num2str(f(2)-f(1))])
       diary off
end


if err ~=0
    disp(['There was an error in building relevant vector components' 'The program will quit'])
    disp(['The error had a value of ',num2str(err)])
    help ArrangeComponents.m
    return
end


for ink=1:length(Tdz)
    if band_index(2:ink0)==5
        Tdz(ink)=0;
    else
        flower=df*(Tf_index(ink)-2);
        fupper=df*(Tf_index(ink)+2);
        Zupper=13*atan(0.76*fupper/1000)+3.5*atan((fupper/7500).^2);
        Zlower=13*atan(0.76*flower/1000)+3.5*atan((flower/7500).^2);
        dzCorrect(ink)=(Zupper-Zlower);
        Tdz(ink)=(Tdz(ink)-dzCorrect(ink)-1);
    end
end

if document=='yes'
    diary Test.txt
    disp(['Corrected Tonal Delta Z, Tdz=' num2str(Tdz)])
    diary off
end
downstairs=Tdz+0.13;
if any(downstairs<0.13)
    for fixx=1:length(downstairs)
        if downstairs(fixx)<0.13
            downstairs(fixx)=0.13;
        end
    end
end


w1=0.13./downstairs;


RR=band_index(4,:);

NAN_Seek=isnan(RR);
NAN_Locations=find(NAN_Seek==1);
RR(NAN_Locations)=999;
PBW=Tdz;
Tf_kHz=Tf_index*df/1000;
zc=21.4*log10(4.37*Tf_kHz+1);
if 1
    a=2.4360;
    b=1.7769;
    c=49.0577;
    d=0.9425;
    w1_mod=(a./(a+PBW)).*(1-exp(-RR./(c+d*zc))).^b;
else
    a=2.4458;
    b=1.7313;
    c=64.5609;
   w1_mod=(a./(a+PBW)).*(1-exp(-RR./c)).^b;
end

Tf1=(df/1000)*(Tf_index-1);
downstairs=(1+0.2*((Tf1./0.7)+0.7./Tf1).^2).^(1/2);
pow=0.29;
disp(['The Exponent for the W2 equation is set to' num2str(pow)])
w2=(1./downstairs).^pow;
MaxW2=(1/sqrt(0.8))^pow;
w2=w2+(1-MaxW2);
w3=(1-exp(-1*TLX/15)).^(pow);
MultWp=(w1.^(1/pow)).*(w2.^(1/pow)).*(w3.^(1/pow));
SqWp=MultWp.^2;
SumWp=sum(SqWp);
Wt=(SumWp)^(1/2);
MultWp_mod=(w1_mod.^(1/pow)).*(w2.^(1/pow)).*(w3.^(1/pow));
SqWp_mod=MultWp_mod.^2;
SumWp_mod=sum(SqWp_mod);
Wt_mod=(SumWp_mod)^(1/2);

clear TONAL
TONAL=zeros(1,size(f,1));
TONAL(RelIndex(2,:))=RelIndex(1,:);

MaxFilt=28;
[H,err]=GenerateFilters(f,fs);

H=0.9639*H;

NoiseZinput=zeros(1,MaxFilt);
for ink=1:MaxFilt
    NoiseZinput(ink)=10*log10(sum((10.^(NOISE/10)).*(abs(H(ink,:).^2))));
end

[NoiseN,NoiseNs,err]=DIN45631(NoiseZinput,field);


YdBZinput=zeros(1,MaxFilt);
for ink=1:MaxFilt
    YdBZinput(ink)=10*log10(sum((10.^(YdB/10)).*(abs(H(ink,:).^2))));
end

[YdBN,YdBNs,err]=DIN45631(YdBZinput,field);

Wgr=1-NoiseN/YdBN;
if Wgr<0
    Wgr=0;
end

if pow==0.29
    
    c=1.0090;
    
elseif pow==0.45
    c=1.0129;
else c=1;
    disp('Warning, final calibration has not been made, c=1')
end

K=c*(Wt^pow)*(Wgr^0.79);



K_mod=c*(Wt_mod^pow)*(Wgr^0.79);


if showplot=='yes'
    
    v=[0 2000 -5 80]
    set(figure,'Position', fullscreen)
    subplot(511)
plot(YdB)
set(get(gcf,'CurrentAxes'),'FontSize',174)
title(['YdB (from PSD) vs Index'])
xlabel('Index, (Frequency/df)')
ylabel('SPL,dB')
subplot(513)
plot(Ytone)
set(get(gcf,'SurrentAxes'),'FontSize',14)
title(['TONAL (Before Excess Level Asjustment) vs Index'])
xlabel('Index, (Frequency/df)')
ylabel('SPL dB')
subplot(515)
plot(NOISE)
set(get(gcf,'CurrentAxes'),'FoneSize',14)
title(['NOISE vs. Index'])
xlabel('Index, (Frequency/df)')
ylabel('SPL,dB')
zoom pm

z=0.1:0.1:24
set(figure,'Position',fullscreen)
subplot(211)
plot(z,YdBNs)
set(get(gcf,'SurrentAxes'),'FontSize',14)
title(['Specific Loudness of Total Spectrum'])
xlabel('Critical Band Rate, Bark')
ylabel('Specific Loudness, Sones/Bark')
subplot(212)
plot(z,NoiseNs)
set(get(gcf,'SurrentAxes'),'FontSize',14)
title(['Specific Loudness of Noise Spectrum'])
xlabel('Critical Band Rate, Bark')
ylabel('Specific Loudness,',10, 'Sones/Bark')
zoom on



set(figure,'Position',fullscreen)
subplot(311)
plot(f,YdB)
set(get(gcf,'SurrentAxes'),'FontSize',14)
title(['Amplitude Spectrum of Signal'])
xlabel('Frequency,Hz')
ylabel('SPL,dB')
grid
subplot(313)
plot(z,YdBNs)
set(get(gcf,'SurrentAxes'),'FontSize',14)
title(['Specific Loudness of Signal',10,'Total Loudness=',num2str(YdBN)])
xlabel('Critical Band Rate, Bark')
ylabel('Specific Loudness,',10, 'Sones/Bark')
grid
zoom on
end

if document=='yes'
    diary Test.txt
    disp(['File Name=',fname])
    disp(['w1=',num2str(w1)])
 disp(['w1_mod=',num2str(w1_mod)])
  disp(['w2=',num2str(w2)])
   disp(['w3=',num2str(w3)])
    disp(['NoiseN=',num2str(NoiseN)])
     disp(['YdBN=',num2str(YdBN)])
 disp(['Wt=',num2str(Wt)])
  disp(['Wgr=',num2str(Wgr)])
   disp(['K=',num2str(K)])
   diary off
end 
return