localpeaks.m

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

if pts_per_bin < 1
    pts_per_bin=1;
end

if pts_per_bin > n1
    pts_per_bin=n1;
end

max_num_bin=floor(n1/pts_per_bin);

% Calculate the number of bins per peak interval 
% 
% Number of points per peak interval divided by
% the number of points in each bin.
bins_per_pk_intvl=floor(num_pts_per_pk_intrl/pts_per_bin);
if bins_per_pk_intvl < 1
    bins_per_pk_intvl=1;
end


% In complex noise, continuous noise and impulsive noise are combined.
%
% Now process SP2 to minimize or remove the continuous noise 
% component of the signal.
%
% After removing the continuous noise, maximize the peakedness of the 
% signal to expose the underlying impulsive peaks.
%
% Apply a threshold to the peakedness of the exposed impulses to select
% impulses that are significantly greater than the remaining noise.  

[gp]=peak_threshhold_function(SP, Fs, pts_per_bin, exhaust_cycle);


[mm1, num_bins]=size(gp);

% Calculate the number of bins per peak interval
if bins_per_pk_intvl > num_bins
    bins_per_pk_intvl=num_bins;
end

% Calculate the number of peak intervals
num_pk_intrl=max([1, floor(num_bins/bins_per_pk_intvl)]);


% find the maximum value of gp for each channel. 
% If sar==1 then at least one impulsive peak must always be returned.
% 
% From the comments on the Philosophoy of Detecting Impulsive Noises
% 
% 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.
%
%     Since at least one peak must be found when sar=1;
%
% Plot the array of peaks and the impulsive peak threshold. 
% Peaks above the threshold line are impulsive peaks. 
%
% Being above the impulsive peak threshold is not sufficient for a peak 
% to be analyzed.  
% 
% To be analyzed, either sar=1 or the corresponding sound pressure peak 
% must be above min_peak.   
% 

[indices, threshold_a ]=calc_impuls_threshold_and_index(gp, percent1, exhaust_cycle, num_pk_intrl);

[buf, num_pk_intrl]=size(threshold_a);
pts_per_pk_intrl=floor(n1/num_pk_intrl);

