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

[Wa, Wc]=AC_weight_NB(f, output_units, one_or_two_sided)
function [Wa, Wc]=AC_weight_NB(f, output_units, one_or_two_sided)
% % AC_weight_NB: Calculates the A and C frequency weights at specified frequencies 
% % 
% % Syntax;  
% % 
% % [Wa, Wc]=AC_weight_NB(f, flag);
% % 
% % ***********************************************************
% % 
% % Description
% %
% % This program calculates the attentuation factors for the A-weighting
% % and C-weighting frequency filters given the frequencies of interest.  
% %
% % [Wa]=AC_weight_NB(f);
% % 
% % Returns Wa (dB) a numeric array of A-weighting frequency filter 
% % attenuation factors at each frequency in the input vector f (Hz). 
% % Wa is the vector of attenuation factor in dB and f is the input 
% % vector of frequencies in (Hz).  
% % 
% % [Wa, Wc]=AC_weight_NB(f);
% % Returns Wc (dB) a numeric array of C-weighting frequency filter 
% % attenuation factors.  
% % 
% % [Wa, Wc]=AC_weight_NB(f, output_units);
% % output_units=0; Returns Wa and Wc in (dB).
% % output_units=1; Returns Wa and Wc as attenuation ratios.  
% % 
% % [Wa, Wc]=AC_weight_NB(f, output_units, one_or_two_sided);  
% % one_or_two_sided=0; Returns Wa and Wc in the same frequencies as f.
% % one_or_two_sided=1; Returns Wa and Wc in a two sided frequency array 
% % f is symmetrically copied about the nyquist frequency assuming that 
% % the nyquist frequency is max(f)/2. 
% % 
% % if one_or_two_sided=1; then the input vector f becomes 
% % f = [f(1:floor(end/2)) fliplr(f(1:ceil(end/2)))];
% % 
% % 
% % ***********************************************************
% % 
% % Input Variables
% % 
% % f=200:     % (Hz) is a numeric array of frequencies to calculate the 
% %            % A and C attenuation factors of the frequency weighting 
% %            %         filter values.  
% %            % 
% %            % default is the third octave bands from 20 Hz to 20 KHz
% % 
% % output_units=0; % is a scalar boolean which specifies the units of the 
% %                 % attenuation factor of the frequency of the frequency weighting 
% %                 % filter values.  
% %                 % 0 is for dB and 1 is for ratio.  
% %                 % 
% %                 % default is output_units=0;
% % 
% % one_or_two_sided=0; % is a scalar boolean which specifies whether to
% %                     % use f or to modify f by making it symmetric about 
% %                     % the nyquist frequency. 
% %                     % 0 just uses f as it is input.
% %                     % 1 modifies f by the formula
% %                     % f = [f(1:floor(end/2)) fliplr(f(1:ceil(end/2)))];
% %                     % default is one_or_two_sided=0;
% % 
% % ***********************************************************
% % 
% % Output Variables
% % 
% % Wa is a numeric array of A-weighting attenuation factors of the  
% %            frequency filter values in specified units at the 
% %            specified frequencies.  
% % 
% % Wc is a numeric array of C-weighting attenuation factors of the  
% %            frequency filter values in specified units at the 
% %            specified frequencies.  
% % 
% % ***********************************************************
% 
%
% Example='1';
%
% % The audiometric frequency range is 20 Hz to 20000 Hz.
% % Make a plot of the A and C weighting attenuation factors.
%
% f=20:1:20000;
% [Wa, Wc]=AC_weight_NB(f);
% 
% figure(1);
% semilogx(f, Wa, 'r');
% hold on;
% semilogx(f, Wc, 'b');
% legend('A-Weighting', 'C-Weighting', 'Location', 'southeast');
% xlabel('Frequency (Hz)', 'Fontsize', 18);
% ylabel('Attenuation Factors (dB ref. 20\muPa)', 'Fontsize', 18);
% title('Frequency Weighting Filter Functions for Sound', 'Fontsize', 20);
% set(gca, 'Fontsize', 14);
% 
% 
% 
% Example='2';
%
% % The frequency range for impulsive noise can reach up to 1 MHz .
% % Make a plot of the A and C weighting attenuation factors.
% % Using 1/12 octave bands from 20 HZ to 1 MHz.
%  
% N=12;        % twelve equal divisions per octave
% min_f=20;    % minimum frequency
% max_f=10^6;  % maximum frequency
% 
% [f]=nth_freq_band(N, min_f, max_f);
% 
% [Wa, Wc]=AC_weight_NB(f);
% 
% figure(2);
% semilogx(f, Wa, 'r');
% hold on;
% semilogx(f, Wc, 'b');
% legend('A-Weighting', 'C-Weighting', 'Location', 'northeast');
% xlabel('Frequency (Hz)', 'Fontsize', 18);
% ylabel('Attenuation Factors (dB ref. 20\muPa)', 'Fontsize', 18);
% title('Frequency Weighting Filter Functions for Sound', 'Fontsize', 20);
% set(gca, 'Fontsize', 14);
% 
% 
% 
% Example='3';
%
% % The frequency range for a narrow band fft at a sampling rate of 50 KHz 
% % is 1 to 25000 Hz.  
% % 
% % Make a plot of the A and C weighting attenuation factors for 
% % simulated data.
%  
% f=1:1:25000;
% % Use the ratio method
% [Wa, Wc]=AC_weight_NB(f, 1);
% 
% Fs=50000;
% x=randn(100000, 1);
% bin_size=50000;
% num_averages=100;
% [SP]=pressure_spectra(x, Fs, bin_size, num_averages );
% 
% figure(3);
% Pref=20*10^(-6);
% semilogx(f, 20*log10(SP./Pref), 'k');
% hold on;
% semilogx(f, 20*log10(SP.*Wa./Pref), 'r');
% semilogx(f, 20*log10(SP.*Wc./Pref), 'b');
% legend('Linear-Weighting', 'A-Weighting', 'C-Weighting', 'Location', 'southeast');
% xlabel('Frequency (Hz)', 'Fontsize', 18);
% ylabel('Sound Pressure (dB ref. 20\muPa)', 'Fontsize', 18);
% title('Frequency Weighting Sound Spectra', 'Fontsize', 20);
% set(gca, 'Fontsize', 14);
% 
% 
% 
% Example='4';
%
% % The frequency range for a narrow band fft at a sampling rate of 50 KHz 
% % is 1 to 25000 Hz.  
% % 
% % Make a plot of the A and C weighting attenuation factors for 
% % simulated data.
% % Using 1/3 octave bands from 20 HZ to 20 KHz.
%  
% N=3;          % three equal divisions per octave
% min_f=20;     % minimum frequency
% max_f=20000;  % maximum frequency
% 
% [f]=nth_freq_band(N, min_f, max_f);
% % Use the dB method
% [Wa, Wc]=AC_weight_NB(f, 0);
% 
% Fs=50000;
% x=randn(100000, 1);
% num_x_filter=2;
% [fc_out, SP_levels, SP_peak_levels, SP_bands]=Nth_oct_time_filter2(x, Fs, num_x_filter, N, f);
% 
% 
% figure(4);
% semilogx(f', SP_levels, 'k');
% hold on;
% semilogx(f', SP_levels+Wa', 'r');
% semilogx(f', SP_levels+Wc', 'b');
% legend('Linear-Weighting', 'A-Weighting', 'C-Weighting', 'Location', 'southeast');
% xlabel('Frequency (Hz)', 'Fontsize', 18);
% ylabel('Sound Pressure (dB ref. 20\muPa)', 'Fontsize', 18);
% title('Frequency Weighting Sound Spectra', 'Fontsize', 20);
% set(gca, 'Fontsize', 14);
% 
% 
% 
% 
% Example='5';
%
% % The frequency range for a narrow band fft at a sampling rate of 50 KHz 
% % is 1 to 25000 Hz.  
% % 
% % Make a plot of the A and C weighting attenuation factors for 
% % simulated data.
% % Using 1/3 octave bands from 20 HZ to 20 KHz.
%  
% N=3;          % three equal divisions per octave
% min_f=20;     % minimum frequency
% max_f=20000;  % maximum frequency
% 
% [f]=nth_freq_band(N, min_f, max_f);
% % Use the dB method
% [Wa, Wc]=AC_weight_NB(f, 0);
% 
% Fs=50000;
% x = spatialPattern([100000,1], -1);
% num_x_filter=2;
% [fc_out, SP_levels, SP_peak_levels, SP_bands]=Nth_oct_time_filter2(x, Fs, num_x_filter, N, f);
% 
% 
% figure(5);
% 
% semilogx(f', SP_levels, 'k');
% hold on;
% semilogx(f', SP_levels+Wa', 'r');
% semilogx(f', SP_levels+Wc', 'b');
% legend('Linear-Weighting', 'A-Weighting', 'C-Weighting', 'Location', 'southeast');
% xlabel('Frequency (Hz)', 'Fontsize', 18);
% ylabel('Sound Pressure (dB ref. 20\muPa)', 'Fontsize', 18);
% title('Frequency Weighting Sound Spectra of Pink Noise', 'Fontsize', 20);
% set(gca, 'Fontsize', 14);
% 
% 
% 
% % ***********************************************************
% % 
% % References
% % 
% % Reference American National Standard for Speicification for Sound Level
% % Meters ANSI S1.4-1983.  
% %    Table IV. Random incidence relative response level as a function of
% %              frequency for various weightings
% %    Appendix C: relative Response of Frequency Weighting Characteristic
% % 
% % ***********************************************************
% % 
% % 
% % Subprograms
% %
% % 
% % List of Dependent Subprograms for 
% % AC_weight_NB
% % 
% % 
% % Program Name   Author   FEX ID#
% %  1) convert_double		
% %  2) estimatenoise		John D'Errico		16683	
% %  3) filter_settling_data		
% %  4) flat_top		
% %  5) moving		Aslak Grinsted		8251	
% %  6) nth_freq_band		
% %  7) Nth_oct_time_filter2		
% %  8) Nth_octdsgn		
% %  9) number_of_averages		
% % 10) pressure_spectra		
% % 11) sd_round		
% % 12) spatialPattern		Jon Yearsley		5091	
% % 13) spectra_estimate		
% % 14) sub_mean		
% % 15) window_correction_factor		
% % 16) wsmooth		Damien Garcia		NA	
% % 
% % 
% % ***********************************************************
% % 
% % Written by William J. Murphy, October 3, 2003
% % 
% % modified by Edward Zechmann 
% % 
% % modified 23 September   2008    Updated Comments
% % 
% % modified  7 January     2008    Updated Comments
% % 
% % ***********************************************************
% % 
% % Please feel free to modify this code.
% % 
% % See also: Leq_all_calc, ACdsgn, adsgn, cdsgn, Aweight_time_filter, Cweight_time_filter, resample
% % 


