% File packer for Dario D'Orazio anechoic recordings 
% into Ambix-16 multich WAV format
% (C) Angelo Farina, 01 July 2022

% the original recordings were in separate files, each is stereo, each 
% containing twice the same microphone signal repeated in the two channels
% there are 12 of these files in each subfolder
% this script removes the duplicate channels and encodes the signals
% into Ambix 3rd order (16 channels), then it saves a single 16-channels
% WAV file

clear all;
% compute the Ambix encoding coefficients taking into account the spherical
% coordinates of the 12 microphones surrounding the sound source as
% described in this web page: http://www.angelofarina.it/SPS-conversion.htm
% considering the doechaedron format and the Ambisonics encoding formulas
% contained in this web page: http://www.angelofarina.it/Aurora/HOA_explicit_formulas.htm

% first define the 12 values of azimut (az) and elevation (el)
az= [-120,0,120,120,60,0,-60,-120,180,-60,180,60];
el= [52.6,52.6,52.6,-10.8,10.8,-10.8,10.8,-10.8,10.8,-52.6,-52.6,-52.6];

% now define the encoding coefficients for the 16 channels and the 12
% microphones
K=zeros(12,16);
for d=1:12
    a=az(d)/180*pi; e=el(d)/180*pi;
    K(d,1)=1;
    K(d,2)=sin(a)*cos(e);
    K(d,3)=sin(e);
    K(d,4)=cos(a)*cos(e);
    K(d,5)=sqrt(3/4)*sin(2*a)*(cos(e))^2;
    K(d,6)=sqrt(3/4)*sin(a)*sin(2*e);
    K(d,7)=(1/2)*(3*(sin(e))^2-1);
    K(d,8)=sqrt(3/4)*cos(a)*sin(2*e);
    K(d,9)=sqrt(3/4)*cos(2*a)*(cos(e))^2;
    K(d,10)=sqrt(5/8)*sin(3*a)*(cos(e))^3;
    K(d,11)=sqrt(15/4)*sin(2*a)*sin(e)*(cos(e))^2;
    K(d,12)=sqrt(3/8)*sin(a)*cos(e)*(5*(sin(e))^2-1);
    K(d,13)=(1/2)*sin(e)*(5*(sin(e))^2-3);
    K(d,14)=sqrt(3/8)*cos(a)*cos(e)*(5*(sin(e))^2-1);
    K(d,15)=sqrt(15/4)*cos(2*a)*sin(e)*(cos(e))^2;
    K(d,16)=sqrt(5/8)*cos(3*a)*(cos(e))^3;
end

rootdir = 'Puccini';
filelist = dir(fullfile(rootdir, '**\*.*'));  %get list of files and folders in any subfolder
filelist = filelist(~[filelist.isdir]);  %remove folders from list
ssz=size(filelist);
sz =ssz(1);
disp(strcat('Number of individual files found: ', string(sz)));
for i=1:sz
    fn=filelist(i).name;
    if(size(fn,2)>6)
        ffn=fn(end-6:end);
    else
        ffn='';
    end
    if strcmp(ffn,'001.wav')
        disp(strcat('Processing file n. ', string(i),' - name:',' ',fn));
        fname=fn(1:end-7);
        path=filelist(i).folder;
        ffname= fullfile(path,fname);
        [w01,fs]=audioread(strcat(ffname,'001.wav'));
        [w02,fs]=audioread(strcat(ffname,'002.wav'));
        [w03,fs]=audioread(strcat(ffname,'003.wav'));
        [w04,fs]=audioread(strcat(ffname,'004.wav'));
        [w05,fs]=audioread(strcat(ffname,'005.wav'));
        [w06,fs]=audioread(strcat(ffname,'006.wav'));
        [w07,fs]=audioread(strcat(ffname,'007.wav'));
        [w08,fs]=audioread(strcat(ffname,'008.wav'));
        [w09,fs]=audioread(strcat(ffname,'009.wav'));
        [w10,fs]=audioread(strcat(ffname,'010.wav'));
        [w11,fs]=audioread(strcat(ffname,'011.wav'));
        [w12,fs]=audioread(strcat(ffname,'012.wav'));
        y=[w01(:,1),w02(:,1),w03(:,1),w04(:,1),w05(:,1),w06(:,1),w07(:,1),w08(:,1),w09(:,1),w10(:,1),w11(:,1),w12(:,1)];
        ln=size(y,1);
        z=zeros(ln,16);
        for d=1:12
            for c=1:16
                 z(:,c)=z(:,c)+y(:,d)*K(d,c);
            end
        end
        outfile=strcat(rootdir,'\',fname,'Ambix16.wav');
        audiowrite(outfile,z,fs,'BitsPerSample',32);
        disp(['Saved a 16ch WAV file named ', outfile]);
    end
end
