d_duration.m

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

function [d, ddd, ee, tau_fit, y2, y4]=D_duration(p, Fs, make_plot)
% % D_duration: Calculates the D-duration for impulsive noise
% % 
% % Syntax:
% % 
% % [d, ddd, ee, tau_fit, y2, y4]=D_duration(p, Fs, make_plot);
% % 
% % *********************************************************************
% % 
% % Description
% % 
% % This program calculates the D-duration for impulsive noise analysis.
% % This program can only process one channel containing one impulsive
% % noise peak.
% % 
% % There is an option to plot the data and the 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 Parameters
% % 
% % p is the data array it should be either a column vector or a row vector
% %    p should contain only one impulsive noise peak to be analyzed.
% %    default value is a simulated impoulsive waveform.
% %  
% % Fs is the sampling rate in Hz.  % Default is 100 KHz.
% % 
% % make_plot    % 0 for no plot
% %              % 1 for make the plot
% %              % Default is to not make any plots.
% % 
% % *********************************************************************
% % 
% % Output variables
% % 
% % d is the interpolated D-duration in seconds
% % 
% % ddd is the array of estimated D-durations
% % 
% % ee is the array of total exponential fit error estimates
% % 
% % tau_fit is the estimated decay constants
% % 
% % y2 is the time averaged instantaneous sound pressure ampltitude
% % 
% % y4 is the smoothed time averaged instantaneous amplitude
% % 
% % 
% % *********************************************************************
% % 
% % 
% Example='';
%
% % Example impulsive data with background noise
%
% Fs=50000; fc=200; td=1; tau=0.1; delay=0.1; A1=3; A2=20; risetime=0.0005;
%  
% [p, t]=analytic_impulse(Fs, fc, td, tau, delay, A1, A2, risetime);
%  
% % p               % Pa sound pressure, single channel data array.
%                   % p should have only one impulse.
%  
% Fs=50000;         % Hz sampling rate, 50000 Hz is the default
%  
% make_plot=0;      % 0 for no plot
%                   % 1 for make the plot
%                   % Default is to not make any plots.
%                   % The plot is the absolute value of the data array.
%                   % blue line for the 10 dB threshold. 
%                   % green circle for each point above the threshold.
%                   % red circle for each signel point above the
%                   % threshold. 
%                   % black circle for each first and last point of a
%                   % multiple point series above the threshold.    
%  
% [d, y2, ddd, ee, A, tau_fit]=D_duration(p, Fs, make_plot);
%
% % *********************************************************************
% % 
% % 
% % List of Dependent Subprograms for 
% % D_duration
% % 
% % 
% % Program Name   Author   FEX ID#
% % 1) analytic_impulse		
% % 2) findextrema		Schuberth Schuberth		3586	
% % 3) geomean2		
% % 4) LMS_trim		
% % 5) LMTSregor		
% % 6) rand_int		
% % 7) sd_round		
% % 8) sub_mean		
% % 9) threshold_bin_peaks	
% % 
% % 
% % *********************************************************************
% % 
% % D_duration was written by Edward L. Zechmann  
% % 
% %      date 11 December 2007
% % 
% %  modified 15 December 2007  added more comments
% % 
% %  modified 15 December 2007  added LMTSregor for robust least trimmed 
% %                             squares regression made the program faster
% %                             updated comments
% % 
% %  modified 15 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, B_Duration, C_Duration
% % 


if nargin < 1
    make_plot=1;
    fc=200; td=1; tau=0.1; delay=0.1; A1=3; A2=20;
    [p, t]=analytic_impulse(Fs, fc, td, tau, delay, A1, A2);
    warning('Not enough input arguments: Running Example Code');
end

if nargin < 2
    Fs=100000;
end

if nargin < 3
    make_plot=0;
end

d=[];
ddd=[];
ee=[];
tau_fit=[];
y2=[];
y4=[];





num_data=length(p);

% y2 is the running average of the instantaneous sound pressure amplitude.
[buf, y2]=sub_mean(abs(hilbert(p)), Fs, 2000);

% y4 is a smoother fit of the running average
[buf2, y4]=sub_mean(y2, Fs, 100);

% maxy2 is the maximum value of y2
[maxy2 maxy2_index]=max(y2);

