ita_sound_power.m 2.92 KB
Newer Older
1 2 3 4
function varargout = ita_sound_power(varargin)
%ITA_SOUND_POWER - calculate the sound power of a source in a specific room
%  This function takes spl-data and RT of the empty room to calculate the
%  equivalent absorption area of the tested room. This is used to calculate
5
%  the sound power of the source according to ISO 3741.
6 7 8 9 10 11 12 13 14 15 16
%
%  Syntax:
%   audioObjOut = ita_sound_power(audioObjIn, T_empty, options)
%
%   Options (default):
%           'room_volume' (124) :  default value fits ITA reverberation
%                                  chamber
%           'room_surface' (181):  default value fits ITA reverberation
%                                  chamber 
%           'T' (20)            :  Temperature in deg Celsius
%           'RH' (0.5)          :  Relative Humidity
17
%           'p_s' (101325)      :  Static pressure in Pa
18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
%
%  Example:
%   audioObjOut = ita_sound_power(audioObjIn)
%
%  See also:
%   ita_sabine
%
%   Reference page in Help browser 
%        <a href="matlab:doc ita_sound_power">doc ita_sound_power</a>

% <ITA-Toolbox>
% This file is part of the ITA-Toolbox. Some rights reserved. 
% You can find the license for this m-file in the license.txt file in the ITA-Toolbox folder. 
% </ITA-Toolbox>

% Author: Christian Haar -- Email: christian.haar@akustik.rwth-aachen.de
% Created:  24-Jun-2010 
35 36
% Re-write to conform with ISO 3741: Markus Mueller-Trapet -- Email: markus.mueller-trapet@nrc.ca
% Modified: 28-Nov-2017
37 38

%% Initialization and Input Parsing
39
sArgs        = struct('pos1_spl', 'itaSuper', 'pos2_T_empty', 'itaSuper', 'room_volume', 124,'room_surface',181,'T',20,'RH',0.5,'p_s',101325, 'freqRange', ita_preferences('freqRange'), 'bandsPerOctave', ita_preferences('bandsPerOctave'), 'bandMode','fft');
40 41
[spl,T_empty,sArgs] = ita_parse_arguments(sArgs,varargin); 

42 43
% acoustic constants
c = 20.05*sqrt(273 + sArgs.T);
44
Z_0 = itaValue(400,'kg/(m^2*s)'); % this is the fixed value to get the reference for 1pW right
45
% correction factors for the meteorological conditions (not in dB!)
46 47
C1 = 101325./double(sArgs.p_s).*sqrt((273.15 + sArgs.T)/314);
C2 = 101325./double(sArgs.p_s).*sqrt((273.15 + sArgs.T)/296).^3;
48

49
%% calculate sound pressure level data
Markus Mueller-Trapet's avatar
Markus Mueller-Trapet committed
50
% first third-octaves then average
51
spl = ita_spk2frequencybands(spl, 'freqRange',sArgs.freqRange , 'bandsPerOctave',sArgs.bandsPerOctave,'mode',sArgs.bandMode);
Markus Mueller-Trapet's avatar
Markus Mueller-Trapet committed
52
spl_m = sqrt(mean(spl^2));
53 54 55 56 57 58
spl_m = itaResult(spl_m,T_empty.freqVector);

%% calculate equivalent absorption area
A = 55.26*itaValue(double(sArgs.room_volume)/c,'s*m^2')/T_empty;

%% calculate sound power (Eq 20 in ISO 3741)
59 60 61
sound_power = spl_m^2 * A /(4*Z_0);
% apply exponent and frequency-dependent corrections
sound_power.freq =  C1.*C2.*sound_power.freq.*exp(A.freq./double(sArgs.room_surface)).*(1 + double(sArgs.room_surface)*c./(8*double(sArgs.room_volume).*sound_power.freqVector));
62 63 64 65 66 67 68 69 70

%% Add history line
sound_power = ita_metainfo_add_historyline(sound_power,mfilename,varargin);

%% Set Output
varargout(1) = {sound_power}; 

%end function
end