# Auralization of High-Order Directional Sources from First-Order RIR Measurements
## Appendix B SRD ARIR  
The Matlab source code of the proposed SRD ARIR method, a source-and-receiver-directional version of the Ambisonic Spatial Decomposition Method (ASDM [22]) can be found in Listing 1. It requires the spherical harmonics implemented in [PERSON] library ([[https://github.com/polarch/Spherical-Harmonic-Transform](https://github.com/polarch/Spherical-Harmonic-Transform)]([https://github.com/polarch/Spherical-Harmonic-Transform](https://github.com/polarch/Spherical-Harmonic-Transform))). The SRD algorithm becomes ASDM when choosing \(N_{\mathrm{S}}=0\) for the source. ```
1%FirstorderB-formatRIR(ARIR)
2[x,fs]=audioread(fname);%sortedin<M,X,Y,Z>fromls1,<M,X,Y,Z>fromls2,...untills6
3%fsisthesamplingfrequency
4X=reshape(x,[size(x,1)4 size(x,2)/4]);%Nsamplesx4(5T45)x6(c.sph.1 spk.) 5
6%ParametersandSettings
7N=7;%Ambisonicorder
8 Nfft=2^ceil(log2(size(x,1)));%fftlength
9 Lsmooth_dirfluct=32;win_dirfluct=hann(Lsmooth_dirfluct);%smoothingoffluctuatingDODS/DOAs
10 Lsmooth_specdecay=2048;win_specdecay=hann(Lsmooth_specdecay);%smoothingforfs=44.1 kHz
11[Nate=700;Wadelate=1000;win_crossf=sin(pi/2+(1:Mradelate)).^2;%fadetospectrallycorrectedlateRIR]
12%WpIVDOAestimation(foromidirectionalsource)
13[b,a]=butter(4,[1002500]/(fs/2));%bandpasswithfl=100 Hzandfh=1.5 kHz
14 xbp=filterilt(b,a,sum(x,3));
15 ix=xbp(:,1).xbp(:,2);ix=circshift{fffftilt(win_dirfluct,ix),-floor(Lsmooth_dirfluct/2));
16z=xbp(:,1).xbp(:,1).xbp(:,4);iz=circshift{fffftilt(win_dirfluct,ix),-floor(Lsmooth_dirfluct/2));
17z=atan2(y,ix);zen=atan2(sqrt(ix.^2+iy.^2),iz);
18 Ydoa=getSH(N,[azi,zen],'real');%from[[https://github.com/polarch/Spherical-Harmonic-Transform/](https://github.com/polarch/Spherical-Harmonic-Transform/)]([https://github.com/polarch/Spherical-Harmonic-Transform/](https://github.com/polarch/Spherical-Harmonic-Transform/))
19
20%rf-magnitude-basedDOestimation(foromidirectionalreceiver)
21[b,a]=butter(4,[100140]/(fs/2));%bandpasswithfl=100 Hzandfh=1.4 kHz
22%wp=diffilt(b,a,squeeze(x:,1:)));
23U=1[100;0110;-100;...%6 arrayloudspeakerdirectionvectors(transposed)
2400-10;00-1;0011;
25E=xbp.^2;E=circshift{fffftilt(win_dirfluct,E),-floor(Lsmooth_dirfluct/2)};
26 rE=E*U;
27 azi=atan2(rE(:,2),rE(:,1));zen=atan2(sqrt(sum(rE(:,:1:2).^2,2)),rE(:,3));
28 Yod=getSH(N,[azi,zen],'real');%from[[https://github.com/polarch/Spherical-Harmonic-Transform/](https://github.com/polarch/Spherical-Harmonic-Transform/)]([https://github.com/polarch/Spherical-Harmonic-Transform/](https://github.com/polarch/Spherical-Harmonic-Transform/))
29
30%SRDupscaling
31X=(sum(x:,1:,3)).*Vdosa). *reshape{Yodod,[size(Yodod,1)1 size(Yodod,2)]};
32
33%SpectralDecayCorrection
3H=thirdoctave_filter_bank_linph(Nfft,fs);
34x_c=zeros(Nfft,(N+1)-2,(N+1)-2);%correctedupscaledARIR
35x_c(1:size(x,1),1,1)=x(:,1);
36 fork=1:size(H,2)
37 xthird0=iffitt(fft(x:,1,1),Nfft). *H(:,k));
38 xthird0=s=sort(circshift(fffftilt(win_specdecay,xthird0.^2),-Lsmooth_specdecay/2));
39 fornrc=1:N
40 nidx_rcv=nrcv2+(1:2=nrcv1);
41 fornsrc=1:N
42 nidx_src=nsrcv2+(1:2=nrcv4);* [46] xthirdn = ifft(fft(x(:,nidx_rcv,nidx_src),Nfft). *H(:,k)); xthirdnrms = sqrt(sum(circshift(ffttilt(win_specdecay,xthirdn(:,:),-Lsmooth_specdecay/2),2)); w.c = xthirdpms./(xthirdnrms+le-6)*sqrt((2=src+1)*(2=src+1)); % correction window x.c(:,nidx_rcv,nidx_src) = x.c(:,nidx_rcv,nidx_src)*xthirdn_sw_c;
* [47] end
* [48] end
* [49] % no spectral correction in early part: idx_xfade = Nlate+(0:Nradelate-1) x(idx_xfade,:) = (1-win_late)
* x(idx_xfade,:) + win_late. *x(idx_xfade,:,:);
* [49] x(idx_xfade(end)+1:end,:) = x(idx_xfade(end)+1:end,:); % x contains N3D normalized unpixed SRD RIR
* [50] % DFT-Domain Filter Bank function H = thirdcave_filter_bank_linph(Nfft,fs)
* f=linspace(0,fs/2,Nfft/2+1);
* f(1)=f(2)/4;
* fc = 25*2. *(0:1/3:9.9); %third-octave vector H = zeros(Nfft/2+1,length(fc));
* for k = 1:length(fc) ntnctwaves = log2(f/(fc)(k)+3; % 3 rd-octave distance from center freq.