[buf2, y4]=sub_mean(y2, Fs, 100);
[maxy4 maxy4_index]=max(y4);
[may4, miy4]=findextrema(y4);

% make sure that indices are integers and are within the range of the 
% array
may4=round(may4);
may4(may4 < 1)=1;
may4(may4 > num_data)=num_data;
maxy4_ix=may4(find(may4 > maxy4_index));

if ~isempty(maxy4_ix)
    maxy4_vals=y4(maxy4_ix);

    num_max=length(maxy4_ix);

    if num_max > 120
        maxy4_vals_fit=maxy4_vals(1:10:num_max);
        maxy4_ix_fit=maxy4_ix(1:10:num_max);
    else
        maxy4_vals_fit=maxy4_vals(1:num_max);
        maxy4_ix_fit=maxy4_ix(1:num_max);
    end
    
    yint=geomean2(maxy4_vals_fit);
    % initialize intercept to the last maxima after the global peak

    xint=maxy4_ix_fit(end);
    % find an earlier crossing of the mean local maxima peak

    x_data_int=maxy2_index:xint;
    y_data_int=y4(x_data_int);

    x_int_ix=x_data_int(find(y_data_int < yint, 1 ));

    if (x_int_ix-maxy2_index)/Fs < 0.001
        x_int_ix=xint;
    end
    
    if isempty(x_int_ix)
        x_int_ix=xint;
    end
else
    x_int_ix=num_data;
    xint=num_data;
    x_data_int=maxy2_index:xint;
    y_data_int=y4(x_data_int);
    yint=y_data_int;
end

% Fit y2 to an exponential curve at five different threshold values
% then select the estimated D-duration with the least error

% array of thresholds dB
% the expoenential fit will be done at 4 dB down ,then 8 dB down, etc...
threshold_a=[ 4, 8, 10, 15, 20];

% number of threshold values
num_thres=length(threshold_a);

% initialize the error array
Error=zeros(num_thres, 1);
dd=zeros(num_thres, 1);
tau_fit=zeros(num_thres, 1);

