b_duration.m

我认为很不错的语音处理的ｍａｔｌａｂ源代码

function [b_min, b, p2]=B_duration(p, Fs)
% % B_duration: Calculates the B-duration for impulsive noise
% % 
% % Syntax:  [b_min, b, p2]=B_duration(p, Fs);
% %
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% % Description 
% % 
% % This program calculates the B-duration for impulsive noise analysis
% % 
% % Reference: Guido F. Smoorenburg, "Damage Risk Criteria for Impulsive
% %            Noise," New Perspectives on Noise Induced Hearing Loss, 
% %            Raven Press, New York, pages(471-490) 1982
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % Input variables
% % 
% % p is the sound pressure in Pa for a single channel data array
% %                 the default value is randn(1, 50000);
% %
% % Fs sampling rate in Hz.  default value is 100000 Hz.  
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% % Output variables
% %
% % b_min is the smallest estimate of the B-duration 
% % It is recommended to use b_min as the estimated B-duration.
% % 
% % b is the array of estimated B-durations for each threshold value
% % thresholds form 1 dB to 20 dB are used to make the esimate more robust.
% % Some thresholds typically will have anomalies and overestimate the 
% % B-duration.  
% %
% % p2 is the processed data array.
% % The processing involvse the running average of the absolute value of 
% % the hilbert transform.
% % 
% % This quantity, p2, is the time average instantaneous sound pressure.
% %
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Example;
%
% % Example impulsive data with background noise
%
% Fs=50000; fc=200; td=1; tau=0.1; delay=0.1; A1=3; A=20;
% [p, t]=analytic_impulse(Fs, fc, td, tau, delay, A1, A2);
% % p               % Pa sound pressure, single channel data array.
%                   % p should have only one impulse.
% Fs=50000;         % Hz sample rate frequency
% 
% [b_min, b, p2]=B_duration(p, Fs);
% 
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% % Sub Programs
% % 
% % sub_mean        Removes the time varying dc offset.  
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% % B_duration.m was originally developed by Chucri Kardous.  
% %
% % This implementation of B_duration was written by Edward L. Zechmann  
% % 
% %      date 11 December  2007
% % 
% %  modified 17 December  2007  Added Comments 
% % 
% %  modified 13 August    2008  Updated Comments
% % 
% %  modified 21 September 2008  Check output for being empty
% %                              Updated Comments
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% % Please feel free to modify this code.
% % 
% % See also: A_Duration, C_Duration, D_Duration
% % 

if (nargin < 1 || isempty(p)) || ~isnumeric(p)
    p=randn(1, 50000);
end

if nargin <2 
    Fs=50000;
end

[buf, p2]=sub_mean(abs(hilbert(p)), Fs, 2000);

% array of threshold values in dB 
threshold_a=[1, 2, 4, 8, 10, 15, 20]; % dB
num_thres=length(threshold_a);

% B-duration (assuming only one impulse)

[maxp maxp_index]=max(abs(p2));
b=zeros(num_thres, 1);
num_data=length(p);

for e2=1:num_thres;
    thres=10^(-threshold_a(e2)/20);
    maxp_thres=maxp*thres;
    x_thres_index=find(p2 > maxp_thres);

    % find start time for B-duration, which is the first 20 dB threshold-crossing before
    % the peak pressure
    flag1=0;
    e1=x_thres_index(1)+1;

    if (x_thres_index(1)+1) > 1 && (x_thres_index(1)+1) < length(p2)
        e1=x_thres_index(1)+1;
    else
        e1=1;
    end

    while (flag1 == 0) && (e1 > 1)
        e1 = e1-1;
        if p2(e1) <= maxp_thres
            flag1=1;
        end
    end

    % interpolate to find better begin time
    if abs(p2(e1)-p2(e1+1)) >= 10^-12
        bt1 = (maxp_thres-p2(e1))/(p2(e1+1)-p2(e1))+e1;
    else
        bt1=e1;
    end

    % find end time for B-duration, i.e. the last 20 dB threshold-crossing after peak
    flag1=0;

    if (x_thres_index(end)-1) > 1 && ((x_thres_index(end)-1) <= length(p))
        e1=(x_thres_index(end)-1);
    else
        e1=2;
    end

    while (flag1 == 0) && (e1 < length(p))
        e1 = e1+1;
        if p2(e1) <= maxp_thres
            flag1=1;
        end
    end

    % interpolate to find better end time
    if abs(p2(e1)-p2(e1-1)) >= 10^-12
        bt2 = (maxp_thres-p2(e1-1))/(p2(e1)-p2(e1-1))+e1-1;
    else
        bt2=e1;
    end

    % B-duration in indices
    b(e2) = (bt2-bt1);

end

% normaize b by the 20 dB threshold  and convert to seconds
b=1/Fs*(20./threshold_a)'.*b;

% Recommended to use the minimum of the estimated B-durations
% Background noise inceases the B-duration 
% Flat or jumpy peaks can increase the B-duration 
% Both of the above affects excessively increase the B-duration
% Also make sure that the B-duraiton is positive
b_min=min(b(b > 0.00000001));

% Make sure that B-duration is not longer than the data array
if b_min > 1/Fs*num_data
    b_min=1/Fs*num_data;
end

if isempty(b)
    b=-1;
    b_min=-1;
end