Code covered by the BSD License  

Highlights from
Nth_Oct_Hand_Arm_&_AC_Filter_Tool_Box

Nth_Oct_Hand_Arm_&_AC_Filter_Tool_Box

by

 

22 Dec 2008 (Updated )

Features Nth octave band, Hand Arm, and A and C weighting filters

[f, A_atten, A_atten2, C_atten, C_atten2]=Test_ACweight(Fs, N, flag2, resample_filter)
function [f, A_atten, A_atten2, C_atten, C_atten2]=Test_ACweight(Fs, N, flag2, resample_filter)
% % Test_ACweight:  Tests the A and C weighting filters using pure tones and tone bursts
% %
% % Syntax;
% %
% % [f, A_atten, A_atten2, C_atten, C_atten2]=Test_ACweight(Fs, N, flag2, resample_filter);
% %
% % **********************************************************************
% %
% % Description
% %
% % Test_ACweight(Fs, N);
% % Returns four figures quantifying the performance of the A and C weight
% % filters and the error of the A and C-weighting filters respectively
% % compared to the Type0, 1, and 2 filter tolerance specifications of
% % ANSI S1.4.  The performance of the filters is calculated at the
% % sampling rate Fs (Hz) for a number of bands, N per octave.
% %
% % [f, A_atten, A_atten2, C_atten, C_atten2]=Test_Aweight2(Fs, N, flag2);
% % Returns the center frequency of each band and the theoretical and
% % actual attenuation for the A and C-weighting filters respectively.
% %
% % There are two options for the downsampling filters to optimize
% % performance for continuous signals or for impulsive signals. 
% % For continuous noise the time domain does not have significant 
% % impulses; however, for impulsive time records there are often very
% % large impulses with distinctive peaks.  
% % 
% % There are two antialiasing filters and interpolation schemes available.
% % The first program is the built-in Matlab "resample" progam which
% % uses a Kaiser window fir filter for antialising and uses an unknown 
% % interpolation method.  The second program available for downsampling 
% % is bessel_down_sample which uses a Bessel filter for antialiasing 
% % and uses interp with the cubic spline option for interpolation.  
% % 
% % The resample function has good antialising up to the Nyquist frequency;
% % however, it has significant ringing effect when there are impulses.  
% % The bessel_down_sample function has good antialising; however, there is
% % excessive attenuation near the Nyquist frequency.  
% % The bessel_down_sample function experiences no ringing due to impulses
% % so it is very useful for peak estimation.  
% %
% % The input and output variables are described in more detail in the
% % respective sections below.
% %
% % ***********************************************************
% %
% % Input Variables
% %
% % Fs=50000;           % (Hz) sampling rate in Hz.  
% %                     % default is 50000 Hz.
% % 
% % N=3;                % Number of frequency bands per octave.
% %                     % Can be any number > 0.  
% %                     % Default is N=3; (third octave bands)
% %
% % flag2=1;            % Boolean which determines whether to test rms 
% %                     % levels with sinusoids or peak levels with 
% %                     % ringing impulses from analytic_impulse.  
% %                     % flag2=1 tests the pure tones.  Otherwise the
% %                     % peak levels aerea tested.
% %                     % default is flag2=1; Testing pure tones.
% %
% % resample_filter=1;  % type of filter to use when resamling
% %                     % 1 resample function Kaiser window fir filter
% %                     % 2 Bessel filter 
% %                     % otherwise resample function Kaiser window fir
% %                     % filter
% %                     % default is resample_filter=1; (Kaiser window)
% %
% %
% %
% % ***********************************************************
% %
% % Output Variables
% %
% % f (Hz) is the center frequency bands of the filters.
% %
% % A_atten (dB) is the theoretical attenuation of the A-weighting filters.
% %
% % A_atten2 (dB) is the actual attenuation of the A-weighting filters.
% %
% % C_atten (dB) is the theoretical attenuation of the C-weighting filters.
% %
% % C_atten2 (dB) is the actual attenuation of the C-weighting filters.
% %
% %
% % **********************************************************************
%
% Example='1';
% % Test the filters for six sampling rates.
% 
% cd..
% mkdir('Resample');
% cd('Resample');
% Fsa=1000*[20 50 100 200 500 1000];
% N=3;
% flag2=1;
% resample_filter=1;
% for e1=1:length(Fsa);
%   close all;
%   Fs=Fsa(e1);
%   Test_ACweight(Fs, N, flag2, resample_filter);
% end
% 
% cd..
% mkdir('Bessel');
% cd('Bessel');
% resample_filter=2;
% for e1=1:length(Fsa);
%   close all;
%   Fs=Fsa(e1);
%   Test_ACweight(Fs, N, flag2, resample_filter);
% end
% 
% 
% 
% Example='2';
% % Test the filters for six sampling rates.
% Fsa=1000*[20 50 100 200 500 1000];
% N=12;
% flag2=1;
% for e1=1:length(Fsa);
%   close all;
%   Fs=Fsa(e1);
%   Test_ACweight(Fs, N, flag2);
% end
%
%
%
% Example='3';
% % Test the filters for six sampling rates.
% Fsa=1000*[20 50 100 200 500 1000];
% N=3;
% flag2=0;
% for e1=1:length(Fsa);
%   close all;
%   Fs=Fsa(e1);
%   Test_ACweight(Fs, N, flag2);
% end
%
%
%
% Example='4';
% % Test the filters for six sampling rates.
% Fsa=1000*[20 50 100 200 500 1000];
% N=12;
% flag2=0;
% resample_filter=1;
% for e1=1:length(Fsa);
%   close all;
%   Fs=Fsa(e1);
%   Test_ACweight(Fs, N, flag2, resample_filter);
% end
%
%
% % **********************************************************************
% % 
% % References
% % 
% % IEC/CD 1672: Electroacoustics-Sound Level Meters, Nov. 1996. 
% % 
% % ANSI S1.4: Specifications for Sound Level Meters, 1983. 
% % 
% % **********************************************************************
% %
% % 
% % 
% % Subprograms
% %
% % This program requires the Matlab Signal Processing Toolbox
% %
% % 
% % 
% % List of Dependent Subprograms for 
% % Test_ACweight
% % 
% % FEX ID# is the File ID on the Matlab Central File Exchange
% % 
% % 
% % Program Name   Author   FEX ID#
% %  1) AC_weight_NB		William Murphy			
% %  2) ACdsgn		Edward L. Zechmann			
% %  3) ACweight_time_filter		Edward L. Zechmann			
% %  4) bessel_antialias		Edward L. Zechmann			
% %  5) bessel_digital		Edward L. Zechmann			
% %  6) bessel_down_sample		Edward L. Zechmann			
% %  7) convert_double		Edward L. Zechmann			
% %  8) file_extension		Edward L. Zechmann			
% %  9) filter_settling_data3		Edward L. Zechmann			
% % 10) flat_top		Edward L. Zechmann			
% % 11) geospace		Edward L. Zechmann			
% % 12) get_p_q2		Edward L. Zechmann			
% % 13) Leq_all_calc		Edward L. Zechmann			
% % 14) LMSloc		Alexandros Leontitsis		801	
% % 15) match_height_and_slopes2		Edward L. Zechmann			
% % 16) moving		Aslak Grinsted		8251	
% % 17) nth_freq_band		Edward L. Zechmann			
% % 18) Nth_oct_time_filter2		Edward L. Zechmann			
% % 19) Nth_octdsgn		Edward L. Zechmann			
% % 20) number_of_averages		Edward L. Zechmann			
% % 21) pressure_spectra		Edward L. Zechmann			
% % 22) remove_filter_settling_data		Edward L. Zechmann			
% % 23) resample_interp3		Edward L. Zechmann			
% % 24) rms_val		Edward L. Zechmann			
% % 25) save_a_plot_reverb_time		Edward L. Zechmann			
% % 26) sd_round		Edward L. Zechmann			
% % 27) spatialPattern		Jon Yearsley		5091	
% % 28) spectra_estimate		Edward L. Zechmann			
% % 29) sub_mean		Edward L. Zechmann			
% % 30) tone_burst		Edward L. Zechmann			
% % 31) window_correction_factor		Edward L. Zechmann										
% %
% %
% % **********************************************************************
% %
% % Program was written by Edward L. Zechmann
% %
% % created                 2007
% %
% % modified 18 November    2008    Updated Comments
% %
% % modified 11 December    2008    Began modifying A and C weighting
% %                                 filters with resampling and filter
% %                                 settling.
% %
% % modified 17 December    2008    Added Type0, 1, and 2 tolerance curves
% %                                 Added ability to test peak levels.
% %
% % modified  2 August      2010    Modified filter settling data program
% %                                 Added resampling using a Bessel Filter.
% %                                 Updated Comments.
% %
% %
% % **********************************************************************
% %
% % Please feel free to modify this code.
% %
% % See Also: filter, filtfilt, resample, ACweight_time_filter,
% %           hand_arm_time_fil, whole_body_time_filter
% %