for e2=1:num_thres;

    thres=10^(-threshold_a(e2)/20);

    maxy2_thresh=maxy2*thres;
    
    % Allow all data points
    % fitpts=1/maxy2*y2(maxy2_index-1+find( y2(maxy2_index:length(y2)) >= maxy2_thresh));
    % Allow all data points up to the temporal threshold set by the
    % background level
    fitpts_ix=maxy2_index-1+find( y2(maxy2_index:x_int_ix) >= maxy2_thresh);
    fitpts=1/maxy2*y2(fitpts_ix);
    
    % get the log10 of the data
    fitpts2=log10(fitpts);
    
    if length(fitpts2) >= 10
        
        y3=y2(maxy2_index:x_int_ix);
        
        [maix]=threshold_bin_peaks(y3, maxy2_thresh, 10, 1000);
        
        % must have at least 1 data point to fit
        if isempty(maix)
            maix=1:length(fitpts2);
        end
        
        % Append the maxima but start at 1  
        % for the exponential fit
        if maix(1) > 1
            maix=[1 maix];
        end
        
        fitpts3=log10(1/maxy2*y2(maxy2_index-1+round(maix)));
        
        ma_val=(-fitpts3(1)+fitpts3)';
        
        if length(ma_val) >= 2
            % Method 1) Best fit line that goes through the global peak
            
            % best fit line that goes through zero
            % LMTSregor stands for
            % Least Median Trimmed of Squares regression through the origin
            [LMSout,blms,Rsq]=LMTSregor(ma_val, [-1+maix]', 1000, 100000 );
            % LMSout is the LMTS estimated values vector.
            % blms is the LMTS slopes vector.
            % Rsq is the squared error term
            Error(e2)=Rsq;
            tau_fit(e2)=log10(exp(1))/(-blms*Fs);
            
            % calculate the best fit exponential curve
            p3=maxy2.*(10.^polyval([blms fitpts2(1)], [0:(-1+length(fitpts2))] ));
            
        else
            maix=1:length(fitpts2);
            [pc,S]=polyfit(1:length(fitpts2),fitpts2,1);
            blms=pc(1);
            Error(e2)=S.normr.^2;
            % calculate the best fit exponential curve
            p3=maxy2.*(10.^polyval(blms,fitpts2 ));
        end
        
        tau_fit(e2)=log10(exp(1))/(-blms*Fs);
        
    elseif length(fitpts2) >= 3
        maix=1:length(fitpts2);
        [pc,S]=polyfit(1:length(fitpts2),fitpts2,1);
        blms=pc(1);
        Error(e2)=S.normr.^2;
        % calculate the best fit exponential curve
        p3=maxy2.*(10.^polyval(blms,fitpts2 ));
        tau_fit(e2)=log10(exp(1))/(-blms*Fs);
    else
        maix=1;
        p3=1;
    end
    
    % Determine the index of the intersection of the 10 dB threshold with
    % the best-fit exponential curve
    
    if length(p3) > 1

        thres=10^(-10/20);
        maxy2_thres=p3(1)*thres;
        flag1=0;

        if (maxy2_index-1) > 1 && ((maxy2_index-1) <= length(p3))
            e1=(maxy2_index-1);
        else
            e1=1;
        end

        while (flag1 == 0) && (e1 < length(p3))
            e1 = e1+1;
            if p3(e1) <= maxy2_thres
                flag1=1;
            end
        end

        % interpolate to find better end time
        if abs(p3(e1)-p3(e1-1)) >= 10^-12
            dt2 = (maxy2_thres-p3(e1-1))/(p3(e1)-p3(e1-1))+e1-1;
        else
            dt2=e1;
        end
    else
        dt2=2;
        Error(e2)=1;
    end
    
    % D-duration in indices
    % dt1 is the index of the peak which is shifted to 1.
    % dt2 is the index of the p3 curve crossing the threshold.
    dt1=1;
    dd(e2)=(dt2-dt1);

    if isequal(make_plot, 1)
        figure(100+e2);
        % plot the running average of the instantaneous amplitude
        plot(1/Fs*(1:num_data), y2, '-b', 'linewidth', 0.5);
        hold on;
        npts=length(p3);
        
        % Plot the Level threshold
        plot(1/Fs*(maxy2_index-1+[1 npts]), maxy2_thresh*[1 1], '--kx', 'linewidth', 2);
        
        % Plot the curve fit data points
        plot(1/Fs*round(maxy2_index-1+round(maix)), y2(maxy2_index-1+round(maix)), 'LineStyle', 'none', 'Marker', 'o', 'MarkerEdgeColor', 'b', 'MarkerSize', 7);
        
        % Plot the exponential curve fit
        plot(1/Fs*(maxy2_index-1+(1:npts)), p3(1:npts), '-k', 'linewidth', 2);
        
        % Plot smoothed time average
        plot(1/Fs*(1:num_data), y4, '-r', 'linewidth', 2 );
        
        % Plot smoothed time average end points
        plot(1/Fs*[1 maxy2_index xint num_data], y4([1 maxy2_index xint num_data]), 'or', 'linewidth', 3);
        
        % Plot data for determining temporal threshold
        plot(1/Fs*x_data_int, y_data_int, 'm', 'linewidth', 2);
        
        plot(1/Fs*x_int_ix, y4(x_int_ix), 'om', 'linewidth', 3);
        xlabel('Time seconds', 'Fontsize', 18);
        ylabel('Sound Pressure Amplitude (Pa)', 'Fontsize', 18);
        title(['D-Duration ', num2str(sd_round(1/Fs*dd(e2), 2)), ' s, tau ', num2str(sd_round(tau_fit(e2), 2)), ' s, Background Amplitude ', num2str(sd_round(yint, 2 )), ' Pa'], 'Fontsize', 18);
        set(gca, 'Fontsize', 16);
        legend('Instantaneous Amplitude', 'Level Threshold', 'Fit Data Points', 'Exponential Fit', 'Smoothed Average', 'Smoothed Average Endpoints', 'Data for Temporal Threshold', 'Temporal Threshold');
    end
    
end

% Calculate the total error
% error times the duration array
ee=dd.*Error; 

% Make sure that the total duration was positive
index5=find(ee>0.00000001);

% find the index of the least total error
[minval minin]=min(ee(index5));

% Get the durations with positive total duration
% Convert the estimated d-duration to seconds
ddd=1/Fs*dd(index5);

% Select the D-duration with the least error*duration
d=ddd(minin);