if isequal(make_plot,1)

    % initialize the figure
    figure(10000);
    delete(10000);
    figure(10000);
    ylim1=0.1*max(max(max(ceil(11*gp))));
    xlim1=1/Fs*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=pts_per_pk_intrl/Fs*(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', [pts_per_pk_intrl*num_bins*0.95/Fs, 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_intrl;
            
            % The peaks that are above the threshold are selected
            pts=pts_per_pk_intrl*(e2-1)+[e2 pts_per_pk_intrl];
            
            % 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(pts./Fs, 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('Time seconds', 'Fontsize', 18);
    %legend('Processed Signal Amplitude', 'Peak Detection Threshold');

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

end
    


SP_peaka=cell(m1,1);
index_peaka=cell(m1,1);
% The cell array index_peaka contains all peaks above the threshold
% The cell array SP_peaka continas the amplitude of those peaks

% set the radius for finding a peak from the center of a bin.
radius=0.5;

for e1=1:m1;
    [SP_peak, index_peak]=peak_index(abs(SP(e1, :)), indices{e1,1}, pts_per_bin, radius*num_pts_per_pk_intrl);
    SP_peaka{e1, 1}=SP_peak;
    index_peaka{e1, 1}=index_peak;
end

% Enforce the min_peak requirement
% index_peaka now only includes peaks above the min_peak amplitude (Pa)
if ~isequal(sar, 1)
    for e1=1:m1;
        index1=find(abs(SP_peaka{e1, 1}) > min_peak);
        SP_peaka{e1, 1}=SP_peaka{e1, 1}(index1);
        index_peaka{e1, 1}=index_peaka{e1, 1}(index1);
    end
end

%SP_local_max=[];
%SP_local_max2={};

low_bin=floor(0.25*num_pts_per_pk_intrl);
high_bin=floor(num_pts_per_pk_intrl-low_bin);

index1=0;
index2=0;

% Align_peaks is a channel number that the peaks are aligned to.  
% 
% align_chans is an array of channels to be aligned to the channel 
% designated by Align_peaks.
% 
if ~isempty(Align_peaks)
    align_chans=setdiff(1:m1, Align_peaks);
else
    align_chans=[];
end

indices2=cell(m1,1);
SP_local_max=[];
SP_local_max2=cell(m1,1);


for e1=1:m1;

    if isempty(Align_peaks)
        e2=e1;
    else
        if isequal(e1,1)
            e2=Align_peaks;
        else
            e2=align_chans(e1-1);
        end
    end
    
    % whether or not to find all of the peaks
    if isempty(Align_peaks) || isequal(e1,1)
        
        % Get the index of the maximum peak 
        % first get teh indices of all of the peaks for the alignment
        % channel or channel 1 
        ix_peak=index_peaka{e2, 1};
        SP_peak=abs(SP_peaka{e2, 1});
        
        % initialize the vector of major peak indices
        major_peak_ix=[];
        
        % Starting with the highest amplitude peak, one by one identify  
        % each peak and discard all other peaks within the radius
        % of the peak.  
        %
        % Each peak has a set number of points before and 
        % after the peak which belong to that peak and only that peak.   
        % 
        while length(ix_peak) >= 1
            
            % get the index of the peak with the greatest amplitude.
            [val, ix]=max(SP_peak);
            cix=ix_peak(ix);

            % Get all of the indices within the umbrella of the maximum peak
            index1=cix-low_bin;
            index2=cix+high_bin;

            if index1 < 1
                index1=1;
                index2=index1+low_bin+high_bin;
            end

            if index2 > n1
                index2=n1;
            end

            % minor_peaks_ix are the indices within the umbrella of the
            % minor local peaks 
            gti1=find(ix_peak >= index1);
            lti2=find(ix_peak <= index2);
            minor_peaks_ix=intersect( gti1, lti2);
            not_m_peaks=setdiff(1:length(ix_peak), minor_peaks_ix);
            ix_peak=ix_peak(not_m_peaks);
            SP_peak=SP_peak(not_m_peaks);

            % Make sure each peak has a zero crossing after the peak
            % find the zero crossing after the peak
            flag1=0;
            e4=cix;
            if SP(cix) < 0
                SP2=-SP;
            else
                SP2=SP;
            end
            
            while (flag1 == 0) && (e4 < n1)
                e4 = e4+1;
                if SP2(e4) <= 0
                    flag1=1;
                end
            end

            % check for zero crossing after max peak
            % interpolate the zero crossing
            % calculate the slope
            %if abs(SP2(e4)-SP2(e4-1)) >= 10^(-12)
            %    at2 = (0-SP2(e4-1))/(SP2(e4)-SP2(e4-1))+e4-1;
            %    m=SP2(e4)-SP2(e4-1);
            %else
            %    m=0;
            %    at2=e4;
            %end
            %
            % There must be a zero crossing with a negative
            % slope for the peak to be accepted.
            %if (at2 <= n1) && m <= 0
                major_peak_ix=[major_peak_ix cix];
            %end
        end
        
        % sort the indices 
        [sorted_mp_ix]=sort(major_peak_ix);
        
        % sorted_mp_ix is the sorted major peak indices
        % sorted_mp_ix will be apended to the cell array (indices2)
        % located after the conclusion of this conditional statemnet 
        
        % Prepare to append the peak amplitudes to the peak amplitude
        % cell array 
        SP_local_max=SP(e2, sorted_mp_ix);
        
    else
        % This code is for the case where teh peaks in channel e2 must
        % be aligned to the peaks in channel Align_peaks.  
        % 
        % Determine the best Alignment of peak indices for the
        % microphone channel (e2) using the peak indices that 
        % were found for the Align_peaks channel.  
        % 
        % All of the impulsive peak indices for channel e2 were
        % previously found.  
        % 
        % Starting with the highest peak from channel Align_peaks the 
        % corresponding peak in channel e2 will be determined.  
        % 
        % The radius of points in channel e2 will be determined by pts_per_bin
        % then the the number of points will be reduced to 
        % PNTS=peak_alignment_tol*num_pts_per_pk_intrl.  The highest peak within the
        % PNTS of the Align_peaks index.  
        % 
        % points of the indices from the aling peaks channel.  
        % 
        % for this bin for this microphone channel
        % using the bins for the Align_peaks microphone
        % 
        indices1=ceil(indices2{Align_peaks}/pts_per_bin);
        indices1(indices1 < 1)=1;
        indices1(indices1 > max_num_bin)=max_num_bin;
        
        % Choose data points for the center of each local peak impulse
        % The peak_index program uses indices rounded to the pts_per_bin,
        % determines the appropriate bin, then finds the peak index.
        
        [SP_maxa, index_maxa]=peak_index(SP(e2, :), indices1, pts_per_bin, peak_alignment_tol*num_pts_per_pk_intrl);
        
        % sort the indices 
        [sorted_mp_ix]=sort(index_maxa);
        SP_local_max=SP(e2, sorted_mp_ix);
    end

    SP_local_max2{e2, 1}=SP_local_max;
    indices2{e2, 1}=shiftdim(sorted_mp_ix);

end