| [SP_local_max2, indices2, SP_peaka, index_peaka, t_SP, gp ]=localpeaks(SP, Fs, num_pts_per_pk_intrl, min_peak, ignore_threshold, Align_peaks, make_plot, peak_alignment_tol)
|
function [SP_local_max2, indices2, SP_peaka, index_peaka, t_SP, gp ]=localpeaks(SP, Fs, num_pts_per_pk_intrl, min_peak, ignore_threshold, Align_peaks, make_plot, peak_alignment_tol)
% % localpeaks: Finds impulsive peaks and returns the amplitudes and indices
% %
% % Syntax:
% %
% % [SP_local_max2 indices2, SP_peaka, index_peaka]=localpeaks(SP, Fs,...
% % num_pts_per_pk_intrl, min_peak, ignore_threshold, Align_peaks, make_plot, ...
% % peak_alignment_tol);
% %
% % ********************************************************************
% %
% % Description
% %
% % This progam finds the impulsive peaks and determines whether the
% % impulsive peak criteria are satisfied.
% %
% % The program can handle large data arrays with multiple channels.
% %
% % ********************************************************************
% %
% % Input Variables
% %
% % SP=randn(4, 100000); % (Pa) is the Sound pressure data in Pa.
% % % default value is randn(4, 100000).
% %
% % Fs=100000; % (Hz) is the sampling rate.
% % % default value is 100000 Hz.
% %
% % num_pts_per_pk_intrl=100000;
% % % is the number of data points used to idntify a single
% % % impulsive peak. Usually 1 second to 1/2 a
% % % second of data points works well.
% % %
% % % default value is 1 second of data points or the
% % % length of the SP whichever is shorter.
% %
% % min_peak=100; % (dB) is the lowest amplitude for a peak to be considered
% % % impulsive.
% % % The default value is min_peak=100; %(dB).
% %
% % ignore_threshold=0; % 1 surpresses the threshold and other requirements
% % % for selecting peaks
% % % ignore_threshold=1; % generally will find more peaks.
% % %
% % % 0 is the default. Returns typical number of
% % % peaks.
% % %
% % % default is ignore_threshold=0;
% %
% % Align_peaks=1; % is the channel number for the impulsive peaks to be
% % % aligned to. default value is [];
% % % Align_peaks=1; peaks are aligned to channel 1.
% % % Align_peaks=[]; peaks are not aligned.
% % %
% % % If Align_peaks > num_channels then
% % % Align_peaks is set to [];
% % %
% % % If Align_peaks < num_channels then
% % % Align_peaks is set to [];
% % % The default value is [];
% %
% % make_plot=1; % 1 Plots the processing data of the peaks.
% % % Otherwise no plots are made.
% % %
% % % default is make_plot=0;
% %
% % peak_alignment_tol=0.25;
% % % This input sets the fraction of the datapoints
% % % of num_pts_per_pk_intrl that a peak can be offset.
% % % When aligning peaks, the channels are
% % % synchronized, but temorally the peaks may be
% % % misaligned by a few milliseconds. The
% % % alignment tolerance allows the program to find
% % % the maximum peak within the alignment
% % % tolerance.
% % %
% % % default is peak_alignment_tol=0.25;
% %
% %
% % % Further explanation of peak_alignment_tol
% %
% % % One-half of the datapoints of num_pts_per_pk_intrl are
% % % used for finding the index of the maximum sound
% % % pressure amplitude.
% %
% % % The impulsive peak is found by selecting the highest
% % % sound pressure amplitude starting at the center of
% % % the bin and working toward the maximum radius which
% % % is one-half of the data points for the impulse.
% %
% % ********************************************************************
% %
% % Output Variables
% %
% % SP_local_max2 is a cell array of peak pressures in Pa for each
% % impulsive peak and for each channel
% %
% % indices2 is the cell arrray of indices for each impulsive peak and for
% % each channel
% %
% % SP_peaka is the cell array of the peak amplitudes that are greater
% % than the threshold and before implementing the min_peak and other
% % selection criterion.
% %
% % index_peaka is the cell array of the indices corresponding to
% % SP_peaka.
% %
% % ********************************************************************
% %
%
% 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);
% % Pa sound pressure time record
%
% Fs=100000; % Hz sampling rate
%
% num_pts_per_pk_intrl=50000; % number of data points used to identify
% % a single impulsive peak usually 1 second to 1/2 a
% % second of data points works well
%
% min_peak=100; % (dB) is the peak value that is the minimum
% % for being considered an impulsive noise.
%
% ignore_threshold=0; % 1 surpresses the threshold and other requirements
% % for selecting peaks
% % ignore_threshold=1; % generally will find more peaks.
% %
% % 0 is the default. Returns typical number of
% % peaks.
% %
% % default is ignore_threshold=0;
%
% Align_peaks=1; % Specify the channel number to align peaks to
% % For the subordinate channels, the highest
% % peak within a radius of 0.25*num_pts_per_pk_intrl
% % data points of the impulsive peak of the
% % align_peaks channel is selected.
%
% make_plot=1; % 1 for make the plot
% % 0 for no plot
% % Default is to not make any plots.
%
% [SP_local_max2 indices2]=localpeaks(SP, Fs, num_pts_per_pk_intrl, min_peak, ignore_threshold, Align_peaks, make_plot);
%
% % ********************************************************************
% %
% %
% % Subprograms
% %
% %
% % List of Dependent Subprograms for
% % localpeaks
% %
% % FEX ID# is the File ID on the Matlab Central File Exchange
% %
% %
% % Program Name Author FEX ID#
% % 1) convert_double Edward L. Zechmann
% % 2) dB_to_Pa Edward L. Zechmann
% % 3) estimatenoise John D'Errico 16683
% % 4) findjobj Yair M. Altman 14317
% % 5) fix_YTick Edward L. Zechmann
% % 6) func_threshold Jing Tian 10462
% % 7) LMSloc Alexandros Leontitsis 801
% % 8) maximize Oliver Woodford 25471
% % 9) moving Aslak Grinsted 8251
% % 10) Pa_to_dB Edward L. Zechmann
% % 11) parseArgs Malcolm Wood 10670
% % 12) peak_index Edward L. Zechmann
% % 13) peak_ix Edward L. Zechmann
% % 14) peak_threshhold_function2 Edward L. Zechmann
% % 15) peakfinder Nate Yoder 25500
% % 16) psuedo_box Edward L. Zechmann
% % 17) resample_plot Edward L. Zechmann
% % 18) sub_mean Edward L. Zechmann
% % 19) subaxis Aslak Grinsted 3696
% % 20) wsmooth Damien Garcia NA
% %
% %
% %
% % ********************************************************************
% %
% % localpeaks is written by Edward L. Zechmann
% %
% % date 10 August 2007
% %
% % modified 18 December 2007 Improved signal processing methods
% % Improved the empty array handling
% % Added align peaks
% % Added analytic_impulse example
% %
% % modified 21 December 2007 Improved plotting the peak threshold
% % function
% % Improved the zero crossing and slope
% % calculation for removing peaks at teh end
% % of the file
% % Added more comments
% %
% % modified 27 December 2007 Changed the process of grouping peaks
% % from marchinig in time to selecting the
% % highest peak to the lowest peak
% % Fixed pts_per_bin to number of points
% % in 0.01 seconds
% %
% % modified 20 February 2008 Updated comments
% %
% %
% % modified 26 February 2008 Changed the peak_index for the align
% % peaks to use a radius of
% % 0.25*num_pts_per_pk_intrl data points.
% % Now the peak circles of the protected
% % signals are located at the peaks.
% %
% % modified 27 February 2008 Made the threshold limit for finding
% % a peak vary with time.
% %
% % modified 11 March 2008 The process of selecting detecting
% % impulsive peaks was greatly improved.
% %
% % modified 12 March 2008 Updated Comments
% %
% % modified 21 March 2008 Updated Comments
% %
% % modified 24 March 2008 Added radius option for
% % selecting the size of the range
% % of data points to find the peak
% % sound pressure amplitude.
% %
% % modified 22 July 2008 Modified Method of finding
% % peaks now uses percentiles.
% %
% % modified 19 August 2008 Added input variable percent1
% % to specify threshold for
% % identifying peaks
% %
% % modified 10 September 2008 Added optional method for finding
% % impulsive Exhaust Noise
% %
% % modified 11 November 2008 Fixed a bug with the exhaust noise
% % feature and with bin size.
% %
% % modified 6 October 2009 Updated comments
% %
% % modified 11 October 2010 Changed peak finder to peak_ix
% % Updated comments.
% %
% % modified 11 February 2011 Updated to new method of finding peaks.
% % Updated Comments.
% %
% %
% % ********************************************************************
% %
% % Please feel free to modify this code.
% %
% % See Also: Impulsive_Noise_Meter, snd_peak_metrics, plot_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);
% number of points in each peak interval
if nargin < 3 || isempty(num_pts_per_pk_intrl) || ~isnumeric(num_pts_per_pk_intrl)
n1=length(SP);
if Fs < n1
num_pts_per_pk_intrl=Fs;
else
num_pts_per_pk_intrl=n1;
end
end
% minimum peak preessure to be an impulsive noise dB
if nargin < 4 || isempty(min_peak) || ~isnumeric(min_peak)
min_peak=100;
end
% ignore_threshold is suspend absolute peak requirement
% a value of zero will not supress any requirements
if nargin < 5 || isempty(ignore_threshold) || ~isnumeric(ignore_threshold)
ignore_threshold=0;
end
% minimum peak preessure to be an impulsive noise
if nargin < 6 || isempty(Align_peaks) || ~isnumeric(Align_peaks)
Align_peaks=[];
end
if ~isempty(Align_peaks)
Align_peaks=round(Align_peaks);
if Align_peaks < 1
Align_peaks=[];
end
if Align_peaks > m1
Align_peaks=[];
end
end
% make a plot of the bin selection criterion
if nargin < 7 || isempty(make_plot) || ~isnumeric(make_plot)
make_plot=0;
end
% set the peak alignment tolerance
if nargin < 8 || isempty(peak_alignment_tol) || ~isnumeric(peak_alignment_tol)
peak_alignment_tol=0.25;
end
% Set no_smoothing to 1.
no_smoothing=1;
sh=0.0;
sv=0.0;
ml=0.14;
mr=0.1;
mt=0.08;
mb=0.12;
% 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 the number of pts in 0.01 seconds
pts_per_bin=floor(Fs/1000);
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
% Find the peaks using peak_ix
if isequal(ignore_threshold, 1)
[indices, gp]=peak_ix(SP, Fs, -inf, 0, 1);
else
[indices, gp]=peak_ix(SP, Fs, min_peak, 0, no_smoothing);
end
[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)]);
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*SP))));
xlim1=n1/Fs;
h2=[];
tot_num_samples=mm1*num_bins;
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=n1/num_bins/Fs*(0:(num_bins-1));
if flag_rs == 1
[t_SP_rs, SP_rs]=resample_plot(t_SP, SP(e1, :));
plot(t_SP_rs, SP_rs);
else
t_SP=1/Fs*(-1+(1:length(SP)));
plot(t_SP, SP(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;
% The peaks that are above the threshold are selected
plot(indices{e1, 1}, SP(e1,indices{e1, 1}), 'linestyle', 'none', 'marker', 'o', 'markeredgecolor', 'k' );
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, :)), floor(indices{e1,1}./pts_per_bin), pts_per_bin, radius*num_pts_per_pk_intrl);
SP_peaka{e1, 1}=SP_peak;
index_peaka{e1, 1}=index_peak;
end
min_peak_Pa=dB_to_Pa(min_peak);
% Enforce the min_peak requirement
% index_peaka now only includes peaks above the min_peak amplitude (Pa)
if ~isequal(ignore_threshold, 1)
for e1=1:m1;
index1=find(abs(SP_peaka{e1, 1}) > min_peak_Pa);
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 the 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
if ~isempty(index_maxa)
[sorted_mp_ix]=sort(index_maxa);
SP_local_max=SP(e2, sorted_mp_ix);
else
SP_local_max=[];
end
end
SP_local_max2{e2, 1}=SP_local_max;
indices2{e2, 1}=shiftdim(sorted_mp_ix);
end
|
|