if (nargin < 1 || isempty(Fs)) || ~isnumeric(Fs)
    Fs=50000;
end

if (nargin < 2 || isempty(N)) || (logical(N < 0) || ~isnumeric(N))
    N=3;
end

if (nargin < 3 || isempty(N)) || (logical(N < 0) || ~isnumeric(N))
    flag2=1;
end

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


close('all');

% List the frequencies and the response tolerances for the Type0, Type 1,
% and Type 2 filters
ftype=[10 12.5 16 20 25 31.5 40 50 63 80 100 125 160 200 250 315 400 500 630 800 1000 1250 1600 2000 2500 3150 4000 5000 6300 8000 10000 12500 16000 20000];
a=[-70.4  -63.4 -56.7 -50.5 -44.7 -39.4 -34.6 -30.2 -26.2 -22.5 -19.1 -16.1 -13.4 -10.9 -8.6 -6.6 -4.8 -3.2 -1.9 -0.8 0 0.6 1.0 1.2 1.3 1.2 1 0.5 -0.1 -1.1 -2.5 -4.3 -6.6 -9.3];
%b=[-38.2 -33.2 -28.5 -24.2 -20.4 -17.1 -14.2 -11.6 -9.3 -7.4 -5.6 -4.2 -3.0 -2.0 -1.3 -0.8 -0.5 -0.3 -0.1 0 0 0 0 -0.1 -0.2 -0.4 -0.7 -1.2 -1.9 -2.9 -4.3 -6.1 -8.4 -11.1];
c=[-14.3 -11.2 -8.5 -6.2 -4.4 -3.0 -2.0 -1.3 -0.8 -0.5 -0.3 -0.2 -0.1 0 0 0 0 0 0 0 0 0 -0.1 -0.2 -0.3 -0.5 -0.8 -1.3 -2.0 -3.0 -4.4 -6.2 -8.5 -11.2];

