%----------- B-format dekoder (Semestralni prace) -------------------------
%............3D regular polygonal periphonic layout loudspeaker feeds......
% 16.11.2008
% verze 1 
%--------------------------------------------------------------------------
function [varargout]=regPol3D(W,X,Y,Z,fc,eps)

    n = nargout;
    a = (2*pi)/(n/2);
    
% shelf filtry   
  
    m = tan(pi*fc/48000);
    
% high-frequency boost k1
    
    G1 = sqrt(2);
    k1 = 10^(G1/20);
    jmenovatel1 = 1 + sqrt(2)*m + m^2;
    b01 = (k1 + sqrt(2*k1)*m + m^2) / jmenovatel1;
    b11 = (2*(m^2 - k1)) / jmenovatel1;
    b21 = (k1 - sqrt(2*k1)*m + m^2) / jmenovatel1;
    a11 = (2*(m^2 - 1)) / jmenovatel1;
    a21 = (1 - sqrt(2)*m + m^2) / jmenovatel1;
    
    b1 = [b01 b11 b21];
    a1 = [1 a11 a21];
    
% low-frequency boost k2
    
    G2 = sqrt(3);
    k2 = 10^(G2/20);
    jmenovatel2 = 1 + sqrt(2)*m + m^2;
    b021 = (1 + sqrt(2*k2)*m + k2*m^2) / jmenovatel2;
    b121 = (2*(k2*m^2 - 1)) / jmenovatel2;
    b221 = (1 - sqrt(2*k2)*m + k2*m^2) / jmenovatel2;
    a121 = (2*(m^2 - 1)) / jmenovatel2;
    a221 = (1 - sqrt(2)*m + m^2) / jmenovatel2;
    
    b21 = [b021 b121 b221];
    a21 = [1 a121 a221];
    
% high-frequency boost k2
    
    G2 = sqrt(2);
    k2 = 10^(G2/20);
    jmenovatel2 = 1 + sqrt(2)*m + m^2;
    b022 = (k2 + sqrt(2*k2)*m + m^2) / jmenovatel2;
    b122 = (2*(m^2 - k2)) / jmenovatel2;
    b222 = (k2 - sqrt(2*k2)*m + m^2) / jmenovatel2;
    a122 = (2*(m^2 - 1)) / jmenovatel2;
    a222 = (1 - sqrt(2)*m + m^2) / jmenovatel2;
    
    b22 = [b022 b122 b222];
    a22 = [1 a122 a222];    

% vysledek

    WF = filter(b1,a1, W);
    
    XF1 = filter(b21,a21, X);
    YF1 = filter(b21,a21, Y);
    ZF1 = filter(b21,a21, Z);
    
    XF = filter(b22,a22, XF1);
    YF = filter(b22,a22, YF1);
    ZF = filter(b22,a22, ZF1);    
    

%  loudspeaker at 0 rad
%           for i=0:(n-1)
%               azim1(i+1) = i*a;
%           end
        
% loudspeaker at 0 + a/2 rad
            
            b = a/2;
            for i=0:((n/2) - 1)
                azim(i+1,1) = b + (i*a);
            end

% loudspeaker feeds

%       down loudspeakers
            
            for k=1:(n/2)
                varargout{k}= 0.5*((sqrt(2)*WF) + (cos(azim(k,1))*cos((-eps/2)*(pi/180))*XF) + (sin(azim(k,1))*cos((-eps/2)*(pi/180))*YF) + (sin((-eps/2)*(pi/180))*ZF));
            end

%       up loudspeakers      
           
            for k=((n/2)+1):n
                varargout{k}= 0.5*((sqrt(2)*WF) + (cos(azim(k-(n/2),1))*cos((eps/2)*(pi/180))*XF) + (sin(azim(k-(n/2),1))*cos((eps/2)*(pi/180))*YF) + (sin((eps/2)*(pi/180))*ZF));
            end

disp('Reproduktory jsou číslovány od předního levého proti směru hod. ručiček, ')
disp('přičemž polovina jich je umístěna pod úhlem -eps/2 a polovina pod úhlem +eps/2 od posluchače.')
disp('(eps je elevační úhel mezi horním a spodním reproduktorovým okruhem)')

% zapis na disk
% 
% mkdir('c:\Dekoder-vystupy','regPol3D');
% 
% %   down loudspeakers
%     
%     for i=1:(n/2)
%         speakerFeed=varargout{i};
%         cislo=num2str(i);
%         text=['c:\Dekoder-vystupy\regPol3D\Dspeaker',cislo];
%         wavwrite(0.999*(speakerFeed/max(abs(speakerFeed))),48000,text);
%     end
%     
% %   up loudspeakers    
%     
%     for i=((n/2)+1):n
%         speakerFeed=varargout{i};
%         cislo=num2str(i-(n/2));
%         text=['c:\Dekoder-vystupy\regPol3D\Uspeaker',cislo];
%         wavwrite(0.999*(speakerFeed/max(abs(speakerFeed))),48000,text);
%     end
    
