calc_impuls_threshold_and_index.m

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

function [indices, threshold_a ]=calc_impuls_threshold_and_index(gp, percent1, exhaust_cycle, num_pk_intrls, make_plot)
% % calc_impuls_threshold_and_index:  Calculates the threshold and finds the indices above the threshold 
% % 
% % Syntax:
% % 
% % [indices, threshold_a ]=calc_impuls_threshold_and_index(gp, percent1, exhaust_cycle);
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % Description  
% % 
% % This program calculates the threshold for determining which data points
% % are impulses and which are background noise.  
% % 
% % After determining the threshold the program finds which points are
% % above the threshold are returns a list of these points in the matrix
% % "indices".  
% % 
% % The input and output variables are described below.
% % 
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % Input Variables
% % 
% % gp=randn(num_channels, numbins);
% %                     % gp is the data array which is optimized for
% %                     % finding peaks
% %
% % percent1=[];        % Is the percentage of the peak levels of the  
% %                     % time record after signal processing to remove
% %                     % continuous noise and sharpen the peak levels.
% %                     % The default is percent1=[];    
% %                     % 
% %                     % If percent1 is empty, < 0, or > 100 then 
% %                     % it is set to [] and not considered during
% %                     % processing;   
% %                     % 
% %                     % percent1=75; is a typical value
% % 
% % exhaust_cycle=1;    % 1 use the threshold method approppriate for 
% %                     % finding exhaust noise type impulses.   
% %                     % 
% %                     % Otherwise the method for impact type impulsive
% %                     % noises is used.
% %                     % 
% %                     % The default value is exhaust_cycle=0;
% % 
% % num_pk_intrls
% % 
% % make_plot
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % Output Variables
% % 
% % indices is a cell array with teh same number of channels as gp.
% %          indices is the set of indices of gp which is above the 
% %          threshold.  
% % 
% % threshold_a is two dimensional array of size [num_channels num_bins]  
% %          each element is the threshold for that particular bin.  
% % 
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% Example='1';
%
% Fs_SP=100000; fc=1000; td=1; tau=0.01; delay=0.1; A1=8; A2=21;
% [SP1, t]=analytic_impulse(Fs_SP, fc, td, tau, delay, A1, A2);
% Fs_SP=100000; fc=1000; td=1; tau=0.1; delay=0.1; A1=2; A2=22;
% [SP2, t]=analytic_impulse(Fs_SP, fc, td, tau, delay, A1, A2);
% Fs_SP=100000; fc=1000; td=2; tau=0.4; delay=0.1; A1=2; A2=17;
% [SP3, t]=analytic_impulse(Fs_SP, fc, td, tau, delay, A1, A2);
% Fs_SP=100000; fc=2000; td=1; tau=0.1; delay=0.1; A1=8; A2=15;
% [SP4, t]=analytic_impulse(Fs_SP, fc, td, tau, delay, A1, A2);
% SP=[[0.9*SP4 SP1 0.3*SP3 SP2 0.5*SP4 0.2*SP1 SP3 5*SP2];...
%     [SP1 0.7*SP4  10*SP2 SP3 0.4*SP2 0.5*SP3 SP1 1*SP4];...
%     [SP3 0.2*SP1  0.5*SP4 2*SP2 0.4*SP2 0.7*SP3 0.3*SP1 10*SP4];...
%     [SP1 0.7*SP2  0.3*SP3 SP4 0.4*SP1 10*SP2 SP3 1*SP4]];
% 
% bin_size=Fs_SP/100; 
% 
% [gp]=peak_threshhold_function(SP, Fs_SP, bin_size, exhaust_cycle);
% 
% percent1=[];
% exhaust_cycle=1;
% num_pk_intrls=50;
% make_plot=1;
% 
% [indices, threshold_a ]=calc_impuls_threshold_and_index(gp, percent1, exhaust_cycle, num_pk_intrls, make_plot);
% 
% 
% 
% Example='2';
% % Use gp from the above code 
% 
% By settign the exhaust cycle to 1 more peaks are found.
% Exhaust cycle conmtrols part of the peak identification function and
% threshold. 
% 
% percent1=[];
% exhaust_cycle=0;
% num_pk_intrls=50;
% make_plot=1;
% [indices, threshold_a ]=calc_impuls_threshold_and_index(gp, percent1, exhaust_cycle, num_pk_intrls, make_plot);
% 
% 
% 
% Example='3';
% % Use gp from the above code 
% 
% By setting percent1 to 90%, instead of empty many teh same number of peaks are found.
%
% percent1=90;
% exhaust_cycle=0;
% num_pk_intrls=50;
% make_plot=1;
% [indices, threshold_a ]=calc_impuls_threshold_and_index(gp, percent1, exhaust_cycle, num_pk_intrls, make_plot); 
%
%
% Example='4';
% % Use gp from the above code 
% 
% By setting percent1 to 50%, more peaks are found.
%
% percent1=90;
% exhaust_cycle=0;
% num_pk_intrls=50;
% make_plot=1;
% [indices, threshold_a ]=calc_impuls_threshold_and_index(gp, percent1, exhaust_cycle, num_pk_intrls, make_plot); 
%
%
% Example='5';
% % Use gp from the above code 
% 
% By setting percent1 to 50%, and using the exhaust_cycle a few more peaks 
% are found.
%
% percent1=50;
% exhaust_cycle=1;
% num_pk_intrls=50;
% make_plot=1;
% [indices, threshold_a ]=calc_impuls_threshold_and_index(gp, percent1, exhaust_cycle, num_pk_intrls, make_plot); 
%
%
% ***********************************************************
% % 
% % 
% % List of Dependent Subprograms for 
% % calc_impuls_threshold_and_index
% % 
% % FEX ID# is the File ID on the Matlab Central File Exchange
% % 
% % Program Name   Author   FEX ID#
% % 1) allstats		Duane Hanselman		NA	
% % 2) func_threshold		Jing Tian		10462	
% % 3) kurtosis2		
% % 4) percentiles		
% % 5) wsmooth		Damien Garcia	
% % 
% ***********************************************************
% % 
% % calc_impuls_threshold_and_index is written by Edward Zechmann 
% % 
% %     date 21 August      2008
% % 
% % modified  3 September   2008    Updated Comemnts
% % 
% % modified 10 September   2008    Added exhaust_cycle finding option
% %                                 Updated Comments
% % 
% % modified 11 November    2008    Fixed a bug with the exhaust noise
% %                                 feature and with bin size.  
% %                                 Addeded more examples.  
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % Please Feel Free to Modify This Program
% % 
% % See Also: wsmooth, percentiles, func_threshold
% % 


