c_duration.m

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

function [c, c2, error_cond]=C_duration(p, Fs, make_plot)
% % C_duration: Calculates the C-duration for impulsive noise 
% % 
% % Syntax:  [c, c2, error_cond]=C_duration(p, Fs, make_plot);
% %
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% % Description 
% % 
% % This program calculates the C-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.  
% %  
% % make_plot=1;    plots the time records and places a circle at each 
% %                 chosen peak. Otherwise no plots are made. 
% %                 default value is 0.
% %  
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % Output variables
% % 
% % c is the interpolated C-duration in seconds
% % 
% % c2 is the rough estimate of the C-duration by counting points without 
% % interpolation
% % 
% % error_cond is 1 if the absolute value of the difference betweeen 
% % c and c2 is greater than the sum of the multiple points and the single
% % points above the threshold.  
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% 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 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.   
%
% [c, c2, error_cond]=C_duration(p, make_plot, Fs);
% 
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% % C_duration.m was originally developed by Chucri Kardous.  
% % 
% % C_duration was modified by Edward L. Zechmann  
% %
% %      date 11 December  2007
% %
% %  modified 15 December  2007  Added more 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 sampling rate is not specified, then 
% assume a data acquisition rate of 50 KHz;
if (nargin < 2 || isempty(Fs)) || ~isnumeric(Fs)
    Fs=100000; 
end

% The default is to not make any plots
if (nargin < 3 || isempty(make_plot)) || ~isnumeric(make_plot)
    make_plot=0; 
end


p2=abs(p);

% C-duration
thres=10^(-10/20);
maxp_thres=max(p2)*thres;
ctime=find(p2 > maxp_thres);

c=0;    
ct1=0; 
ct2=0;
c2=length(ctime);

dct=diff(ctime);
    
flag1=0;
mpa=[];
spa=[];

if isempty(dct)
    % C-duration for a single very large impulse
    thres=10^(-10/20);
    [maxp maxp_index]=max(p2);
    maxp_thres=maxp*thres;
    x_thres_index=find(p2 > maxp_thres);
    
    % find start time for C-duration, which is the first 10 dB 
    % threshold-crossing before the peak pressure
    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
    flag1=0;
    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
        ct1 = (maxp_thres-p2(e1))/(p2(e1+1)-p2(e1))+e1;
    else
        ct1=e1;
    end
    
    % find end time for C-duration, i.e. the last 10 dB threshold-crossing
    % after peak
    if (x_thres_index(end)-1) > 1 && ((x_thres_index(end)-1) <= length(p2))
        e1=(x_thres_index(end)-1);
    else
        e1=2;
    end
    
    flag1=0;
    while (flag1 == 0) && (e1 < length(p2))
        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
        ct2 = (maxp_thres-p2(e1-1))/(p2(e1)-p2(e1-1))+e1-1;
    else
        ct2=e1;
    end
    
    % C-duration in indices
    c = (ct2-ct1);  
    c2=c;        
    error_cond=0;
        
else
    % mpa multiple points above threshold
    % spa single point above threshold
    
    mpa=zeros(c2, 1);
    spa=zeros(c2, 1);
    
    mpc=0; % multiple point counter
    spc=0; % single point counter
    
    for e1=1:(length(dct));
        
        if dct(e1) > 1
            if flag1 == 1
                mpc=mpc+1;
                mpa(mpc)=ctime(e1);
            else
                spc=spc+1;
                spa(spc)=ctime(e1);
            end 
            flag1=0;
        else
            if flag1 == 0
                mpc=mpc+1;
                mpa(mpc)=ctime(e1);
            end 
            flag1=1;   
        end
    end
    
    mpa=mpa(1:mpc);
    spa=spa(1:spc);
    
    t_array=[];
    
    % Every pair of multiple points above threshold 
    % must have an end point 
    % if no last end point, then the last 
    % point becomes the last end point
    if mod(length(mpa), 2) > 0
        mpc=mpc+1;
        mpa(mpc)=ctime(end);
    end
    
    nmpa=floor(length(mpa)/2);
    t_array=zeros(nmpa+spc, 1);
    
    for e1=1:nmpa;
        
        % get height of previous data point
        if (mpa(2*e1-1)-1) > 1
            pt1=p2(mpa(2*e1-1)-1);
        else
            pt1=p2(1);
        end
        % get height of data point
        pt2=p2(mpa(2*e1-1));
       
        % get height of next data point
        if (mpa(2*e1)+1) < length(p2)
            pt3=p2(mpa(2*e1)+1);
        else
            pt3=p2(end);
        end
        
        % get height of data point
        pt22=p2(mpa(2*e1));
        
        % interpolate for data point before crossing threshold
        if abs(pt2-pt1) >= 10^-12
            t3L=1*(pt2-maxp_thres)/(pt2-pt1)+0;
        else
            t3L=0;
        end
        
        % interpolate for data point after crossing threshold
        if abs(pt22-pt3) >= 10^-12
            t3U=1*(pt22-maxp_thres)/(pt22-pt3)+0;
        else
            t3U=0;
        end
        
        t_array(e1)=t3U+t3L+mpa(2*e1)-mpa(2*e1-1);
        
    end
    
    for e1=1:spc;
        
        % get height of previous data point
        if (spa(e1)-1) > 1
            pt1=p2(spa(e1)-1);
        else
            pt1=p2(1);
        end
        
        % get height of next data point
        if (spa(e1)+1) < length(p2)
            pt3=p2(spa(e1)+1);
        else
            pt3=p2(end);
        end
        
        % get height of data point
        pt2=p2(spa(e1));
        
        % interpolate for data point before crossing threshold
        if abs(pt2-pt1) >= 10^-12
            t3L=1*(pt2-maxp_thres)/(pt2-pt1)+0; 
        else
            t3L=0;
        end

        % interpolate for data point after crossing threshold
        if abs(pt2-pt3) >= 10^-12
            t3U=1*(pt2-maxp_thres)/(pt2-pt3)+0; 
        else
            t3U=0;
        end
        
        t_array(e1+nmpa)=t3U+t3L;
        
    end
    
    % calculate the interpolated sum
    c=sum(t_array);
    % calculate the rough estimated sum
    c2=length(ctime);
    
    % Calculate if the maximum error is satisfied
    error_cond=0;
    max_error=spc + nmpa;
    if abs(c2 - c) > max_error
        error_cond=1;
    end
    
    if make_plot == 1
        close all;
    
        figure(2);
        % plot the absolute value of the data array
        plot(p2);
        hold on;
        
        % draw the threshold line
        plot( [1 length(p2)], maxp_thres*[1 1]);
        
        % put a green circle for each point above the threshold
        for e1=1:length(ctime); 
            plot( ctime(e1), p2(ctime(e1)), 'og', 'linestyle', 'none', 'markersize', 7);
        end
        
        % put a red circle for each single point above the threshold
        for e1=1:length(spa); 
            plot( spa(e1), p2(spa(e1)), 'or', 'linestyle', 'none', 'markersize', 7);
        end
        
        % put a black circle for the beginning and end point of each series of multiple points above the threshold
        for e1=1:length(mpa); 
            plot( mpa(e1), p2(mpa(e1)), 'ok', 'linestyle', 'none', 'markersize', 7);
        end
    
    end
    
end


% convert c and c2 to seconds
c=1/Fs*c;
c2=1/Fs*c2;

if isempty(c)
    c=-1;
end

if isempty(c2)
    c2=-1;
end