function varargout = ita_roomacoustics_analytic_FRF_pdi(varargin) %ITA_ROOMACOUSTICS_ANALYTIC_FRF - After Kuttruff Room Acoustics pp.66 % This function calculates the ideal FRF of an rectangular room of size % (Lx, Ly, Lz) for source position (rs_x, rs_y, rs_z) and receiver % positions ... PLUS: cartesian differentian of the source or receiver % % Syntax: % audioObjOut = ita_roomacoustics_analytic_FRF_differentiation(geometryCoordinates, sourceCoordinates, receiverCoordinates, options) % % Options (default): % 'T' (2) : description % 'f_max' (10000) : description % 'fft_degree' (18) : description % 'sourceDiff' ([0 0 0]): monopole, dipole in x ([1 0 0]), etc. % 'receiverDiff' ([0 0 0]): monopole, dipole in x ([1 0 0]), etc. % % Example: % audioObjOut = ita_roomacoustics_analytic_FRF(audioObjIn) % % See also: % ita_toolbox_gui, ita_read, ita_write, ita_generate % % Reference page in Help browser % doc ita_roomacoustics_analytic_FRF % % This file is part of the application RoomAcoustics for the ITA-Toolbox. All rights reserved. % You can find the license for this m-file in the application folder. % % Author: Pascal -- Email: pdi@akustik.rwth-aachen.de % Created: 13-Jul-2011 %% Initialization and Input Parsing sArgs = struct('pos1_data','itaCoordinates', 'pos2_data','itaCoordinates','pos3_data','itaCoordinates',... 'fft_degree',18,'T',2,'f_max',10000,'c',ita_constants('c'),'rho_0',ita_constants('rho_0'),... 'receiverDiff',[0 0 0], 'sourceDiff', [0 0 0] ); [geometry, source_pos, receiver_pos ,sArgs] = ita_parse_arguments(sArgs,varargin); %% Inits sArgs.f_max = sArgs.f_max; c = double(sArgs.c); delta_n_raw = 3*log(10)/sArgs.T; L = geometry.cart; r_source = source_pos.cart; r_receiver = receiver_pos.cart; % TODO: pdi: what about multiple receiving points? %% check loop warning off %#ok idx = 0; for nx = 0:floor(2*sArgs.f_max/c * L(1)) for ny = 0:floor( real(sqrt( (2*sArgs.f_max/c)^2 - (nx/L(1))^2 ) * L(2) )) if ~isempty(ny) && isreal(ny) idx = idx + floor( real(sqrt( (2*sArgs.f_max/c)^2 - (nx/L(1))^2 - (ny/L(2))^2) * L(3))) + 1; end end end nModes = idx; dummy = zeros(1,nModes); %% n = zeros(3,nModes); idx = 1; nx_max = floor(2*sArgs.f_max/c * L(1)); for nx = 0:nx_max ny_max = floor( real(sqrt( (2*sArgs.f_max/c)^2 - (nx/L(1))^2 ) * L(2) )); for ny = 0:ny_max nz_max = floor( real(sqrt( (2*sArgs.f_max/c)^2 - (nx/L(1))^2 - (ny/L(2))^2) * L(3))); idx_end = idx + nz_max; n(1,idx:idx_end) = nx; n(2,idx:idx_end) = ny; n(3,idx:idx_end) = 0:nz_max; idx = idx + nz_max + 1; end end ita_verbose_info([num2str(idx - nz_max - 1) ' modes have been calculated...'],1); %% Eigenfrequencies f_n = c / 2 * sqrt( sum ( bsxfun(@rdivide,n,L').^2 ,1) ); %% Source Receiver Coefficients coeff_r = potential_diff(r_receiver,L,n,sArgs.receiverDiff); coeff_s = potential_diff(r_source,L,n,sArgs.sourceDiff); %% new with SourceFactor [kg/s^2] K_n = itaValue(prod(L),'m3') * 0.5.^(sum(n > 0,1)); factor = sArgs.c^2 / K_n /itaValue('Hz2'); coeff = coeff_r .* coeff_s .* double(factor).'; %unit = []; corr = 2; % why 2??? %pdi:HUHU: not totally correct, should be j*omega after modal superposition coeff = coeff ./ corr; coeff(corr == 0) = 0; %avoid NaN for strange correction coefficient (zero) if isa(delta_n_raw,'itaSuper') delta = delta_n_raw.freq2value(f_n).'; else delta = dummy + delta_n_raw; end delta = -abs(delta); %% Set Output if nargout == 3 varargout(1) = {f_n}; varargout(2) = {coeff}; varargout(3) = {delta}; else pz = itaPZ; pz.f = f_n; pz.C = coeff; pz.sigma = delta; pz.exp_s = -1; % pdi: currently it seems to be s/m !!! Jul 2011 - this is probably due % to the douple pole handling here and the single pole handling with % mirrowing positive poles in the itaPZ.freqresp pz.unit = itaValue('1/m') * itaValue('1/m') ^ (sum(sArgs.receiverDiff) + sum(sArgs.sourceDiff)); varargout{1} = pz; end %end function end function res = potential_diff(pos,L,n,derivation) % pos: coordinates (not normalized) % L: geometry of room, length in meters % n: eigennumber % derivation: [0 0 0] - monopole, [1 0 0] - first dipole res = 1; for idx = 1:3 if mod(derivation(idx),2) % 1 3 5... res_part = sin(pi * n(idx,:) * pos(idx)/L(idx)); else % 0 2 4 res_part = cos(pi * n(idx,:) * pos(idx)/L(idx)); end res_part = res_part .* (pi * n(idx,:) / L(idx)).^derivation(idx) * (-1).^(1+mod(floor((derivation(idx)-1)/2),2)); res = res.*res_part; end end