if (nargin < 1 || isempty(gp)) || ~isnumeric(gp) 
    gp=rand(1, 2000);
end

[m1, num_bins]=size(gp);

flag2=1;
% set the percentile for finding impulsive peaks
if (nargin < 2 || isempty(percent1)) || ((~isnumeric(percent1) || logical(percent1 > 100)) || logical(percent1 <= 0))
    percent1=[];
    flag2=0;
end

if nargin < 3 || isempty(exhaust_cycle) || ~isnumeric(exhaust_cycle)
    exhaust_cycle=0;
end

if nargin < 4 || isempty(num_pk_intrls) || ~isnumeric(num_pk_intrls)
    num_pk_intrls=1;
end

if nargin < 5 || isempty(make_plot) || ~isnumeric(make_plot)
    make_plot=0;
end


% determine threshold for each channel

%
% ***********************************************************
%
% Philoshophy of Calculating the Threshold for Detecting Impulsive
% Peaks
%
% I)  Absolute Threshold Requirement
%
% An absolute threshold requirement is necessary to avoid
% mistaking a pure tone or a pure random noise
% as an impulsive noise.
%
% This is accomplished with the code
% threshold=max(0.0004, threshold);
%
%
% II)  Relative Threshold Requirement
%
% A relative threshold requirement is used to distinguish
% between the background noise and impulsive noise.
%
% 1)  global_max
%     The global maximum peak for the channel is found.
%     Only the absolute threshold can be greater than the global
%     maximum peak value.  This ensures that at least one peak
%     will be found if sar=1.
%
% 2 ) First the array of the peaks is broken up into data_bins
%     containing all of the peaks within num_pts_impulse.
%     Steps 3 through 9 are repeated for each bin of 100 peaks.
%
% 3)  The maxima of the array of peaks are found.
%     The indicies are assigned to the variable ma.
%
% 5)  max
%     Calculate the maximum value of the array of the
%     peaks.  This value is called max.
%
% 8)  Calculate of the Threshold Value
%
% relative threshold value is based largely on the scaled kurtosis value
%
% relative threshold value = min([max/kurtosis*2/3)+mean+std,global_max]);
%
% 9) Calculate the absolute threshold requirement
%
%    threshold=max(0.0004, threshold);
%

[m1, num_bins]=size(gp);

pts_per_pk_intrl=floor(num_bins/num_pk_intrls);

threshold_a=zeros(m1,num_pk_intrls);

bbmaxgp=max(gp, [], 2);

indices=cell(m1,1);