% Response limits of the Type0 filters
type0_h= [2 2 2 2 1.5 1 1 1 1 1 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 1 1   1 2 2 2 2 ];
type0_l=-[5 4 3 2 1.5 1 1 1 1 1 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 1 1.5 2 3 3 3 3 ];

% Response limits of the Type1 filters
type1_h= [4 3.5 3 2.5 2 1.5 1.5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1.5 1.5 1.5 2 3 3   3];
type1_l=-[4 3.5 3 2.5 2 1.5 1.5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1.5 2   3   4 6 inf inf];

% Response limits of the Type2 filters
type2_h= [5   5   5   3 3 3 2 2 2 2 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5  2 2 2.5 2.5 3 3.5 4.5 5 5   5   5   5];
type2_l=-[inf inf inf 3 3 3 2 2 2 2 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5  2 2 2.5 2.5 3 3.5 4.5 5 inf inf inf inf];

% Calculate the center band frequencies of the signals
min_f=10;
max_f=min([Fs/3, 50000]);
[f] = nth_freq_band(N, min_f, max_f);

% Calculate the powers of 10 for plotting the frequency in the xtick
[f_plot, buf2, buf3, f_str] = nth_freq_band(3/10, min_f, max_f);


cf=1;
settling_time=0.1;

num_bands=length(f);

