Aufgrund einer Wartung wird GitLab am 28.09. zwischen 10:00 und 11:00 Uhr kurzzeitig nicht zur Verfügung stehen. / Due to maintenance, GitLab will be temporarily unavailable on 28.09. between 10:00 and 11:00 am.

ita_diffusion_coefficient.m 3.49 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
function varargout = ita_diffusion_coefficient(varargin)
%ITA_DIFFUSION_COEFFICIENT - diffusion coefficient according to AES-4id-2001
%  This function calculates the diffusion coefficient as documented in the
%  AES-4id-2001 document.
%  Input arguments are the reflection directivity of the test sample and
%  the microphone array positions.
%
%  If an optional third argument is given, the direction of the different
%  plane waves are expected. They will be used to determine the
%  random-incidence value using Paris' formula.
%
%  The output is the diffusion coefficient for each plane wave direction.
%  If a second output argument is requested, the random-incidence value is
%  also returned.
%
%  Syntax:
%   audioObjOut = ita_diffusion_coefficient(p_sample,mics,'directions',plane_waves)
%
%   Options (default):
%           'directions' ([]) : direction of the plane waves
%
%  See also:
%   ita_toolbox_gui, ita_read, ita_write, ita_generate
%
%   Reference page in Help browser 
%        <a href="matlab:doc ita_diffusion_coefficient">doc ita_diffusion_coefficient</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: Markus Mueller Trapet -- Email: mmt@akustik.rwth-aachen.de
% Created:  27-Mar-2011 


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

%% Initialization and Input Parsing
sArgs        = struct('pos1_data','itaSuper', 'pos2_mics', 'itaCoordinates','directions',[]);
[data,mics,sArgs] = ita_parse_arguments(sArgs,varargin);

%% rearrange data
if numel(data) > 1
    nDirections = numel(data);
    p_sample = data;
else
    nDirections = data.nChannels/mics.nPoints;
    if round(nDirections) ~= nDirections
        error([thisFuncStr 'strange number of plane waves, what did you give me?']);
    end
    p_sample = repmat(data.ch(1:mics.nPoints),[nDirections,1]);
    for i = 1:nDirections
        p_sample(i) = data.ch((i-1)*mics.nPoints+(1:mics.nPoints));
    end
end

clear data;

%% diffusion coefficient by autocorrelation
d_tmp = zeros(p_sample(1).nBins,nDirections);
% N = mics.nPoints;
if isempty(mics.weights) || all(mics.weights == 1)
    theta   = mics.theta;
    phi     = mics.phi;
    dtheta  = mean(diff(unique(theta)));
    dphi    = mean(diff(unique(phi)));
    weights = dtheta*dphi.*sin(theta(:));
else
    weights = mics.weights;
end
weights = weights./min(weights);

% calculate including weights according to standard
for iWave = 1:nDirections
    E = abs(p_sample(iWave).freq).^2;
    d_tmp(:,iWave) = (sum(bsxfun(@times,E,weights(:).'),2).^2 - sum(bsxfun(@times,E.^2,weights(:).'),2))./((sum(weights)-1).*sum(bsxfun(@times,E.^2,weights(:).'),2));
end

d = p_sample(1).ch(1);
d.freq = d_tmp;
d.channelNames = cellstr([repmat('Direction ',[nDirections,1]) num2str((1:nDirections).','%02d')]);

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

%% Set Output
varargout(1) = {d};
if nargout > 1
    if ~isempty(sArgs.directions) && isa(sArgs.directions,'itaCoordinates')
        d_weights = sin(2*sArgs.directions.theta);
        d_weights = itaValue(d_weights./sum(d_weights));
        varargout(2) = {sum(d*d_weights)};
    else
        ita_verbose_info([thisFuncStr 'I do not have any info on the incoming plane waves, just using mean'],0);
        varargout(2) = {mean(d)};
    end
    if nargout > 2
        varargout(3) = {p_sample};
    end
end

%end function
end