%----------- B-format dekoder (Semestralni prace) -------------------------
%............BHJ......................
% 9.11.2008
% verze 3 
%--------------------------------------------------------------------------
function [WF,XF,YF] = BHJ(SUMA,DELTA,fc)
    
disp('Konverze B-formatu do formatu BHJ')

    WE= 0.982*SUMA + i*0.164*DELTA;
    XE= 0.419*SUMA - i*0.828*DELTA;
    YE= i*0.385*SUMA + 0.763*DELTA;
    BE= -i*0.694*SUMA + 0.116*DELTA;
    
k = 0.6;
    
% shelf filtry   

    m = tan(pi*fc/48000);
    
% low-frequency cut k1    
    
    G1 = 20*log10(0.646);
    k1 = 10^(-G1/20);
    jmenovatel1 = 1 + sqrt(2*k1)*m + k1*m^2;
    b01 = (1 + sqrt(2)*m + m^2) / jmenovatel1;
    b11 = (2*(m^2 - 1)) / jmenovatel1;
    b21 = (1 - sqrt(2)*m + m^2) / jmenovatel1;
    a11 = (2*(k1*m^2 - 1)) / jmenovatel1;
    a21 = (1 - sqrt(2*k1)*m + k1*m^2) / jmenovatel1;
    
    b1 = [b01 b11 b21];
    a1 = [1 a11 a21];
    
 
% low-frequency boost k2
    
    G2 = 20*log10(1.263);
    k2 = 10^(G2/20);
    jmenovatel2 = 1 + sqrt(2)*m + m^2;
    b02 = (1 + sqrt(2*k2)*m + k2*m^2) / jmenovatel2;
    b12 = (2*(k2*m^2 - 1)) / jmenovatel2;
    b22 = (1 - sqrt(2*k2)*m + k2*m^2) / jmenovatel2;
    a12 = (2*(m^2 - 1)) / jmenovatel2;
    a22 = (1 - sqrt(2)*m + m^2) / jmenovatel2;
    
    b2 = [b02 b12 b22];
    a2 = [1 a12 a22];
    
    
    k2YE = filter(b2,a2, YE);
    
% low-frequency cut k3    
    
    G3 = 20*log10(0.775);
    k3 = 10^(-G3/20);
    jmenovatel3 = 1 + sqrt(2*k3)*m + k3*m^2;
    b03 = (1 + sqrt(2)*m + m^2) / jmenovatel3;
    b13 = (2*(m^2 - 1)) / jmenovatel3;
    b23 = (1 - sqrt(2)*m + m^2) / jmenovatel3;
    a13 = (2*(k3*m^2 - 1)) / jmenovatel3;
    a23 = (1 - sqrt(2*k3)*m + k3*m^2) / jmenovatel3;
    
    b3 = [b03 b13 b23];
    a3 = [1 a13 a23];
    
    BF = filter(b3,a3, BE);
    kk3BF = k*BF;

% vysledek

    WF = filter(b1,a1, WE);
    XF = filter(b2,a2, XE);
    YF = k2YE + kk3BF;
    
    
% zobrazeni kmitoctove charakteristiky
% figure(1);
% 
% [H1,w1] = freqz(b1, a1, 1024, 48000);
% subplot(3,1,1);
% semilogx (w1, 20*log10(abs(H1)));
% handle = title('Modulová kmitočtová charakteristika');
% set(handle, 'FontName', 'Times New Roman CE');
% xlabel('{\it f} [Hz] \rightarrow');
% ylabel('|{\it H }| [dB] \rightarrow');
% grid on;
% 
% [H2,w2] = freqz(b2, a2, 1024, 48000);
% subplot(3,1,2);
% semilogx (w2, 20*log10(abs(H2)));
% handle = title('Modulová kmitočtová charakteristika');
% set(handle, 'FontName', 'Times New Roman CE');
% xlabel('{\it f} [Hz] \rightarrow');
% ylabel('|{\it H }| [dB] \rightarrow');
% grid on;
% 
% [H3,w3] = freqz(b3, a3, 1024, 48000);
% subplot(3,1,3);
% semilogx (w3, 20*log10(abs(H3)));
% handle = title('Modulová kmitočtová charakteristika');
% set(handle, 'FontName', 'Times New Roman CE');
% xlabel('{\it f} [Hz] \rightarrow');
% ylabel('|{\it H }| [dB] \rightarrow');
% grid on;
    
% .......................faze.......................

% figure(2)
% 
% subplot(3,1,1);
% semilogx (w1, 180/pi*unwrap(angle(H1)));
% handle = title('Fázová kmitočtová charakteristika');
% set(handle, 'FontName', 'Times New Roman CE');
% xlabel('{\it f} [Hz] \rightarrow');
% ylabel('angle({\itH}) [\circ] \rightarrow');
% grid on;
% 
% subplot(3,1,2);
% semilogx (w2, 180/pi*unwrap(angle(H2)));
% handle = title('Fázová kmitočtová charakteristika');
% set(handle, 'FontName', 'Times New Roman CE');
% xlabel('{\it f} [Hz] \rightarrow');
% ylabel('angle({\itH}) [\circ] \rightarrow');
% grid on;
% 
% subplot(3,1,3);
% semilogx (w3, 180/pi*unwrap(angle(H3)));
% handle = title('Fázová kmitočtová charakteristika');
% set(handle, 'FontName', 'Times New Roman CE');
% xlabel('{\it f} [Hz] \rightarrow');
% ylabel('angle({\itH}) [\circ] \rightarrow');
% grid on;