for e1=1:m1;

    buf=[];
    bb=wsmooth(gp(e1, :), 2);
    if isequal(flag2, 1)
        [percts]=percentiles(gp(e1, :), 2, percent1);
    else
        percts=bbmaxgp(e1);
    end
    
    calc_thres3 = func_threshold(bb);
    
    for e2=1:num_pk_intrls;

        % Select the peaks that fall within the num_pts_impulse
        % depends on the length of the array and the number of peaks
        % available for analysis
        pts=pts_per_pk_intrl*(e2-1)+(1:pts_per_pk_intrl);

        if max(pts) > num_bins || isequal(e2, num_pk_intrls)
            pts=min(pts):num_bins;
        end
        
        calc_thres2=max(max(gp(pts)),bbmaxgp(e1)/num_pk_intrls)/(kurtosis2(bb(pts))/3)+percentiles(gp(e1, pts), 2, 75);
        
        if isequal(exhaust_cycle, 1)
            calc_thres=0.5*calc_thres3;
        else
            calc_thres=calc_thres3;
        end
        
        threshold=min([percts, calc_thres, calc_thres2, bbmaxgp(e1)]);

        % an absolute threshold requirement is necessary to avoid
        % mistaking a pure tone or a pure random noise
        % as an impulsive noise.
        % 
        % The nunmber 0.0004 is based upon simulations with sin waves 
        % and noise generated with randn.  
        
        threshold_a(e1, e2)=max([0.0004, threshold]);

        buf=[buf (pts_per_pk_intrl*(e2-1)+(find(gp(e1, pts) >= threshold)))];
        
    end
 
    indices{e1,1}=buf;
    
end


sh=0.0;
sv=0.0;
ml=0.14;
mr=0.1;
mt=0.08;
mb=0.12;
n1=length(gp);

if isequal(make_plot,1)

    % initialize the figure
    figure(10000);
    delete(10000);
    figure(10000);
    ylim1=0.1*max(max(max(ceil(11*gp))));
    xlim1=n1;
    h2=[];

    tot_num_samples=m1*n1;
    flag_rs=0;
    if tot_num_samples > 1000000
        flag_rs=1;
    end

    t_SP_rs=[];
    SP_rs=[];
    
    
    for e1=1:m1;

        subaxis(m1, 1, e1, 'sh', sh, 'sv', sv , 'pl', 0, 'pr', 0, 'pt', 0, 'pb', 0, 'ml', ml, 'mr', mr, 'mt', mt, 'mb', mb);
        t_SP=(0:(length(gp(e1, :))-1));
        
        if flag_rs == 1
            [t_SP_rs, SP_rs]=resample_plot(t_SP, gp(e1, :));
            plot(t_SP_rs, SP_rs);
        else
            plot(t_SP,gp(e1, :));
        end
        
        if isequal(e1, 1)
            title({'', 'Impulsive Noise Detection Routine', 'Processed Signal and Peak Threshold Line'}, 'Fontsize', 18);
            ylabel('Amplitude', 'Fontsize', 18);
        end
        
        h2=[h2 gca];
        set(gca, 'Fontsize', 12, 'box', 'off' );
        
        if e1 ~= m1
            set(gca, 'xtick', [], 'XTickLabel', '');
        end

        text( 'Position', [num_bins*0.95, 0.95*ylim1], 'String', ['Microphone ', num2str(e1)], 'Fontsize', 20, 'Color', [0 0 0], 'HorizontalAlignment', 'right', 'VerticalAlignment', 'top', 'BackgroundColor', [1 1 1] );

        hold on;
        
        
        for e2=1:num_pk_intrls;
            
            % Select the peaks that fall within the num_pts_per_pk_intrl
            % depends on the length of the array and the number of peaks
            % available for analysis
            pts=pts_per_pk_intrl*(e2-1)+(1:pts_per_pk_intrl);
            
            if max(pts) > num_bins || isequal(e2, num_pk_intrls)
                pts=min(pts):num_bins;
            end
            
            % an absolute threshold requirement is necessary to avoid
            % mistaking a pure tone or a pure random noise
            % as an impulsive noise.
            threshold=threshold_a(e1, e2);
            plot([min(pts) max(pts)], threshold*[1 1], 'r');
        end
    
        set(gca, 'Fontsize', 16);
        ylim([0 ylim1]);
        xlim([0 xlim1]);
        [ytick_m, YTickLabel1, ytick_good, ytick_new, yticklabel_new]=fix_YTick(0);
        
    end

    xlabel('Data Points', 'Fontsize', 18);
    %legend('Processed Signal Amplitude', 'Peak Detection Threshold');

    if length(h2) >= 1
        psuedo_box(h2);
    end

end