A_atten2=zeros(num_bands, 1);
C_atten2=zeros(num_bands, 1);


for e1=1:length(f);

    td=max([10/Fs, 10/f(e1), 0.2]);

    % Create the pure tone or a tone burst time record
    if isequal(flag2, 1)
        % produce a sinusoid with 10 full periods.
        t_SP=(0:(1/Fs):td);
        y=sin(2*pi*t_SP*f(e1));
    else
        % produce a pure tone.
        delay=0.1;
        A1=0;
        A2=1;
        num_waves=10;
        [y]=tone_burst(Fs, f(e1), td, delay, num_waves, A1, A2);
    end
    
    [LeqA, LeqA8, LeqC, LeqC8, Leq, Leq8, peak_dB, peak_dBA, peak_dBC, peak_Pa, peak_PaA, peak_PaC]=Leq_all_calc(y, Fs, cf, settling_time, resample_filter);

    if isequal(flag2, 1)
        A_atten2(e1)=LeqA-Leq;
        C_atten2(e1)=LeqC-Leq;
    else
        A_atten2(e1)=peak_dBA-peak_dB;
        C_atten2(e1)=peak_dBC-peak_dB;
    end

end

[Wa, Wc]=AC_weight_NB(f, 0);
A_atten=Wa;
C_atten=Wc;

LW=2;
if isequal(flag2, 1)
    title_blurb='RMS';
else
    title_blurb='Peak';
end

% Plot the A-weight Filter performance curve
figure(1);
semilogx(f, A_atten2, 'r', 'LineWidth', LW );
hold on;
semilogx(f, A_atten, 'b', 'LineWidth', LW );
semilogx(ftype, a+type0_h, '--g', 'LineWidth', LW);
semilogx(ftype, a+type1_h, '-.m', 'LineWidth', LW);
semilogx(ftype, a+type2_h, ':r', 'LineWidth', LW);
semilogx(ftype, a+type0_l, '--g', 'LineWidth', LW);
semilogx(ftype, a+type1_l, '-.m', 'LineWidth', LW);
semilogx(ftype, a+type2_l, ':r', 'LineWidth', LW);

legend({'Actual', 'Theoretical', 'Type-0', 'Type-1', 'Type-2'}, 'Location', 'South');
xlabel('Frequency Hz', 'Fontsize', 28);
ylabel('Attnuation dB', 'Fontsize', 28);
title( {[title_blurb, ' A-weight Filter'], [num2str(Fs/1000) , ' KHz Sampling Rate']},  'Fontsize', 40 );
set(gca, 'Fontsize', 28);
set( gca, 'XTickMode', 'manual', 'XTick', f_plot, 'XTickLabel', f_str);
xlim([min(f), 50000]);
ylim([-80 20]);


% Plot the C-weight Filter performance curve
figure(2);
semilogx(f, C_atten2, 'color', 'b', 'LineWidth', LW);
hold on;
semilogx(f, C_atten, 'k', 'LineWidth', LW);
semilogx(ftype, c+type0_h, '--g', 'LineWidth', LW);
semilogx(ftype, c+type1_h, '-.m', 'LineWidth', LW);
semilogx(ftype, c+type2_h, ':r', 'LineWidth', LW);
semilogx(ftype, c+type0_l, '--g', 'LineWidth', LW);
semilogx(ftype, c+type1_l, '-.m', 'LineWidth', LW);
semilogx(ftype, c+type2_l, ':r', 'LineWidth', LW);

legend({'Actual', 'Theoretical', 'Type-0', 'Type-1', 'Type-2'}, 'Location', 'South');
xlabel('Frequency Hz', 'Fontsize', 28);
ylabel('Attnuation dB', 'Fontsize', 28);
title( {[title_blurb, ' C-weight Filter'], [num2str(Fs/1000) , ' KHz Sampling Rate']},  'Fontsize', 40 );
set(gca, 'Fontsize', 28);
set( gca, 'XTickMode', 'manual', 'XTick', f_plot, 'XTickLabel', f_str);
xlim([min(f), 50000]);
ylim([-80 20]);


