peak_threshhold_function.m

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

function gp=peak_threshhold_function(SP, Fs, pts_per_bin, exhaust_cycle)
% % peak_threshhold_function: Calculates a threshold for finding peaks
% % 
% % Syntax:
% % 
% % gp=peak_threshhold_function(SP, Fs, pts_per_bin);
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % Description
% %  
% % In complex noise, continuous noise and impulsive noise are combined.
% % To make it easier to identify impulsive peaks.  The time record is
% % processed in severl steps. 
% %
% % Step  Description
% % 
% % 1)    Remove the continuous noise.   
% % 2)    Subtract the Minimum Value from Each Channel.
% % 3)    Break up time record into bins and select the peak from each bin
% % 4)    Apply the gradient operator twice to sharpen the peaks
% % 
% % This process makes it easier to apply a threshold to identify peaks.  
% % The impulses are typically significantly greater than the 
% % remaining noise.  
% % 
% % There is difficulty in identifying impulses generated by 
% % white noise sources such as intermittent exhaust noise.  The first step
% % attempts to eliminate the any contrbution due to continuous noise.  
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % Input Variables
% %
% % SP is the array of sound pressure data in Pa.  
% %      should have size [num channels length_time_record]
% %      assumes the length of the tiem record is greater than the number
% %      of channels and tranposes SP if necessary
% %      default value is randn(4, 100000).
% % 
% % Fs sampling rate in Hz.  default value is 100000 Hz.  
% % 
% % pts_per_bin is the number of points in the sub bins for speedings up 
% %      the process of analyzing the data array.  The pts_per_bin is set 
% %      so that there are 100 bins peak peak interval. 
% %      The default value is pts_per_bin=floor(Fs/100);.
% % 
% % 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;
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % Output Variables
% % 
% % gp is a matrix of size [num_channels pts_per_bin].  gp is formed by 
% %      processing SP.  The purpose is to maximize the likelyhood of 
% %      detecting a peak while minimizing teh likelihood of falsely 
% %      detecting a peak.   
% %      gp makes identifying impulsive peaks much easier than 
% %      applying a threshold on SP.   
% %      The highest values of gp correspond to impulsive noise peaks.  
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% 
% Example
%
% Fs=100000; fc=1000; td=1; tau=0.01; delay=0.1; A1=5; A2=20;
% [SP, t]=analytic_impulse(Fs, fc, td, tau, delay, A1, A2);
% pts_per_bin=floor(Fs/100);
%
% gp=peak_threshhold_function(SP, Fs, pts_per_bin);
%
% 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % 
% % List of Dependent Subprograms for 
% % peak_threshhold_function
% % 
% % FEX ID# is the File ID on the Matlab Central File Exchange
% % 
% % Program Name   Author   FEX ID#
% % 1) convert_double		
% % 2) hilbert2		
% % 3) sub_mean	
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% % peak_threshhold_function is written by Edward L. Zechmann
% % 
% %     date     August     2008
% % 
% % modified  16 August     2008  Updated Comments
% % 
% % modified 10 September   2008  Added exhaust_cycle finding option
% %                               Updated Comemnts
% % 
% % modified 11 November   2008   Fixed a bug with the exhaust noise
% %                               feature and with bin size.  
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% % Please feel free to modify this code.
% % 
% % See Also: Impulsive_Noise_Meter, plot_peaks, local_peaks
% % 

if nargin < 1 || isempty(SP) || ~isnumeric(SP)
    SP=randn(4, 100000);
end

% make sure that input data is double precision.
[SP]=convert_double(SP);

% Make sure the data has the channels and data array in the
% correct order, transpose if necessary.
[m1 n1]=size(SP);

if m1 > n1
    SP=SP';
    [m1 n1]=size(SP);
end

if nargin < 2 || isempty(Fs) || ~isnumeric(Fs)
    Fs=100000;
end

[Fs]=convert_double(Fs);


if nargin < 3 || isempty(Fs) || ~isnumeric(Fs)
    pts_per_bin=floor(Fs/100);

    if pts_per_bin < 1
        pts_per_bin=1;
    end

    if pts_per_bin > n1
        pts_per_bin=n1;
    end
end

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


% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% STEP 1 Remove the Continuous Noise. 
% 
% SP2 is the instantaneous amplitude of the sound pressure
SP=SP-mean(SP, 2)*ones(1, n1);
if ~isequal(exhaust_cycle, 1)
    try
        SP2 = abs(hilbert(SP')');
    catch
        SP2 = abs(hilbert2(SP')');
    end
    
end

% SP2 is the time average instantaneous sound pressure amplitude
% with 2000 averages per second
if ~isequal(exhaust_cycle, 1)
    [buf, SP2]=sub_mean(SP2, Fs, 2000);
else
    SP2=abs(SP);
end


% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% STEP 2 Subtract the Minimum Value from Each Channel.
%
% Now the continuous noise has been greatly reduced.
% Subtract the minimum value from each channel. 
% So that the minimum value for each channel is now 0. 
% This will make the plots look better.
% Also it is better for locating peaks
SP2=SP2-min(SP2, [], 2)*ones(1, n1);

% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% STEP 3 Break up time record into bins and select the peak from each bin
% 
% Now the peaks are of primary interest.
% 
% Break up the processed time record into 100 bins per second
% and calculate some statistics for each bin.  
% The most important metric is the maximum peak in each bin.
n_bins=floor(n1/pts_per_bin);
SP_bins_peak=zeros(m1, n_bins);

for e1=1:n_bins;
    for e2=1:m1;
        
        SP_bins_peak( e2, e1)=max(max(abs(SP2( e2, (e1-1)*pts_per_bin+(1:pts_per_bin) ))));
        
    end
end


% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% STEP 4 Apply the gradient operator twice to sharpen the peaks
% 
% Now the peakedness of the signal must be maximized to expose all 
% underlying impulsive peaks.
% 
% The exponentially decaying impulses can have very long tails which make 
% it difficult to detect other impulses.  
% 
% To make it easier to identify impulsive peaks, further process the data 
% to narrow the exponentially decaying tails by taking the derivative with 
% respect to time.  The derivative is caculated as many times as possible 
% currently the gradient program only supports taking the derivative 
% twice.  
% 
% gp is the variable that is used to detect peaks
% gp is tha array of the numerical laplacian of the peaks and is 
% referred to as the "array of the peaks" in the rest of the comments.  
if isequal(exhaust_cycle, 1)
    gp=abs(SP_bins_peak);
else
    gp=abs(gradient(gradient(SP_bins_peak)));
end