ita_scattering_coefficient_diffuse.m 7.4 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212
function varargout = ita_scattering_coefficient_diffuse(varargin)
%ITA_SCATTERING_COEFFICIENT_DIFFUSE - random-incidence scattering coefficient
% This function calculate the scattering coefficient from the reverberation
% times calculated in four different conditions:
%
%                Test sample         |   Turning table
% - RT1     |        NO              |      NO
% - RT2     |        YES             |      NO
% - RT3     |        NO              |      YES
% - RT4     |        YES             |      YES
%
%  atm_mat is a matrix that contains information about the atmospheric
%  condition during each of the four measurements. It must be structured as
%  follow:
%
%  atm_mat = [T1  T2  T3  T4;
%             RH1 RH2 RH3 RH4]
%
%  where T stands for temperature (in degree Celsius) and RH stands for
%  relative humidity in percent
%
%  RT contains the performed measurements and is structured like this:
%
%  RT = [mean([RT1_p1, RT1_p2]);
%        mean([RT2_p1, RT2_p2]);
%        mean([RT3_p1, RT3_p2]);
%        mean([RT4_p1, RT4_p2])];
%
%  where RTX_pY is the reverberation time evaluated at position Y
%
%
%  Syntax:
%   s = ita_scattering_coefficient_diffuse(RT,atm_vec)
%
%   Options (default):
%           'plot' (false)         : plot all results
%           'alphaMode' ('Eyring') : Eyring or Sabine for alpha
%           'scaleFactor' (5)      : small-scale factor
%           'sampleArea' (pi*0.4^2): surface area of sample
%           'RT_std' (itaResult()) : std dev of RT
%           'S_room' (9.05)        : surface area of empty room
%           'V_room' (1.67)        : volume of empty room
%
%
%   Reference page in Help browser
%        <a href="matlab:doc ita_scattering_coefficient_diffuse">doc ita_scattering_coefficient_diffuse</a>

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


% Author: MMT -- Email: mmt@akustik.rwth-aachen.de
% Created:  21-Jul-2010

%% Get Function String
thisFuncStr  = [upper(mfilename) ':']; %Use to show warnings or infos in this functions

%% Initialization and Input Parsing
sArgs        = struct('pos1_RT','itaResult', 'pos2_atm_mat','double','plot',false,'alphaMode','Eyring','sampleArea',pi*(0.4^2),'scaleFactor',5,'RT_std',itaResult(),'S_room',9.05,'V_room',1.67);
[RT,atm_mat,sArgs] = ita_parse_arguments(sArgs,varargin);

channelNames = {'no sample; no rotation','with sample; no rotation',...
    'no sample; with rotation','with sample; with roation'};
if numel(RT) > 1
    RT = merge(RT);
end
if RT.nChannels ~= 4
    error([thisFuncStr '4 reverberation times are required for this parameter']);
end
RT.comment = 'Reverberation Time';
RT.channelNames = channelNames;

%% Frequencies
scale_factor = sArgs.scaleFactor;
% Frequencies under analysis
f = RT(1).freqVector;
f0 = f./scale_factor;

%% Constants
c = cell(4,1);
m = cell(4,1);
for i = 1:4
    [c{i},m{i}] = ita_constants({'c','m'},'T',atm_mat(1,i),'phi',atm_mat(2,i)/100,'f',f);
end

S_baseplate = pi*(0.43^2);  % Area of the baseplate in square meters
S_sample    = sArgs.sampleArea; % Area of the test sample in square meters
S_room      = sArgs.S_room;
V_room      = sArgs.V_room;

%% Calculate the virtual absorption coefficient in each of the four cases
% empty - no rotation
alpha(1) = ita_sabine('c',double(c{1}),'v',V_room,'s',S_room,'m',double(m{1}),'t60',RT.ch(1),'mode',sArgs.alphaMode);
% with sample - no rotation
alpha(2) = ita_sabine('c',double(c{2}),'v',V_room,'s',S_room,'m',double(m{2}),'t60',RT.ch(2),'mode',sArgs.alphaMode);
% empty - with rotation
alpha(3) = ita_sabine('c',double(c{3}),'v',V_room,'s',S_room,'m',double(m{3}),'t60',RT.ch(3),'mode',sArgs.alphaMode);
% with sample - with rotation
alpha(4) = ita_sabine('c',double(c{4}),'v',V_room,'s',S_room,'m',double(m{4}),'t60',RT.ch(4),'mode',sArgs.alphaMode);

RT.freqVector = f0;
RT.channelNames = channelNames;