% This code is added to help find the dependent functions to run the
% examples
void_call=0;
if void_call
    pressure_spectra(0);
    spatialPattern(0);
    nth_freq_band(0);
    Nth_oct_time_filter2(0);
end


if (nargin < 1 || isempty(f)) || ~isnumeric(f)
     f=[20 25 31 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];
end

if (nargin < 2 || isempty(output_units)) || ~isnumeric(output_units) 
    output_units=0;
end

if ~isequal(output_units, 1)
    output_units=0;
end
    
if (nargin < 3 || isempty(one_or_two_sided)) || ~isnumeric(one_or_two_sided) 
    one_or_two_sided=0;
end

if ~isequal(one_or_two_sided, 1)
    one_or_two_sided=0;
end


% Define the filter constants

K1 =     2.242881e16;
K2 =     1.025119;
K3 =     1.562339;
f1 =    20.598997;
f2 =   107.65265;
f3 =   737.86223;
f4 = 12194.22;
f5 =   158.48932;

% Extend the frequency vector to double sided if necessary
if isequal(one_or_two_sided, 1)
    [m1, n1]=size(f);
    
    if m1 > n1
        f=f';
    end
    
    f = [f(1:floor(end/2)) fliplr(f(1:ceil(end/2)))];
end


switch output_units
    
    case 0
        
        Wc = 10.*log10((K1.*f.^4)./((f.^2 + f1.^2).^2 .*(f.^2 +f4.^2).^2));
        Wa = 10.*log10((K3.*f.^4)./((f.^2 + f2.^2).*(f.^2 + f3.^2))) + Wc;
        
    case 1
        
        Wc = (K1.*f.^4)./((f.^2 + f1.^2).^2 .*(f.^2 +f4.^2).^2);
        Wa = (K3.*f.^4)./((f.^2 + f2.^2).*(f.^2 + f3.^2)).*Wc;
        
    otherwise
        
        Wc = 10.*log10((K1.*f.^4)./((f.^2 + f1.^2).^2 .*(f.^2 +f4.^2).^2));
        Wa = 10.*log10((K3.*f.^4)./((f.^2 + f2.^2).*(f.^2 + f3.^2))) + Wc; 
        
end

Contact us