% Plot the A-weight Filter performance error curve
figure(3);
semilogx(f, -A_atten+A_atten2, 'k', 'LineWidth', LW, 'marker', 'p', 'markersize', 7);
hold on;
semilogx(ftype, type0_h, '--g', 'LineWidth', LW);
semilogx(ftype, type1_h, '-.m', 'LineWidth', LW);
semilogx(ftype, type2_h, ':r', 'LineWidth', LW);
semilogx(ftype, type0_l, '--g', 'LineWidth', LW);
semilogx(ftype, type1_l, '-.m', 'LineWidth', LW);
semilogx(ftype, type2_l, ':r', 'LineWidth', LW);

legend({'Actual error', 'Type-0', 'Type-1', 'Type-2'}, 'Location', 'South');
xlabel('Frequency Hz', 'Fontsize', 28);
ylabel('Absolute Error dB', 'Fontsize', 28);
title( {[title_blurb, ' A-weight Filter Error'], [num2str(Fs/1000) , ' KHz Sampling Rate']},  'Fontsize', 40 );
set(gca, 'Fontsize', 28);
set( gca, 'XTickMode', 'manual', 'XTick', f_plot, 'XTickLabel', f_str);
xlim([min(f), 50000]);
ylim(10*[-1 1]);

% Plot the C-weight Filter performance error curve
figure(4);
semilogx(f, -C_atten+C_atten2, 'k', 'LineWidth', LW, 'marker', 'd', 'markersize', 7);
hold on;
semilogx(ftype, type0_h, '--g', 'LineWidth', LW);
semilogx(ftype, type1_h, '-.m', 'LineWidth', LW);
semilogx(ftype, type2_h, ':r', 'LineWidth', LW);
semilogx(ftype, type0_l, '--g', 'LineWidth', LW);
semilogx(ftype, type1_l, '-.m', 'LineWidth', LW);
semilogx(ftype, type2_l, ':r', 'LineWidth', LW);

legend({'Actual error', 'Type-0', 'Type-1', 'Type-2'}, 'Location', 'South');
xlabel('Frequency Hz', 'Fontsize', 28);
ylabel('Absolute Error dB', 'Fontsize', 28);
title( {[title_blurb, ' C-weight Filter Error'], [num2str(Fs/1000) , ' KHz Sampling Rate']},  'Fontsize', 40 );
set(gca, 'Fontsize', 28);
set( gca, 'XTickMode', 'manual', 'XTick', f_plot, 'XTickLabel', f_str);
xlim([min(f), 50000]);
ylim(10*[-1 1]);

% Save the plots to Fig and Tiff File formats
figure(1); save_a_plot_reverb_time(2, [title_blurb, '_A_curve_', num2str(Fs/1000), 'KHZ'], 1);
figure(1); save_a_plot_reverb_time(6, [title_blurb, '_A_curve_', num2str(Fs/1000), 'KHZ'], 1);

figure(2); save_a_plot_reverb_time(2, [title_blurb, '_C_curve_', num2str(Fs/1000), 'KHZ'], 1);
figure(2); save_a_plot_reverb_time(6, [title_blurb, '_C_curve_', num2str(Fs/1000), 'KHZ'], 1);

figure(3); save_a_plot_reverb_time(2, [title_blurb, '_A_error_', num2str(Fs/1000), 'KHZ'], 1);
figure(3); save_a_plot_reverb_time(6, [title_blurb, '_A_error_', num2str(Fs/1000), 'KHZ'], 1);

figure(4); save_a_plot_reverb_time(2, [title_blurb, '_C_error_', num2str(Fs/1000), 'KHZ'], 1);
figure(4); save_a_plot_reverb_time(6, [title_blurb, '_C_error_', num2str(Fs/1000), 'KHZ'], 1);

Contact us