%ITA_SPH_MODAL_STRENGTH - Calculate the modal strength function for a spherical array% This function calculates the modal strength for a plane wave incident % onto a spherical microphone aray. The output is a diagonal matrix if the sampling % has a unique radius. If not, the output is a matrix with full rank. In case the wave % number k is given as a vector the output is of size N x M x nBins.

%% See also:% ita_toolbox_gui, ita_read, ita_write, ita_generate%% Reference page in Help browser % <a href="matlab:doc ita_sph_modal_strength">doc ita_sph_modal_strength</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>

% Created: 26-Feb-2016 %% Initialization and Input ParsingsArgs=struct('pos1_sampling','itaCoordinates',...'pos2_maxOrder','integer',...'pos3_k','double',...'pos4_type','string',...'transducer','microphone',...'scatterer',[],...'hankelKind',2,...'dist',2);[sampling,Nmax,kVec,type,sArgs]=ita_parse_arguments(sArgs,varargin);% check if the sampling has a unique radiusifnumel(unique(sampling.r))>1uniRSort=sort(unique(sampling.r));ifuniRSort(1)>uniRSort(end)*(1-2*eps)&&uniRSort(end)<uniRSort(1)*(1+2*eps)sampling.r=repmat(mean(uniRSort),size(sampling.r));uniqueRad=mean(uniRSort);B=zeros((Nmax+1)^2,(Nmax+1)^2,numel(kVec));else[uniqueRad,~,idxSamplingPos]=unique(sampling.r);B=zeros(sampling.nPoints,(Nmax+1)^2,numel(kVec));endelseuniqueRad=unique(sampling.r);end% calculate impedance only once to avoid dealing with unitsif(strcmp(sArgs.transducer,'ls'))Z_air=double(ita_constants('z_0'));endn=(0:Nmax).';foridxRad=1:numel(uniqueRad)switchtypecase{'open','interior'}bn=ita_sph_besselj(n,uniqueRad(idxRad)*kVec);case'rigid'ifisempty(sArgs.scatterer)% if the scatterer and the sampling share the same radius% use the wronskian for speed reasonsbn=1./ita_sph_besselh_diff(n,sArgs.hankelKind,uniqueRad(idxRad)*kVec);else% if a spherical scatterer with a radius different than the% radius of the sampling is used we cannot use the% wronskian relationbn=ita_sph_besselj(n,uniqueRad(idxRad)*kVec)-...ita_sph_besselj_diff(n,sArgs.scatterer*kVec)./...ita_sph_besselh_diff(n,sArgs.hankelKind,sArgs.scatterer*kVec).*...ita_sph_besselh(n,sArgs.hankelKind,uniqueRad(idxRad)*kVec);endcase'cardioid'bn=ita_sph_besselj(n,uniqueRad(idxRad).*kVec)-...1i.*ita_sph_besselj_diff(n,1,uniqueRad(idxRad).*kVec);otherwiseita_verbose_info('This is not a valid design type.',0);varargout{1}=[];returnendswitchsArgs.transducercase{'microphone','mic'}ifstrcmp(type,'rigid')&&isempty(sArgs.scatterer)

elsebn=diag(4*pi*1i.^n)*bn;endcase{'loudspeaker','ls'}ifstrcmp(type,'rigid')&&isempty(sArgs.scatterer)bn=diag(-1i.^n*Z_air)*bn.*ita_sph_besselh(n,sArgs.hankelKind,sArgs.dist*kVec);elseita_verbose_info('This design type makes no sense for a loudspeaker array.',0);varargout{1}=[];returnendotherwiseita_verbose_info('I do not know this kind of transducer.',0);varargout{1}=[];returnendifkVec(1)==0bn(:,1)=real(bn(:,2));endbn=ita_sph_eye(Nmax,'n-nm').'*bn;ifnumel(uniqueRad)>1B(idxSamplingPos==idxRad,:,:)=repmat(permute(bn,[3,1,2]),[numel(idxSamplingPos(idxSamplingPos==idxRad)),1,1]);elseforidxFreq=1:numel(kVec)B(:,:,idxFreq)=diag(bn(:,idxFreq));endendendvarargout{1}=B;end