alpha = merge(alpha);
alpha.freqVector = f0;
alpha.comment = ['Absorption Coefficient (according to ' sArgs.alphaMode ')'];
alpha.channelNames = channelNames;
A_empty = itaValue(S_room,'m^2') * alpha.ch(1);
A_empty.comment = ['Equivalent Absorption Area (according to ' sArgs.alphaMode ')'];

%% Calculations
% calculate the random incidence absorption coefficient alpha_s
alpha_s = S_room/S_sample * (alpha.ch(2) - alpha.ch(1)) + alpha.ch(1);

% calculate the specular absorption coefficient alpha_spec
alpha_spec = S_room/S_sample * (alpha.ch(4) - alpha.ch(3)) + alpha.ch(3);

% Calculate the scattering coefficient of the sample
s = 1 - (1 - alpha_spec)/(1 - alpha_s);
s.freq = max(s.freq,0);

% Calculate the scattering coefficient of the base plate
s_baseplate = S_room/S_baseplate * (alpha.ch(3) - alpha.ch(1));
s_baseplate.freq = max(s_baseplate.freq,0);

%% Meta Infos
s_baseplate.channelNames = {'scattering coefficient of the base plate'};
s_baseplate.comment = 'Scattering Coefficient';
s.comment = 'Scattering Coefficient';
s.channelNames = {'scattering coefficient of the sample'};

%% determine accuracy with error propagation
if ~isempty(sArgs.RT_std)
    stdRT  = merge(sArgs.RT_std);
    stdRT.freqVector = f0;
    
    K = itaValue(24.*log(10).*V_room./mean(cellfun(@double,c))./S_sample,'s');
    stdS1 = sum((merge([(1-s) (1-s) 0*s+1 0*s+1])*stdRT/RT^2)^2);
    
    r0 = 0*s; % zero ...
    r1 = 1 + r0; % and one dummy coefficients
    r13 = r1;
    r24 = r1;
    r24.freq = max(0,1 - alpha_spec.freq);
    
    % error propagation with correlation
    stdS2 = (s-1)*(r13*stdRT.ch(1)*stdRT.ch(3)/(RT.ch(1)*RT.ch(3))^2 + r24*stdRT.ch(2)*stdRT.ch(4)/(RT.ch(2)*RT.ch(4))^2);
    
    stdS = K/(1-alpha_s)*sqrt(stdS1 + 2*stdS2);
    stdS.channelNames = {'Standard deviation of the scattering coefficient'};
else
    stdRT = itaResult();
    stdS = itaResult();
end

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

%% Set Output
varargout(1) = {s};
if nargout > 1
    varargout(2) = {s_baseplate};
    if nargout > 2
        varargout(3) = {alpha};
        if nargout > 3
            varargout(4) = {RT};
            if nargout > 4
                varargout(5) = {stdRT};
                if nargout > 5
                    varargout(6) = {stdS};
                end
            end
        end
    end
end

%% do plotting
if sArgs.plot
    % Max ISO limits
    ISO_limits_354 = ita_ISO_limits_absorption(V_room*(scale_factor^3),'freqVector',f0)/scale_factor^2;
    ISO_limits_17497 = ita_ISO_limits_scattering(V_room,'freqVector',f0);
    alpha_s_max = ISO_limits_17497(1);
    A_empty_ISO = merge(ISO_limits_354,ISO_limits_17497(2));
    s_baseplate_ISO = ISO_limits_17497(3);
    
    % reverberation times
    ita_plot_freq(RT,'ylim',[0 3]);
    ylabel('Reverberation Time [s]');
    
    % absorption coefficient of the sample alpha
    ita_plot_freq(merge(alpha_s,alpha_s_max),'ylim',[0 1.05]);
    
    % maximum alowable absorption area of the emtpy room
    ita_plot_freq(merge(A_empty,A_empty_ISO),'ylim',[0 max(1.05,ceil(max(A_empty_ISO.freq)))]);
    ylabel('Absorption Area [m^2]');
    
    % maximum allowable scattering coefficient of the baseplate
    ita_plot_freq(merge(s_baseplate,s_baseplate_ISO),'ylim',[0 0.5]);
    title('Scattering Coefficient - Base Plate');
    ylabel('Random Incidence Scattering Coefficient [-]');
    
    % scattering coefficient of the sample
    ita_plot_freq(s,'ylim',[0 1.05]);
    title('Scattering Coefficient - Sample');
    ylabel('Random Incidence Scattering Coefficient [-]');
    
end

%end function
end