| [metrics, metric_str, metric_units, AHAAH_str, round_kind, round_digits, dB_or_linear]=snd_peak_metrics(SP, Fs, indices2, num_pts_per_pk_intrl, N, min_f, max_f, ear_model, Polarized_or_ICP)
|
function [metrics, metric_str, metric_units, AHAAH_str, round_kind, round_digits, dB_or_linear]=snd_peak_metrics(SP, Fs, indices2, num_pts_per_pk_intrl, N, min_f, max_f, ear_model, Polarized_or_ICP)
% % snd_peak_metrics: Calculates impulsive sound metrics including Nth octave band Leqs and peak levels
% %
% % Syntax:
% %
% % [metrics]=snd_peak_metrics(SP, Fs, indices2, num_pts_per_pk_intrl, N, min_f, max_f, ear_model, Polarized_or_ICP);
% %
% % ********************************************************************
% %
% % Description
% %
% % This program calculates metrics for impulsive sound including
% % peak, peakA, peakC, spl_peak, spl_peakA, spl_peakC, LeqA, LeqA8, LeqC,
% % LeqC8, Leq, Leq8, num_pts_per_pk_intrl/Fs, [a, b_mil, c, d, e], kurt1,
% % Nth octave band peak levels, and time average sound pressure levels.
% %
% % The input and output variables are described below
% %
% % ********************************************************************
% %
% % Input Variables
% %
% % SP=randn(10, 50000); % (Pa) is the time record of the sound pressure
% % % default value is rand(1, 50000);.
% %
% % Fs=50000; % sampling rate in Hz.
% % % default value is 100000 Hz.
% %
% % indices2=ones(10,1); % cell array of indices of the peaks. Use
% % % local_peaks.m to find the peaks.
% % %
% % % default value is
% % % indices2=cell(m1, 1);
% % % for e1=1:m1;
% % % indices2{e1, 1}=1;
% % % end
% %
% % num_pts_per_pk_intrl=100000;
% % % Is the number of data points used to
% % % characterize each 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.
% % % default is all of the datapoints.
% %
% % N=3; % number of bands per octave
% % % Program Now supports Nth octave band filters
% % %
% % % N=3; is one-third octave
% % % N=1; is full octave
% % % N=12; is one-twelveth octave
% % % default is 3;
% %
% % min_f=200; % (Hz) this is the center frequency of the
% % % lowest frequency to analyze with the
% % % third octave peaks.
% % % default is 200;
% %
% % max_f=20000; % (Hz) this is the center frequency of the
% % % highest frequency to analyze with the
% % % third octave peaks.
% % % default is 20000;
% %
% % ear_model=ones(size(SP, 1)); % Not Supported
% % % AHAAH Model is not Currently
% % % Available.
% % %
% % % ear-model is a vector of integers.
% % % this vector specifies which transfer
% % % function the ahaah model will compute
% % % each of the transfer functions starts
% % % with part of the ear and finishes at
% % % the hair cells.
% % %
% % % 1 outer ear to hair cell
% % % 2 middle ear to hair cell
% % % 3 inner ear to hair cell
% %
% % Polarized_or_ICP
% %
% %
% % ********************************************************************
% %
% % Output Variables
% %
% % metrics includes all of the sound pressure metrics calculated by this
% % program the metrics are described below
% %
% % There are 20 metrics + the third octave band metrics + Ahaaha
% % Warned and Unwarned will be added soon.
% %
% % 'Peak Index', (index)
% % 'Peak Time', (s)
% %
% % The Peak Amplitude Pressure is recorded in Pa with the metric
% % names
% % 'Peak Pres Linear', (Pa)
% % 'Peak Pres A-Weight' (Pa)
% % 'Peak Pres C-Weight' (Pa)
% %
% % The Peak Level Pressure is recorded in Pa with the metric names
% % 'Peak Level Linear' (dB)
% % 'Peak Level A-Weight' (dB)
% % 'Peak Level C-Weight' (dB)
% %
% % The Time Avereaged Pressure Level is recorded in dB with the
% % metric names
% % 'LeqA' (dB)
% % 'LeqA8' (dB)
% % 'LeqC' (dB)
% % 'LeqC8' (dB)
% % 'Leq' (dB)
% % 'Leq8' (dB)
% %
% % The total time span of the data points used to characterize the
% % impulse is recorded in (s) with the metric name
% % 'Time Span' (s)
% %
% % rt is the rise time in (s)
% %
% % The Impulsive time durations are recorded in (s) with the
% % metric names
% % 'A-Duration' (s)
% % 'B-MIL-STD-1474D-Duration' (s)
% % 'C-Duration' (s)
% % 'D-Duration' (s)
% % 'E-Duration' (s)
% %
% % 'Kurtosis' (unitless)
% %
% % Nth octave peaks (dB) from min_f to max_f (Hz)
% %
% % Nth octave levels (dB) from min_f to max_f (Hz)
% %
% %
% % AHAAH Model is not currently Available.
% % Ahahha model selected model type for each channel units are in
% % Auditory Hazard Units (AHU)
% %
% % AHAAH_str is a cell array of strings of size [num_channels, 1].
% % Each cell describes which the contents of
% %
% %
% % round_kind % Array of values one element for the rta array
% % % and one element for each varargin array
% % % (see example)
% % % 1 round to specified number of significant
% % % digits
% % %
% % % 0 round to specified digits place
% %
% % round_digits % Array of values one element for the rta array
% % % and one element for each varargin array
% % % (see example)% Type of rounding depends on round_kind
% % %
% % % if round_kind==1 number of significant digits
% % % if round_kind==0 specified digits place
% %
% % dB_or_linear % 1 if in logarithmic units 0 if in linear units
% %
% % ********************************************************************
% %
% % Ahaah Note:
% %
% % The Ahaah model is a lumped parameter model of how the human hearing
% % mechanism. It works from the outer ear to the inner ear hair cell
% % stimulations.
% %
% % Dick Price and Joel Calb co-invented the model in the 1980's.
% % Their research was funded by the US Army.
% %
% % ********************************************************************
%
%
% Example='1';
%
% m1=10; n1=50000;
% SP=randn(m1, n1);
% Fs=50000;
% indices2=cell(m1, 1);
% for e1=1:m1;
% indices2{e1, 1}=1;
% end
% num_pts_per_pk_intrl=n1;
% N=3;
% min_f=200;
% max_f=20000;
% ear_model=[1,2,3,1,2,3,1,2,3,1];
%
% [metrics]=snd_peak_metrics(SP, Fs, indices2, num_pts_per_pk_intrl, N, ...
% min_f, max_f, ear_model);
%
% 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, 10*A2);
% [bb ix]=max(abs(SP1));
% num_pts_per_pk_intrl=100000; N=3; min_f=200; max_f=20000; ear_model=1;
%
% Polarized_or_ICP=0; % assume an ICP microphone
%
% [metrics, metric_str, metric_units, AHAAH_str]=snd_peak_metrics(SP1, Fs_SP, {ix}, num_pts_per_pk_intrl, N, min_f, max_f, ear_model, Polarized_or_ICP);
%
%
% % ********************************************************************
% %
% %
% % References:
% %
% %
% % ANSI S1.4-1983 American National Standard
% % Specificaitons for Sound Level Meters
% %
% % ANSI S1.6-1984 Preferred Frequencies, Frequency Levels, and Band
% % Numbers for Acoustical Measurements
% %
% % ANSI S1.43-1997 American National Standard
% % Specificaitons for Integrating-Averaging Sound Level
% % Meters
% %
% % ANSI S1.11-1986 American National Standard
% % Specification for Octave-band and Fractional-Octave-
% % band Analog and Digital Filters
% %
% % ANSI S3.44-1986 American National Standard
% % Determination of Occupational Noise Exposure and
% % Estimation of Noise Induced Hearing Impairment
% %
% %
% % Mil Standard 1474D, DEPARTMENT OF DEFENSE
% % DESIGN CRITERIA STANDARDm, NOISE LIMITS, 12 February 1997
% % www.silencertests.com/docs/mil-std-1474d.pdf
% %
% % Guido F. Smoorenburg, "Damage Risk Criteria for Impulsive
% % Noise," New Perspectives on Noise Induced Hearing Loss,
% % Raven Press, New York, pages(471-490) 1982
% %
% %
% % ********************************************************************
% %
% %
% % Subprograms
% %
% % This program requires the Matlab Signal Processing Toolbox
% % Some subprograms are based on the Octave Toolbox by Christophe Couvreur
% % Matlab Central File Exchange ID 69
% %
% %
% %
% %
% % List of Dependent Subprograms for
% % snd_peak_metrics
% %
% % FEX ID# is the File ID on the Matlab Central File Exchange
% %
% %
% % Program Name Author FEX ID#
% % 1) A_duration Edward L. Zechmann
% % 2) abcd_durations Edward L. Zechmann
% % 3) ACdsgn Edward L. Zechmann
% % 4) ACweight_time_filter Edward L. Zechmann
% % 5) analytic_impulse Edward L. Zechmann
% % 6) B_Duration_second_flucts Edward L. Zechmann
% % 7) B_mil_1474D_duration Edward L. Zechmann
% % 8) bessel_antialias Edward L. Zechmann
% % 9) bessel_digital Edward L. Zechmann
% % 10) bessel_down_sample Edward L. Zechmann
% % 11) C_duration Edward L. Zechmann
% % 12) convert_double Edward L. Zechmann
% % 13) D_duration Edward L. Zechmann
% % 14) E_duration Edward L. Zechmann
% % 15) filter_settling_data3 Edward L. Zechmann
% % 16) find_previous_crossing Edward L. Zechmann
% % 17) findextrema Schuberth Schuberth 3586
% % 18) geomean2 Edward L. Zechmann
% % 19) geospace Edward L. Zechmann
% % 20) get_p_q2 Edward L. Zechmann
% % 21) kurtosis2 William Murphy
% % 22) Leq_all_calc Edward L. Zechmann
% % 23) LMS_trim Edward L. Zechmann
% % 24) LMSloc Alexandros Leontitsis 801
% % 25) LMTSregor Edward L. Zechmann
% % 26) m_round Edward L. Zechmann
% % 27) match_height_and_slopes2 Edward L. Zechmann
% % 28) moving Aslak Grinsted 8251
% % 29) nth_freq_band Edward L. Zechmann
% % 30) Nth_oct_time_filter2 Edward L. Zechmann
% % 31) Nth_octdsgn Edward L. Zechmann
% % 32) pow10_round Edward L. Zechmann
% % 33) rand_int Edward L. Zechmann
% % 34) remove_filter_settling_data Edward L. Zechmann
% % 35) resample_interp3 Edward L. Zechmann
% % 36) rise_time Edward L. Zechmann
% % 37) rms_val Edward L. Zechmann
% % 38) sd_round Edward L. Zechmann
% % 39) sub_mean Edward L. Zechmann
% % 40) sub_mean2 Edward L. Zechmann
% % 41) threshold_bin_peaks Edward L. Zechmann
% %
% %
% %
% % ********************************************************************
% %
% %
% % Program was written by Edward L. Zechmann
% %
% % date 16 October 2007
% %
% % modified 13 January 2008 Updated comments
% % changed from filter.m to filtfilt.m
% % Added impulsive noise method
% %
% % modified 14 January 2008 Reduced number of points allowed for
% % the robust line fit
% % Removed the time truncation
% % updated comments
% %
% % modified 13 August 2008 Added third octave peak levels
% % Added third octave time average levels
% % Updated comments
% %
% % modified 18 August 2008 Added Ahaah model developed by US Army
% % Updated comments
% %
% % modified 10 September 2008 Updated comments
% %
% % modified 10 December 2008 Upgraded the Nth octave band
% % filtering programs and the A and C-
% % weighting filter programs to include
% % filter settling and resampling.
% % Previously only suppoerted
% % 3rd octave band filters.
% %
% % modified 11 December 2008 Upgraded the A and C-
% % weighting filter programs,
% % to include iterative filtering.
% % The filters are now very stable.
% %
% % Removed filter coefficients from input
% % and output;
% % Peaks pressures and Levels are output.
% %
% % modified 16 December 2008 Use convolution to make filter
% % coefficients (b and a) into
% % arrays from cell arrays.
% %
% % modified 11 January 2009 Added rounding.
% %
% % modified 21 March 2009 Updated comments
% %
% % modified 6 October 2009 Updated comments
% %
% % modified 5 August 2010 Added resampling using Bessel
% % antialiasing filter
% % Updated Comments
% %
% % modified 10 August 2010 Added rise time
% % Updated Comments
% %
% % modified 11 February 2011 Fixed bug in outputting the E-Duration
% % Updated Comments
% %
% %
% %
% % ********************************************************************
% %
% %
% % Please feel free to modify this code.
% %
% % See Also: Impulsive_Noise_Meter, plot_peaks, local_peaks, octave, oct3dsgn
% %
% Check for erroneous input and use reasonable values
if (nargin < 1 || isempty(SP)) || ~isnumeric(SP)
SP=randn(1, 50000);
end
[SP]=convert_double(SP);
[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(indices2))
indices2=cell(m1, 1);
for e1=1:m1;
indices2{e1, 1}=1;
end
end
if (nargin < 4 || isempty(num_pts_per_pk_intrl)) || ~isnumeric(num_pts_per_pk_intrl)
num_pts_per_pk_intrl=n1;
end
% Program Now Supports Nth Octave Band Filtering
if (nargin < 5 || isempty(N)) || ~isnumeric(N)
N=3;
end
if (nargin < 6 || isempty(min_f)) || ~isnumeric(min_f)
min_f=200;
end
if (nargin < 7 || isempty(max_f)) || ~isnumeric(max_f)
max_f=20000;
end
if (nargin < 8 || isempty(max_f)) || ~isnumeric(max_f)
ear_model=ones(m1, 1);
end
low_bin=floor(0.25*num_pts_per_pk_intrl);
high_bin=num_pts_per_pk_intrl-low_bin;
metrics=cell(m1, 1);
SI_prefixes=0;
[fc, fc_l, fc_u, fc_str] = nth_freq_band(N, min_f, max_f, SI_prefixes);
array_names= {'Peak Index', 'Peak Time', 'Peak Pres. Linear', 'Peak Pres. A-Weight', 'Peak Pres. C-Weight', 'Peak Level Linear ', 'Peak Level A-Weight ', 'Peak Level C-Weight ', 'LeqA', 'LeqA8', 'LeqC', 'LeqC8', 'Leq', 'Leq8', 'Time Span', 'Rise Time', 'A-Duration', 'B-MIL-STD-1474D-Duration', 'C-Duration', 'D-Duration', 'E-Duration', 'Kurtosis'};
array_units= {'(Indices)', '(s)', '(Pa)', '(Pa)', '(Pa)', '(dB)', '(dBA)', '(dBC)', '(dBA)', '(dBA)', '(dBC)', '(dBC)', '(dB)', '(dB)', '(s)', '(s)', '(s)', '(s)', '(s)', '(s)', '(s)', '(No units)'};
num_snd_simple_metrics=length(array_names);
num_metrics=num_snd_simple_metrics+2*length(fc);
metric_str=cell(num_metrics, 1);
metric_units=cell(num_metrics, 1);
dB_or_linear=zeros(1, num_metrics);
dB_or_linear(1, 1:num_snd_simple_metrics)=[0 0 0 0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 ];
% append the strings for the simple sound metrics
for e1=1:num_snd_simple_metrics;
metric_str{e1}=array_names{e1};
metric_units{e1}=array_units{e1};
end
% append the strings for the third octave band sound metrics
for e1=1:length(fc);
metric_str{num_snd_simple_metrics+e1}=['Peak ', fc_str{e1}, ' Hz'];
metric_str{num_snd_simple_metrics+length(fc)+e1}=['Level ', fc_str{e1}, ' Hz'];
metric_units{num_snd_simple_metrics+e1}='(dB)';
metric_units{num_snd_simple_metrics+length(fc)+e1}='(dB)';
dB_or_linear(num_snd_simple_metrics+e1)=1;
dB_or_linear(num_snd_simple_metrics+length(fc)+e1)=1;
end
% append the strings for the third octave band sound metrics
% metric_str{num_snd_simple_metrics+2*length(fc)+1}='AHAAH Warned';
% metric_str{num_snd_simple_metrics+2*length(fc)+2}='AHAAH UnWarned';
% metric_units{num_snd_simple_metrics+2*length(fc)+1}='(AHU)';
% metric_units{num_snd_simple_metrics+2*length(fc)+2}='(AHU)';
ahaah_strs={'Outer Ear', 'Middle Ear', 'Inner Ear'};
AHAAH_str=cell(m1, 1);
for e2=1:m1; % number of microphone channels
max_size=length(indices2{e2}); % number of impulsive peaks found
buffer4=zeros(max_size, num_metrics);
if max_size >= 1
for e1=1:max_size;
index1=indices2{e2}(e1)-low_bin+1;
if index1 < 0
index1=1;
end
index2=index1+num_pts_per_pk_intrl-1;
if (index2 > n1)
index2 = n1;
end
if index2 <= n1
p=SP(e2, index1:index2);
% **********************************************************
%
% Calculate the A-weighted, C-weighted, and Linear-Weighted
% sound pressure levels (dB) and peak amplitudes (Pa)
%
% Set the settling time to 0.1 seconds which is a typical
% value.
settling_time=0.1;
% Set the calibration factor to 1.
cf=1;
resample_filter=2;
[LeqA, LeqA8, LeqC, LeqC8, Leq, Leq8, spl_peak, spl_peakA, spl_peakC, peak, peakA, peakC]=Leq_all_calc(p, Fs, cf, settling_time, resample_filter);
% **********************************************************
%
%
% rise time
percent_low=10; % percentage of the range from 0 to the peak level
% 10 percent is the typical value from the Federal
% Telecommunications Glossary on rise time.
percent_high=90; % percentage of the range from 0 to the peak level
% 90 percent is the typical value from the Federal
% Telecommunications Glossary on rise time.
peak_index=low_bin;
make_plot=0; % 1 Plots the processing data of the peaks.
[rt]=rise_time(p, Fs, percent_low, percent_high, peak_index, make_plot);
% **********************************************************
%
%
% Calculate the durations for each peak
[a, b_mil, c, d, e]=abcd_durations(p, Fs, 0, Polarized_or_ICP);
kurt1=kurtosis2(p, 2);
if isempty(kurt1)
kurt1=-1;
end
% **********************************************************
%
%
% Calculate the one-third octave linear time average levels and peak levels.
% Run the third octave band fitler twice.
num_x_filter=1;
sensor=1;
settling_time=0.1;
filter_program=2;
resample_filter=2;
[fc_out, SP_levels, SP_peak_levels]=Nth_oct_time_filter2(p, Fs, num_x_filter, N, fc, sensor, settling_time, filter_program, resample_filter);
num_fc=length(fc_out);
% **********************************************************
%
%
% AHAAH model metrics are not currently available.
%[err, WAHU, UwAHU] = ahaah(p, int32(fs), int32(ear_model(e2)));
AHAAH_str{e2}=ahaah_strs(ear_model(e2));
% Specify the kind of rounding.
round_kind= [0 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 zeros(1, num_fc) zeros(1, num_fc)];
% Specify the number of significant digits or the digits place.
round_digits=[0 3 3 3 3 0 0 0 -1 -1 -1 -1 -1 -1 4 3 3 3 3 3 3 3 zeros(1, num_fc) -ones(1, num_fc)];
% Concatenate all the metrics into one variable.
all_metrics=[indices2{e2}(e1), 1/Fs*indices2{e2}(e1), peak, peakA, peakC, spl_peak, spl_peakA, spl_peakC, LeqA, LeqA8, LeqC, LeqC8, Leq, Leq8, num_pts_per_pk_intrl/Fs, rt, a, b_mil, c, d, e, kurt1, SP_peak_levels, SP_levels];
% Apply te especified rounding to the metrics.
buffer4(e1, :)=m_round(all_metrics, round_kind, round_digits);
end
end
else
num_fc=length(fc);
% Specify the kind of rounding.
round_kind= [0 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 zeros(1, num_fc) zeros(1, num_fc)];
% Specify the number of significant digits or the digits place.
round_digits=[0 3 3 3 3 0 0 0 -1 -1 -1 -1 -1 -1 4 3 3 3 3 3 3 3 zeros(1, num_fc) -ones(1, num_fc)];
% Concatenate all the metrics into one variable.
all_metrics=[-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -ones(1, num_fc), -ones(1, num_fc)];
% Apply te especified rounding to the metrics.
buffer4(e1, :)=m_round(all_metrics, round_kind, round_digits);
end
metrics{e2, 1}=buffer4;
end
|
|
